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