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 8 PETSC_INTERN PetscErrorCode SNESDiffParameterCreate_More(SNES, Vec, void **); 9 PETSC_INTERN PetscErrorCode SNESDiffParameterCompute_More(SNES, void *, Vec, Vec, PetscReal *, PetscReal *); 10 PETSC_INTERN PetscErrorCode SNESDiffParameterDestroy_More(void *); 11 12 typedef struct { /* default context for matrix-free SNES */ 13 SNES snes; /* SNES context */ 14 Vec w; /* work vector */ 15 MatNullSpace sp; /* null space context */ 16 PetscReal error_rel; /* square root of relative error in computing function */ 17 PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */ 18 PetscBool jorge; /* flag indicating use of Jorge's method for determining the differencing parameter */ 19 PetscReal h; /* differencing parameter */ 20 PetscBool need_h; /* flag indicating whether we must compute h */ 21 PetscBool need_err; /* flag indicating whether we must currently compute error_rel */ 22 PetscBool compute_err; /* flag indicating whether we must ever compute error_rel */ 23 PetscInt compute_err_iter; /* last iter where we've computer error_rel */ 24 PetscInt compute_err_freq; /* frequency of computing error_rel */ 25 void *data; /* implementation-specific data */ 26 } MFCtx_Private; 27 28 PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat) { 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(0); 38 } 39 40 /* 41 SNESMatrixFreeView2_Private - Views matrix-free parameters. 42 */ 43 PetscErrorCode SNESMatrixFreeView2_Private(Mat J, PetscViewer viewer) { 44 MFCtx_Private *ctx; 45 PetscBool iascii; 46 47 PetscFunctionBegin; 48 PetscCall(MatShellGetContext(J, &ctx)); 49 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 50 if (iascii) { 51 PetscCall(PetscViewerASCIIPrintf(viewer, " SNES matrix-free approximation:\n")); 52 if (ctx->jorge) PetscCall(PetscViewerASCIIPrintf(viewer, " using Jorge's method of determining differencing parameter\n")); 53 PetscCall(PetscViewerASCIIPrintf(viewer, " err=%g (relative error in function evaluation)\n", (double)ctx->error_rel)); 54 PetscCall(PetscViewerASCIIPrintf(viewer, " umin=%g (minimum iterate parameter)\n", (double)ctx->umin)); 55 if (ctx->compute_err) PetscCall(PetscViewerASCIIPrintf(viewer, " freq_err=%" PetscInt_FMT " (frequency for computing err)\n", ctx->compute_err_freq)); 56 } 57 PetscFunctionReturn(0); 58 } 59 60 /* 61 SNESMatrixFreeMult2_Private - Default matrix-free form for Jacobian-vector 62 product, y = F'(u)*a: 63 y = (F(u + ha) - F(u)) /h, 64 where F = nonlinear function, as set by SNESSetFunction() 65 u = current iterate 66 h = difference interval 67 */ 68 PetscErrorCode SNESMatrixFreeMult2_Private(Mat mat, Vec a, Vec y) { 69 MFCtx_Private *ctx; 70 SNES snes; 71 PetscReal h, norm, sum, umin, noise; 72 PetscScalar hs, dot; 73 Vec w, U, F; 74 MPI_Comm comm; 75 PetscInt iter; 76 PetscErrorCode (*eval_fct)(SNES, Vec, Vec); 77 78 PetscFunctionBegin; 79 /* We log matrix-free matrix-vector products separately, so that we can 80 separate the performance monitoring from the cases that use conventional 81 storage. We may eventually modify event logging to associate events 82 with particular objects, hence alleviating the more general problem. */ 83 PetscCall(PetscLogEventBegin(MATMFFD_Mult, a, y, 0, 0)); 84 85 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 86 PetscCall(MatShellGetContext(mat, &ctx)); 87 snes = ctx->snes; 88 w = ctx->w; 89 umin = ctx->umin; 90 91 PetscCall(SNESGetSolution(snes, &U)); 92 eval_fct = SNESComputeFunction; 93 PetscCall(SNESGetFunction(snes, &F, NULL, NULL)); 94 95 /* Determine a "good" step size, h */ 96 if (ctx->need_h) { 97 /* Use Jorge's method to compute h */ 98 if (ctx->jorge) { 99 PetscCall(SNESDiffParameterCompute_More(snes, ctx->data, U, a, &noise, &h)); 100 101 /* Use the Brown/Saad method to compute h */ 102 } else { 103 /* Compute error if desired */ 104 PetscCall(SNESGetIterationNumber(snes, &iter)); 105 if ((ctx->need_err) || ((ctx->compute_err_freq) && (ctx->compute_err_iter != iter) && (!((iter - 1) % ctx->compute_err_freq)))) { 106 /* Use Jorge's method to compute noise */ 107 PetscCall(SNESDiffParameterCompute_More(snes, ctx->data, U, a, &noise, &h)); 108 109 ctx->error_rel = PetscSqrtReal(noise); 110 111 PetscCall(PetscInfo(snes, "Using Jorge's noise: noise=%g, sqrt(noise)=%g, h_more=%g\n", (double)noise, (double)ctx->error_rel, (double)h)); 112 113 ctx->compute_err_iter = iter; 114 ctx->need_err = PETSC_FALSE; 115 } 116 117 PetscCall(VecDotBegin(U, a, &dot)); 118 PetscCall(VecNormBegin(a, NORM_1, &sum)); 119 PetscCall(VecNormBegin(a, NORM_2, &norm)); 120 PetscCall(VecDotEnd(U, a, &dot)); 121 PetscCall(VecNormEnd(a, NORM_1, &sum)); 122 PetscCall(VecNormEnd(a, NORM_2, &norm)); 123 124 /* Safeguard for step sizes too small */ 125 if (sum == 0.0) { 126 dot = 1.0; 127 norm = 1.0; 128 } else if (PetscAbsScalar(dot) < umin * sum && PetscRealPart(dot) >= 0.0) dot = umin * sum; 129 else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin * sum) dot = -umin * sum; 130 h = PetscRealPart(ctx->error_rel * dot / (norm * norm)); 131 } 132 } else h = ctx->h; 133 134 if (!ctx->jorge || !ctx->need_h) PetscCall(PetscInfo(snes, "h = %g\n", (double)h)); 135 136 /* Evaluate function at F(u + ha) */ 137 hs = h; 138 PetscCall(VecWAXPY(w, hs, a, U)); 139 PetscCall(eval_fct(snes, w, y)); 140 PetscCall(VecAXPY(y, -1.0, F)); 141 PetscCall(VecScale(y, 1.0 / hs)); 142 if (mat->nullsp) PetscCall(MatNullSpaceRemove(mat->nullsp, y)); 143 144 PetscCall(PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0)); 145 PetscFunctionReturn(0); 146 } 147 148 /*@C 149 MatCreateSNESMFMore - Creates a matrix-free matrix 150 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 151 the Jacobian argument for the routine `SNESSetJacobian()`. 152 153 Input Parameters: 154 + snes - the `SNES` context 155 - x - vector where `SNES` solution is to be stored. 156 157 Output Parameter: 158 . J - the matrix-free matrix 159 160 Options Database Keys: 161 + -snes_mf_err <error_rel> - see `MatCreateSNESMF()` 162 . -snes_mf_umin <umin> - see `MatCreateSNESMF()` 163 . -snes_mf_compute_err - compute the square root or relative error in function 164 . -snes_mf_freq_err <freq> - set the frequency to recompute the parameters 165 - -snes_mf_jorge - use the method of Jorge More 166 167 Level: advanced 168 169 Notes: 170 This is an experimental approach, use `MatCreateSNESMF()`. 171 172 The matrix-free matrix context merely contains the function pointers 173 and work space for performing finite difference approximations of 174 Jacobian-vector products, J(u)*a, via 175 176 $ J(u)*a = [J(u+h*a) - J(u)]/h, 177 $ where by default 178 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 179 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 180 $ where 181 $ error_rel = square root of relative error in 182 $ function evaluation 183 $ umin = minimum iterate parameter 184 $ Alternatively, the differencing parameter, h, can be set using 185 $ Jorge's nifty new strategy if one specifies the option 186 $ -snes_mf_jorge 187 188 The user can set these parameters via `MatMFFDSetFunctionError()`. 189 See Users-Manual: ch_snes for details 190 191 The user should call `MatDestroy()` when finished with the matrix-free 192 matrix context. 193 194 .seealso: `SNESCreateMF`(), `MatCreateMFFD()`, `MatDestroy()`, `MatMFFDSetFunctionError()` 195 @*/ 196 PetscErrorCode MatCreateSNESMFMore(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 MatSNESMFMoreSetParameters - Sets the parameters for the approximation of 264 matrix-vector products using finite differences, see `MatCreateSNESMFMore()` 265 266 Input Parameters: 267 + mat - the matrix 268 . error_rel - relative error (should be set to the square root of the relative error in the function evaluations) 269 . umin - minimum allowable u-value 270 - h - differencing parameter 271 272 Options Database Keys: 273 + -snes_mf_err <error_rel> - see `MatCreateSNESMF()` 274 . -snes_mf_umin <umin> - see `MatCreateSNESMF()` 275 . -snes_mf_compute_err - compute the square root or relative error in function 276 . -snes_mf_freq_err <freq> - set the frequency to recompute the parameters 277 - -snes_mf_jorge - use the method of Jorge More 278 279 Level: advanced 280 281 Note: 282 If the user sets the parameter h directly, then this value will be used 283 instead of the default computation as discussed in `MatCreateSNESMFMore()` 284 285 .seealso: `MatCreateSNESMF()`, `MatCreateSNESMFMore()` 286 @*/ 287 PetscErrorCode MatSNESMFMoreSetParameters(Mat mat, PetscReal error, PetscReal umin, PetscReal h) { 288 MFCtx_Private *ctx; 289 290 PetscFunctionBegin; 291 PetscCall(MatShellGetContext(mat, &ctx)); 292 if (ctx) { 293 if (error != PETSC_DEFAULT) ctx->error_rel = error; 294 if (umin != PETSC_DEFAULT) ctx->umin = umin; 295 if (h != PETSC_DEFAULT) { 296 ctx->h = h; 297 ctx->need_h = PETSC_FALSE; 298 } 299 } 300 PetscFunctionReturn(0); 301 } 302 303 PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes) { 304 MFCtx_Private *ctx; 305 Mat mat; 306 307 PetscFunctionBegin; 308 PetscCall(SNESGetJacobian(snes, &mat, NULL, NULL, NULL)); 309 PetscCall(MatShellGetContext(mat, &ctx)); 310 if (ctx) ctx->need_h = PETSC_TRUE; 311 PetscFunctionReturn(0); 312 } 313