Actual source code: snesmfj2.c
1: #include <petsc/private/snesimpl.h>
2: /* matimpl.h is needed only for logging of matrix operation */
3: #include <petsc/private/matimpl.h>
5: PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES);
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 *);
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;
27: static PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat)
28: {
29: MFCtx_Private *ctx;
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: }
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 iascii;
48: PetscFunctionBegin;
49: PetscCall(MatShellGetContext(J, &ctx));
50: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
51: if (iascii) {
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: }
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);
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));
87: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
88: PetscCall(MatShellGetContext(mat, &ctx));
89: snes = ctx->snes;
90: w = ctx->w;
91: umin = ctx->umin;
93: PetscCall(SNESGetSolution(snes, &U));
94: eval_fct = SNESComputeFunction;
95: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
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));
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));
111: ctx->error_rel = PetscSqrtReal(noise);
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));
115: ctx->compute_err_iter = iter;
116: ctx->need_err = PETSC_FALSE;
117: }
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));
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;
136: if (!ctx->jorge || !ctx->need_h) PetscCall(PetscInfo(snes, "h = %g\n", (double)h));
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));
146: PetscCall(PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
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()`.
155: Input Parameters:
156: + snes - the `SNES` context
157: - x - vector where `SNES` solution is to be stored.
159: Output Parameter:
160: . J - the matrix-free matrix
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
169: Level: advanced
171: Notes:
172: This is an experimental approach, use `MatCreateSNESMF()`.
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
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
191: The user can set these parameters via `MatMFFDSetFunctionError()`.
193: The user should call `MatDestroy()` when finished with the matrix-free
194: matrix context.
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];
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;
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) {
232: PetscCall(SNESDiffParameterCreate_More(snes, x, &mfctx->data));
233: } else mfctx->data = NULL;
235: PetscCall(PetscOptionsHasHelp(((PetscObject)snes)->options, &flg));
236: PetscCall(PetscStrncpy(p, "-", sizeof(p)));
237: if (((PetscObject)snes)->prefix) PetscCall(PetscStrlcat(p, ((PetscObject)snes)->prefix, sizeof(p)));
238: if (flg) {
239: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " Matrix-free Options (via SNES):\n"));
240: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_err <err>: set sqrt of relative error in function (default %g)\n", p, (double)mfctx->error_rel));
241: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_umin <umin>: see users manual (default %g)\n", p, (double)mfctx->umin));
242: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_jorge: use Jorge More's method\n", p));
243: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_compute_err: compute sqrt or relative error in function\n", p));
244: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n", p));
245: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_noise_file <file>: set file for printing noise info\n", p));
246: }
247: PetscCall(VecDuplicate(x, &mfctx->w));
248: PetscCall(PetscObjectGetComm((PetscObject)x, &comm));
249: PetscCall(VecGetSize(x, &n));
250: PetscCall(VecGetLocalSize(x, &nloc));
251: PetscCall(MatCreate(comm, J));
252: PetscCall(MatSetSizes(*J, nloc, n, n, n));
253: PetscCall(MatSetType(*J, MATSHELL));
254: PetscCall(MatShellSetContext(*J, mfctx));
255: PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))SNESMatrixFreeMult2_Private));
256: PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))SNESMatrixFreeDestroy2_Private));
257: PetscCall(MatShellSetOperation(*J, MATOP_VIEW, (void (*)(void))SNESMatrixFreeView2_Private));
258: PetscCall(MatSetUp(*J));
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: /*@
263: MatSNESMFMoreSetParameters - Sets the parameters for the approximation of
264: matrix-vector products using finite differences, see `MatCreateSNESMFMore()`
266: Input Parameters:
267: + mat - the matrix
268: . error - 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
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
279: Level: advanced
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()`
285: .seealso: [](ch_snes), `SNES`, `MatCreateSNESMF()`, `MatCreateSNESMFMore()`
286: @*/
287: PetscErrorCode MatSNESMFMoreSetParameters(Mat mat, PetscReal error, PetscReal umin, PetscReal h)
288: {
289: MFCtx_Private *ctx;
291: PetscFunctionBegin;
292: PetscCall(MatShellGetContext(mat, &ctx));
293: if (ctx) {
294: if (error != (PetscReal)PETSC_DEFAULT) ctx->error_rel = error;
295: if (umin != (PetscReal)PETSC_DEFAULT) ctx->umin = umin;
296: if (h != (PetscReal)PETSC_DEFAULT) {
297: ctx->h = h;
298: ctx->need_h = PETSC_FALSE;
299: }
300: }
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes)
305: {
306: MFCtx_Private *ctx;
307: Mat mat;
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(PETSC_SUCCESS);
314: }