Actual source code: snesmfj2.c


  2: #include <petsc/private/snesimpl.h>
  3: /* matimpl.h is needed only for logging of matrix operation */
  4: #include <petsc/private/matimpl.h>

  6: PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES);

  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 *);

 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;

 28: PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat)
 29: {
 30:   MFCtx_Private *ctx;

 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: }

 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;

 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: }

 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);

 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));

 88:   PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
 89:   PetscCall(MatShellGetContext(mat, &ctx));
 90:   snes = ctx->snes;
 91:   w    = ctx->w;
 92:   umin = ctx->umin;

 94:   PetscCall(SNESGetSolution(snes, &U));
 95:   eval_fct = SNESComputeFunction;
 96:   PetscCall(SNESGetFunction(snes, &F, NULL, NULL));

 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));

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));

112:         ctx->error_rel = PetscSqrtReal(noise);

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));

116:         ctx->compute_err_iter = iter;
117:         ctx->need_err         = PETSC_FALSE;
118:       }

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));

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;

137:   if (!ctx->jorge || !ctx->need_h) PetscCall(PetscInfo(snes, "h = %g\n", (double)h));

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));

147:   PetscCall(PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0));
148:   PetscFunctionReturn(PETSC_SUCCESS);
149: }

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()`.

156:    Input Parameters:
157: +  snes - the `SNES` context
158: -  x - vector where `SNES` solution is to be stored.

160:    Output Parameter:
161: .  J - the matrix-free matrix

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

170:    Level: advanced

172:    Notes:
173:    This is an experimental approach, use `MatCreateSNESMF()`.

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

179: .vb
180:        J(u)*a = [J(u+h*a) - J(u)]/h,
181:    where by default
182:         h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
183:           = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
184:    where
185:         error_rel = square root of relative error in 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: .ve

192:    The user can set these parameters via `MatMFFDSetFunctionError()`.

194:    The user should call `MatDestroy()` when finished with the matrix-free
195:    matrix context.

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];

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;

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;

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));

261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }

264: /*@C
265:    MatSNESMFMoreSetParameters - Sets the parameters for the approximation of
266:    matrix-vector products using finite differences, see  `MatCreateSNESMFMore()`

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

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

281:    Level: advanced

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()`

287: .seealso: `MatCreateSNESMF()`, `MatCreateSNESMFMore()`
288: @*/
289: PetscErrorCode MatSNESMFMoreSetParameters(Mat mat, PetscReal error, PetscReal umin, PetscReal h)
290: {
291:   MFCtx_Private *ctx;

293:   PetscFunctionBegin;
294:   PetscCall(MatShellGetContext(mat, &ctx));
295:   if (ctx) {
296:     if (error != (PetscReal)PETSC_DEFAULT) ctx->error_rel = error;
297:     if (umin != (PetscReal)PETSC_DEFAULT) ctx->umin = umin;
298:     if (h != (PetscReal)PETSC_DEFAULT) {
299:       ctx->h      = h;
300:       ctx->need_h = PETSC_FALSE;
301:     }
302:   }
303:   PetscFunctionReturn(PETSC_SUCCESS);
304: }

306: PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes)
307: {
308:   MFCtx_Private *ctx;
309:   Mat            mat;

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: }