1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 2 /* matimpl.h is needed only for logging of matrix operation */ 3 #include <petsc/private/matimpl.h> 4 5 PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES); 6 7 PETSC_INTERN PetscErrorCode SNESDiffParameterCreate_More(SNES, Vec, void **); 8 PETSC_INTERN PetscErrorCode SNESDiffParameterCompute_More(SNES, void *, Vec, Vec, PetscReal *, PetscReal *); 9 PETSC_INTERN PetscErrorCode SNESDiffParameterDestroy_More(void *); 10 11 typedef struct { /* default context for matrix-free SNES */ 12 SNES snes; /* SNES context */ 13 Vec w; /* work vector */ 14 MatNullSpace sp; /* null space context */ 15 PetscReal error_rel; /* square root of relative error in computing function */ 16 PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */ 17 PetscBool jorge; /* flag indicating use of Jorge's method for determining the differencing parameter */ 18 PetscReal h; /* differencing parameter */ 19 PetscBool need_h; /* flag indicating whether we must compute h */ 20 PetscBool need_err; /* flag indicating whether we must currently compute error_rel */ 21 PetscBool compute_err; /* flag indicating whether we must ever compute error_rel */ 22 PetscInt compute_err_iter; /* last iter where we've computer error_rel */ 23 PetscInt compute_err_freq; /* frequency of computing error_rel */ 24 void *data; /* implementation-specific data */ 25 } MFCtx_Private; 26 27 static PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat) 28 { 29 MFCtx_Private *ctx; 30 31 PetscFunctionBegin; 32 PetscCall(MatShellGetContext(mat, &ctx)); 33 PetscCall(VecDestroy(&ctx->w)); 34 PetscCall(MatNullSpaceDestroy(&ctx->sp)); 35 if (ctx->jorge || ctx->compute_err) PetscCall(SNESDiffParameterDestroy_More(ctx->data)); 36 PetscCall(PetscFree(ctx)); 37 PetscFunctionReturn(PETSC_SUCCESS); 38 } 39 40 /* 41 SNESMatrixFreeView2_Private - Views matrix-free parameters. 42 */ 43 static PetscErrorCode SNESMatrixFreeView2_Private(Mat J, PetscViewer viewer) 44 { 45 MFCtx_Private *ctx; 46 PetscBool isascii; 47 48 PetscFunctionBegin; 49 PetscCall(MatShellGetContext(J, &ctx)); 50 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 51 if (isascii) { 52 PetscCall(PetscViewerASCIIPrintf(viewer, " SNES matrix-free approximation:\n")); 53 if (ctx->jorge) PetscCall(PetscViewerASCIIPrintf(viewer, " using Jorge's method of determining differencing parameter\n")); 54 PetscCall(PetscViewerASCIIPrintf(viewer, " err=%g (relative error in function evaluation)\n", (double)ctx->error_rel)); 55 PetscCall(PetscViewerASCIIPrintf(viewer, " umin=%g (minimum iterate parameter)\n", (double)ctx->umin)); 56 if (ctx->compute_err) PetscCall(PetscViewerASCIIPrintf(viewer, " freq_err=%" PetscInt_FMT " (frequency for computing err)\n", ctx->compute_err_freq)); 57 } 58 PetscFunctionReturn(PETSC_SUCCESS); 59 } 60 61 /* 62 SNESMatrixFreeMult2_Private - Default matrix-free form for Jacobian-vector 63 product, y = F'(u)*a: 64 y = (F(u + ha) - F(u)) /h, 65 where F = nonlinear function, as set by SNESSetFunction() 66 u = current iterate 67 h = difference interval 68 */ 69 static PetscErrorCode SNESMatrixFreeMult2_Private(Mat mat, Vec a, Vec y) 70 { 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(PETSC_SUCCESS); 148 } 149 150 /*@ 151 MatCreateSNESMFMore - Creates a matrix-free matrix 152 context for use with a `SNES` solver that uses the More method to compute an optimal h based on the noise of the function. 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 Options Database Keys: 163 + -snes_mf_err <error_rel> - see `MatCreateSNESMF()` 164 . -snes_mf_umin <umin> - see `MatCreateSNESMF()` 165 . -snes_mf_compute_err - compute the square root or relative error in function 166 . -snes_mf_freq_err <freq> - set the frequency to recompute the parameters 167 - -snes_mf_jorge - use the method of Jorge More 168 169 Level: advanced 170 171 Notes: 172 This is an experimental approach, use `MatCreateSNESMF()`. 173 174 The matrix-free matrix context merely contains the function pointers 175 and work space for performing finite difference approximations of 176 Jacobian-vector products, J(u)*a, via 177 178 .vb 179 J(u)*a = [J(u+h*a) - J(u)]/h, 180 where by default 181 h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 182 = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 183 where 184 error_rel = square root of relative error in function evaluation 185 umin = minimum iterate parameter 186 Alternatively, the differencing parameter, h, can be set using 187 Jorge's nifty new strategy if one specifies the option 188 -snes_mf_jorge 189 .ve 190 191 The user can set these parameters via `MatMFFDSetFunctionError()`. 192 193 The user should call `MatDestroy()` when finished with the matrix-free 194 matrix context. 195 196 .seealso: [](ch_snes), `SNESCreateMF()`, `MatCreateMFFD()`, `MatDestroy()`, `MatMFFDSetFunctionError()` 197 @*/ 198 PetscErrorCode MatCreateSNESMFMore(SNES snes, Vec x, Mat *J) 199 { 200 MPI_Comm comm; 201 MFCtx_Private *mfctx; 202 PetscInt n, nloc; 203 PetscBool flg; 204 char p[64]; 205 206 PetscFunctionBegin; 207 PetscCall(PetscNew(&mfctx)); 208 mfctx->sp = NULL; 209 mfctx->snes = snes; 210 mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON; 211 mfctx->umin = 1.e-6; 212 mfctx->h = 0.0; 213 mfctx->need_h = PETSC_TRUE; 214 mfctx->need_err = PETSC_FALSE; 215 mfctx->compute_err = PETSC_FALSE; 216 mfctx->compute_err_freq = 0; 217 mfctx->compute_err_iter = -1; 218 mfctx->compute_err = PETSC_FALSE; 219 mfctx->jorge = PETSC_FALSE; 220 221 PetscCall(PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_err", &mfctx->error_rel, NULL)); 222 PetscCall(PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_umin", &mfctx->umin, NULL)); 223 PetscCall(PetscOptionsGetBool(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_jorge", &mfctx->jorge, NULL)); 224 PetscCall(PetscOptionsGetBool(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_compute_err", &mfctx->compute_err, NULL)); 225 PetscCall(PetscOptionsGetInt(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_freq_err", &mfctx->compute_err_freq, &flg)); 226 if (flg) { 227 if (mfctx->compute_err_freq < 0) mfctx->compute_err_freq = 0; 228 mfctx->compute_err = PETSC_TRUE; 229 } 230 if (mfctx->compute_err) mfctx->need_err = PETSC_TRUE; 231 if (mfctx->jorge || mfctx->compute_err) PetscCall(SNESDiffParameterCreate_More(snes, x, &mfctx->data)); 232 else mfctx->data = NULL; 233 234 PetscCall(PetscOptionsHasHelp(((PetscObject)snes)->options, &flg)); 235 PetscCall(PetscStrncpy(p, "-", sizeof(p))); 236 if (((PetscObject)snes)->prefix) PetscCall(PetscStrlcat(p, ((PetscObject)snes)->prefix, sizeof(p))); 237 if (flg) { 238 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " Matrix-free Options (via SNES):\n")); 239 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_err <err>: set sqrt of relative error in function (default %g)\n", p, (double)mfctx->error_rel)); 240 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_umin <umin>: see users manual (default %g)\n", p, (double)mfctx->umin)); 241 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_jorge: use Jorge More's method\n", p)); 242 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_compute_err: compute sqrt or relative error in function\n", p)); 243 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n", p)); 244 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_noise_file <file>: set file for printing noise info\n", p)); 245 } 246 PetscCall(VecDuplicate(x, &mfctx->w)); 247 PetscCall(PetscObjectGetComm((PetscObject)x, &comm)); 248 PetscCall(VecGetSize(x, &n)); 249 PetscCall(VecGetLocalSize(x, &nloc)); 250 PetscCall(MatCreate(comm, J)); 251 PetscCall(MatSetSizes(*J, nloc, n, n, n)); 252 PetscCall(MatSetType(*J, MATSHELL)); 253 PetscCall(MatShellSetContext(*J, mfctx)); 254 PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))SNESMatrixFreeMult2_Private)); 255 PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))SNESMatrixFreeDestroy2_Private)); 256 PetscCall(MatShellSetOperation(*J, MATOP_VIEW, (void (*)(void))SNESMatrixFreeView2_Private)); 257 PetscCall(MatSetUp(*J)); 258 PetscFunctionReturn(PETSC_SUCCESS); 259 } 260 261 /*@ 262 MatSNESMFMoreSetParameters - Sets the parameters for the approximation of 263 matrix-vector products using finite differences, see `MatCreateSNESMFMore()` 264 265 Input Parameters: 266 + mat - the matrix 267 . error - relative error (should be set to the square root of the relative error in the function evaluations) 268 . umin - minimum allowable u-value 269 - h - differencing parameter 270 271 Options Database Keys: 272 + -snes_mf_err <error_rel> - see `MatCreateSNESMF()` 273 . -snes_mf_umin <umin> - see `MatCreateSNESMF()` 274 . -snes_mf_compute_err - compute the square root or relative error in function 275 . -snes_mf_freq_err <freq> - set the frequency to recompute the parameters 276 - -snes_mf_jorge - use the method of Jorge More 277 278 Level: advanced 279 280 Note: 281 If the user sets the parameter `h` directly, then this value will be used 282 instead of the default computation as discussed in `MatCreateSNESMFMore()` 283 284 .seealso: [](ch_snes), `SNES`, `MatCreateSNESMF()`, `MatCreateSNESMFMore()` 285 @*/ 286 PetscErrorCode MatSNESMFMoreSetParameters(Mat mat, PetscReal error, PetscReal umin, PetscReal h) 287 { 288 MFCtx_Private *ctx; 289 290 PetscFunctionBegin; 291 PetscCall(MatShellGetContext(mat, &ctx)); 292 if (ctx) { 293 if (error != (PetscReal)PETSC_DEFAULT) ctx->error_rel = error; 294 if (umin != (PetscReal)PETSC_DEFAULT) ctx->umin = umin; 295 if (h != (PetscReal)PETSC_DEFAULT) { 296 ctx->h = h; 297 ctx->need_h = PETSC_FALSE; 298 } 299 } 300 PetscFunctionReturn(PETSC_SUCCESS); 301 } 302 303 PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes) 304 { 305 MFCtx_Private *ctx; 306 Mat mat; 307 308 PetscFunctionBegin; 309 PetscCall(SNESGetJacobian(snes, &mat, NULL, NULL, NULL)); 310 PetscCall(MatShellGetContext(mat, &ctx)); 311 if (ctx) ctx->need_h = PETSC_TRUE; 312 PetscFunctionReturn(PETSC_SUCCESS); 313 } 314