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