1 2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 3 /* matimpl.h is needed only for logging of matrix operation */ 4 #include <petsc/private/matimpl.h> 5 6 extern PetscErrorCode SNESDiffParameterCreate_More(SNES,Vec,void**); 7 extern PetscErrorCode SNESDiffParameterCompute_More(SNES,void*,Vec,Vec,PetscReal*,PetscReal*); 8 extern PetscErrorCode SNESDiffParameterDestroy_More(void*); 9 10 typedef struct { /* default context for matrix-free SNES */ 11 SNES snes; /* SNES context */ 12 Vec w; /* work vector */ 13 MatNullSpace sp; /* null space context */ 14 PetscReal error_rel; /* square root of relative error in computing function */ 15 PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */ 16 PetscBool jorge; /* flag indicating use of Jorge's method for determining the differencing parameter */ 17 PetscReal h; /* differencing parameter */ 18 PetscBool need_h; /* flag indicating whether we must compute h */ 19 PetscBool need_err; /* flag indicating whether we must currently compute error_rel */ 20 PetscBool compute_err; /* flag indicating whether we must ever compute error_rel */ 21 PetscInt compute_err_iter; /* last iter where we've computer error_rel */ 22 PetscInt compute_err_freq; /* frequency of computing error_rel */ 23 void *data; /* implementation-specific data */ 24 } MFCtx_Private; 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "SNESMatrixFreeDestroy2_Private" 28 PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat) 29 { 30 PetscErrorCode ierr; 31 MFCtx_Private *ctx; 32 33 PetscFunctionBegin; 34 ierr = MatShellGetContext(mat,(void**)&ctx);CHKERRQ(ierr); 35 ierr = VecDestroy(&ctx->w);CHKERRQ(ierr); 36 ierr = MatNullSpaceDestroy(&ctx->sp);CHKERRQ(ierr); 37 if (ctx->jorge || ctx->compute_err) {ierr = SNESDiffParameterDestroy_More(ctx->data);CHKERRQ(ierr);} 38 ierr = PetscFree(ctx);CHKERRQ(ierr); 39 PetscFunctionReturn(0); 40 } 41 42 #undef __FUNCT__ 43 #define __FUNCT__ "SNESMatrixFreeView2_Private" 44 /* 45 SNESMatrixFreeView2_Private - Views matrix-free parameters. 46 */ 47 PetscErrorCode SNESMatrixFreeView2_Private(Mat J,PetscViewer viewer) 48 { 49 PetscErrorCode ierr; 50 MFCtx_Private *ctx; 51 PetscBool iascii; 52 53 PetscFunctionBegin; 54 ierr = MatShellGetContext(J,(void**)&ctx);CHKERRQ(ierr); 55 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 56 if (iascii) { 57 ierr = PetscViewerASCIIPrintf(viewer," SNES matrix-free approximation:\n");CHKERRQ(ierr); 58 if (ctx->jorge) { 59 ierr = PetscViewerASCIIPrintf(viewer," using Jorge's method of determining differencing parameter\n");CHKERRQ(ierr); 60 } 61 ierr = PetscViewerASCIIPrintf(viewer," err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);CHKERRQ(ierr); 62 ierr = PetscViewerASCIIPrintf(viewer," umin=%g (minimum iterate parameter)\n",(double)ctx->umin);CHKERRQ(ierr); 63 if (ctx->compute_err) { 64 ierr = PetscViewerASCIIPrintf(viewer," freq_err=%D (frequency for computing err)\n",ctx->compute_err_freq);CHKERRQ(ierr); 65 } 66 } 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "SNESMatrixFreeMult2_Private" 72 /* 73 SNESMatrixFreeMult2_Private - Default matrix-free form for Jacobian-vector 74 product, y = F'(u)*a: 75 y = (F(u + ha) - F(u)) /h, 76 where F = nonlinear function, as set by SNESSetFunction() 77 u = current iterate 78 h = difference interval 79 */ 80 PetscErrorCode SNESMatrixFreeMult2_Private(Mat mat,Vec a,Vec y) 81 { 82 MFCtx_Private *ctx; 83 SNES snes; 84 PetscReal h,norm,sum,umin,noise; 85 PetscScalar hs,dot; 86 Vec w,U,F; 87 PetscErrorCode ierr,(*eval_fct)(SNES,Vec,Vec); 88 MPI_Comm comm; 89 PetscInt iter; 90 91 PetscFunctionBegin; 92 /* We log matrix-free matrix-vector products separately, so that we can 93 separate the performance monitoring from the cases that use conventional 94 storage. We may eventually modify event logging to associate events 95 with particular objects, hence alleviating the more general problem. */ 96 ierr = PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);CHKERRQ(ierr); 97 98 ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr); 99 ierr = MatShellGetContext(mat,(void**)&ctx);CHKERRQ(ierr); 100 snes = ctx->snes; 101 w = ctx->w; 102 umin = ctx->umin; 103 104 ierr = SNESGetSolution(snes,&U);CHKERRQ(ierr); 105 eval_fct = SNESComputeFunction; 106 ierr = SNESGetFunction(snes,&F,NULL,NULL);CHKERRQ(ierr); 107 108 /* Determine a "good" step size, h */ 109 if (ctx->need_h) { 110 111 /* Use Jorge's method to compute h */ 112 if (ctx->jorge) { 113 ierr = SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);CHKERRQ(ierr); 114 115 /* Use the Brown/Saad method to compute h */ 116 } else { 117 /* Compute error if desired */ 118 ierr = SNESGetIterationNumber(snes,&iter);CHKERRQ(ierr); 119 if ((ctx->need_err) || ((ctx->compute_err_freq) && (ctx->compute_err_iter != iter) && (!((iter-1)%ctx->compute_err_freq)))) { 120 /* Use Jorge's method to compute noise */ 121 ierr = SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);CHKERRQ(ierr); 122 123 ctx->error_rel = PetscSqrtReal(noise); 124 125 ierr = PetscInfo3(snes,"Using Jorge's noise: noise=%g, sqrt(noise)=%g, h_more=%g\n",(double)noise,(double)ctx->error_rel,(double)h);CHKERRQ(ierr); 126 127 ctx->compute_err_iter = iter; 128 ctx->need_err = PETSC_FALSE; 129 } 130 131 ierr = VecDotBegin(U,a,&dot);CHKERRQ(ierr); 132 ierr = VecNormBegin(a,NORM_1,&sum);CHKERRQ(ierr); 133 ierr = VecNormBegin(a,NORM_2,&norm);CHKERRQ(ierr); 134 ierr = VecDotEnd(U,a,&dot);CHKERRQ(ierr); 135 ierr = VecNormEnd(a,NORM_1,&sum);CHKERRQ(ierr); 136 ierr = VecNormEnd(a,NORM_2,&norm);CHKERRQ(ierr); 137 138 139 /* Safeguard for step sizes too small */ 140 if (sum == 0.0) { 141 dot = 1.0; 142 norm = 1.0; 143 } else if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum; 144 else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum; 145 h = PetscRealPart(ctx->error_rel*dot/(norm*norm)); 146 } 147 } else h = ctx->h; 148 149 if (!ctx->jorge || !ctx->need_h) {ierr = PetscInfo1(snes,"h = %g\n",(double)h);CHKERRQ(ierr);} 150 151 /* Evaluate function at F(u + ha) */ 152 hs = h; 153 ierr = VecWAXPY(w,hs,a,U);CHKERRQ(ierr); 154 ierr = eval_fct(snes,w,y);CHKERRQ(ierr); 155 ierr = VecAXPY(y,-1.0,F);CHKERRQ(ierr); 156 ierr = VecScale(y,1.0/hs);CHKERRQ(ierr); 157 if (ctx->sp) {ierr = MatNullSpaceRemove(ctx->sp,y);CHKERRQ(ierr);} 158 159 ierr = PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "SNESDefaultMatrixFreeCreate2" 165 /*@C 166 SNESMatrixFreeCreate2 - Creates a matrix-free matrix 167 context for use with a SNES solver. This matrix can be used as 168 the Jacobian argument for the routine SNESSetJacobian(). 169 170 Input Parameters: 171 . snes - the SNES context 172 . x - vector where SNES solution is to be stored. 173 174 Output Parameter: 175 . J - the matrix-free matrix 176 177 Level: advanced 178 179 Notes: 180 The matrix-free matrix context merely contains the function pointers 181 and work space for performing finite difference approximations of 182 Jacobian-vector products, J(u)*a, via 183 184 $ J(u)*a = [J(u+h*a) - J(u)]/h, 185 $ where by default 186 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 187 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 188 $ where 189 $ error_rel = square root of relative error in 190 $ function evaluation 191 $ umin = minimum iterate parameter 192 $ Alternatively, the differencing parameter, h, can be set using 193 $ Jorge's nifty new strategy if one specifies the option 194 $ -snes_mf_jorge 195 196 The user can set these parameters via MatMFFDSetFunctionError(). 197 See Users-Manual: ch_snes for details 198 199 The user should call MatDestroy() when finished with the matrix-free 200 matrix context. 201 202 Options Database Keys: 203 $ -snes_mf_err <error_rel> 204 $ -snes_mf_unim <umin> 205 $ -snes_mf_compute_err 206 $ -snes_mf_freq_err <freq> 207 $ -snes_mf_jorge 208 209 .keywords: SNES, default, matrix-free, create, matrix 210 211 .seealso: MatDestroy(), MatMFFDSetFunctionError() 212 @*/ 213 PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES snes,Vec x,Mat *J) 214 { 215 MPI_Comm comm; 216 MFCtx_Private *mfctx; 217 PetscErrorCode ierr; 218 PetscInt n,nloc; 219 PetscBool flg; 220 char p[64]; 221 222 PetscFunctionBegin; 223 ierr = PetscNewLog(snes,&mfctx);CHKERRQ(ierr); 224 mfctx->sp = 0; 225 mfctx->snes = snes; 226 mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON; 227 mfctx->umin = 1.e-6; 228 mfctx->h = 0.0; 229 mfctx->need_h = PETSC_TRUE; 230 mfctx->need_err = PETSC_FALSE; 231 mfctx->compute_err = PETSC_FALSE; 232 mfctx->compute_err_freq = 0; 233 mfctx->compute_err_iter = -1; 234 mfctx->compute_err = PETSC_FALSE; 235 mfctx->jorge = PETSC_FALSE; 236 237 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_mf_err",&mfctx->error_rel,NULL);CHKERRQ(ierr); 238 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_mf_umin",&mfctx->umin,NULL);CHKERRQ(ierr); 239 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_mf_jorge",&mfctx->jorge,NULL);CHKERRQ(ierr); 240 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_mf_compute_err",&mfctx->compute_err,NULL);CHKERRQ(ierr); 241 ierr = PetscOptionsGetInt(((PetscObject)snes)->prefix,"-snes_mf_freq_err",&mfctx->compute_err_freq,&flg);CHKERRQ(ierr); 242 if (flg) { 243 if (mfctx->compute_err_freq < 0) mfctx->compute_err_freq = 0; 244 mfctx->compute_err = PETSC_TRUE; 245 } 246 if (mfctx->compute_err) mfctx->need_err = PETSC_TRUE; 247 if (mfctx->jorge || mfctx->compute_err) { 248 ierr = SNESDiffParameterCreate_More(snes,x,&mfctx->data);CHKERRQ(ierr); 249 } else mfctx->data = 0; 250 251 ierr = PetscOptionsHasName(NULL,"-help",&flg);CHKERRQ(ierr); 252 ierr = PetscStrcpy(p,"-");CHKERRQ(ierr); 253 if (((PetscObject)snes)->prefix) PetscStrcat(p,((PetscObject)snes)->prefix); 254 if (flg) { 255 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," Matrix-free Options (via SNES):\n");CHKERRQ(ierr); 256 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_err <err>: set sqrt of relative error in function (default %g)\n",p,(double)mfctx->error_rel);CHKERRQ(ierr); 257 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_umin <umin>: see users manual (default %g)\n",p,(double)mfctx->umin);CHKERRQ(ierr); 258 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_jorge: use Jorge More's method\n",p);CHKERRQ(ierr); 259 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p);CHKERRQ(ierr); 260 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p);CHKERRQ(ierr); 261 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_noise_file <file>: set file for printing noise info\n",p);CHKERRQ(ierr); 262 } 263 ierr = VecDuplicate(x,&mfctx->w);CHKERRQ(ierr); 264 ierr = PetscObjectGetComm((PetscObject)x,&comm);CHKERRQ(ierr); 265 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 266 ierr = VecGetLocalSize(x,&nloc);CHKERRQ(ierr); 267 ierr = MatCreate(comm,J);CHKERRQ(ierr); 268 ierr = MatSetSizes(*J,nloc,n,n,n);CHKERRQ(ierr); 269 ierr = MatSetType(*J,MATSHELL);CHKERRQ(ierr); 270 ierr = MatShellSetContext(*J,mfctx);CHKERRQ(ierr); 271 ierr = MatShellSetOperation(*J,MATOP_MULT,(void (*)(void))SNESMatrixFreeMult2_Private);CHKERRQ(ierr); 272 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void (*)(void))SNESMatrixFreeDestroy2_Private);CHKERRQ(ierr); 273 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void (*)(void))SNESMatrixFreeView2_Private);CHKERRQ(ierr); 274 275 ierr = PetscLogObjectParent((PetscObject)*J,(PetscObject)mfctx->w);CHKERRQ(ierr); 276 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)*J);CHKERRQ(ierr); 277 PetscFunctionReturn(0); 278 } 279 280 #undef __FUNCT__ 281 #define __FUNCT__ "SNESDefaultMatrixFreeSetParameters2" 282 /*@C 283 SNESDefaultMatrixFreeSetParameters2 - Sets the parameters for the approximation of 284 matrix-vector products using finite differences. 285 286 $ J(u)*a = [J(u+h*a) - J(u)]/h where 287 288 either the user sets h directly here, or this parameter is computed via 289 290 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 291 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 292 $ 293 294 Input Parameters: 295 + mat - the matrix 296 . error_rel - relative error (should be set to the square root of 297 the relative error in the function evaluations) 298 . umin - minimum allowable u-value 299 - h - differencing parameter 300 301 Level: advanced 302 303 Notes: 304 If the user sets the parameter h directly, then this value will be used 305 instead of the default computation indicated above. 306 307 .keywords: SNES, matrix-free, parameters 308 309 .seealso: MatCreateSNESMF() 310 @*/ 311 PetscErrorCode SNESDefaultMatrixFreeSetParameters2(Mat mat,PetscReal error,PetscReal umin,PetscReal h) 312 { 313 MFCtx_Private *ctx; 314 PetscErrorCode ierr; 315 316 PetscFunctionBegin; 317 ierr = MatShellGetContext(mat,(void**)&ctx);CHKERRQ(ierr); 318 if (ctx) { 319 if (error != PETSC_DEFAULT) ctx->error_rel = error; 320 if (umin != PETSC_DEFAULT) ctx->umin = umin; 321 if (h != PETSC_DEFAULT) { 322 ctx->h = h; 323 ctx->need_h = PETSC_FALSE; 324 } 325 } 326 PetscFunctionReturn(0); 327 } 328 329 PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes) 330 { 331 MFCtx_Private *ctx; 332 PetscErrorCode ierr; 333 Mat mat; 334 335 PetscFunctionBegin; 336 ierr = SNESGetJacobian(snes,&mat,NULL,NULL,NULL);CHKERRQ(ierr); 337 ierr = MatShellGetContext(mat,(void**)&ctx);CHKERRQ(ierr); 338 if (ctx) ctx->need_h = PETSC_TRUE; 339 PetscFunctionReturn(0); 340 } 341 342