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