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: MatShellGetContext(mat, &ctx);
33: VecDestroy(&ctx->w);
34: MatNullSpaceDestroy(&ctx->sp);
35: if (ctx->jorge || ctx->compute_err) SNESDiffParameterDestroy_More(ctx->data);
36: PetscFree(ctx);
37: return 0;
38: }
40: /*
41: SNESMatrixFreeView2_Private - Views matrix-free parameters.
42: */
43: PetscErrorCode SNESMatrixFreeView2_Private(Mat J, PetscViewer viewer)
44: {
45: MFCtx_Private *ctx;
46: PetscBool iascii;
48: MatShellGetContext(J, &ctx);
49: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
50: if (iascii) {
51: PetscViewerASCIIPrintf(viewer, " SNES matrix-free approximation:\n");
52: if (ctx->jorge) PetscViewerASCIIPrintf(viewer, " using Jorge's method of determining differencing parameter\n");
53: PetscViewerASCIIPrintf(viewer, " err=%g (relative error in function evaluation)\n", (double)ctx->error_rel);
54: PetscViewerASCIIPrintf(viewer, " umin=%g (minimum iterate parameter)\n", (double)ctx->umin);
55: if (ctx->compute_err) PetscViewerASCIIPrintf(viewer, " freq_err=%" PetscInt_FMT " (frequency for computing err)\n", ctx->compute_err_freq);
56: }
57: return 0;
58: }
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: {
70: MFCtx_Private *ctx;
71: SNES snes;
72: PetscReal h, norm, sum, umin, noise;
73: PetscScalar hs, dot;
74: Vec w, U, F;
75: MPI_Comm comm;
76: PetscInt iter;
77: PetscErrorCode (*eval_fct)(SNES, Vec, Vec);
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: PetscLogEventBegin(MATMFFD_Mult, a, y, 0, 0);
85: PetscObjectGetComm((PetscObject)mat, &comm);
86: MatShellGetContext(mat, &ctx);
87: snes = ctx->snes;
88: w = ctx->w;
89: umin = ctx->umin;
91: SNESGetSolution(snes, &U);
92: eval_fct = SNESComputeFunction;
93: SNESGetFunction(snes, &F, NULL, NULL);
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: SNESDiffParameterCompute_More(snes, ctx->data, U, a, &noise, &h);
101: /* Use the Brown/Saad method to compute h */
102: } else {
103: /* Compute error if desired */
104: 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: SNESDiffParameterCompute_More(snes, ctx->data, U, a, &noise, &h);
109: ctx->error_rel = PetscSqrtReal(noise);
111: PetscInfo(snes, "Using Jorge's noise: noise=%g, sqrt(noise)=%g, h_more=%g\n", (double)noise, (double)ctx->error_rel, (double)h);
113: ctx->compute_err_iter = iter;
114: ctx->need_err = PETSC_FALSE;
115: }
117: VecDotBegin(U, a, &dot);
118: VecNormBegin(a, NORM_1, &sum);
119: VecNormBegin(a, NORM_2, &norm);
120: VecDotEnd(U, a, &dot);
121: VecNormEnd(a, NORM_1, &sum);
122: VecNormEnd(a, NORM_2, &norm);
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;
134: if (!ctx->jorge || !ctx->need_h) PetscInfo(snes, "h = %g\n", (double)h);
136: /* Evaluate function at F(u + ha) */
137: hs = h;
138: VecWAXPY(w, hs, a, U);
139: eval_fct(snes, w, y);
140: VecAXPY(y, -1.0, F);
141: VecScale(y, 1.0 / hs);
142: if (mat->nullsp) MatNullSpaceRemove(mat->nullsp, y);
144: PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0);
145: return 0;
146: }
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()`.
153: Input Parameters:
154: + snes - the `SNES` context
155: - x - vector where `SNES` solution is to be stored.
157: Output Parameter:
158: . J - the matrix-free matrix
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
167: Level: advanced
169: Notes:
170: This is an experimental approach, use `MatCreateSNESMF()`.
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
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
188: The user can set these parameters via `MatMFFDSetFunctionError()`.
189: See Users-Manual: ch_snes for details
191: The user should call `MatDestroy()` when finished with the matrix-free
192: matrix context.
194: .seealso: `SNESCreateMF`(), `MatCreateMFFD()`, `MatDestroy()`, `MatMFFDSetFunctionError()`
195: @*/
196: PetscErrorCode MatCreateSNESMFMore(SNES snes, Vec x, Mat *J)
197: {
198: MPI_Comm comm;
199: MFCtx_Private *mfctx;
200: PetscInt n, nloc;
201: PetscBool flg;
202: char p[64];
204: PetscNew(&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;
218: PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_err", &mfctx->error_rel, NULL);
219: PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_umin", &mfctx->umin, NULL);
220: PetscOptionsGetBool(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_jorge", &mfctx->jorge, NULL);
221: PetscOptionsGetBool(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_compute_err", &mfctx->compute_err, NULL);
222: 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: SNESDiffParameterCreate_More(snes, x, &mfctx->data);
230: } else mfctx->data = NULL;
232: PetscOptionsHasHelp(((PetscObject)snes)->options, &flg);
233: PetscStrncpy(p, "-", sizeof(p));
234: if (((PetscObject)snes)->prefix) PetscStrlcat(p, ((PetscObject)snes)->prefix, sizeof(p));
235: if (flg) {
236: PetscPrintf(PetscObjectComm((PetscObject)snes), " Matrix-free Options (via SNES):\n");
237: PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_err <err>: set sqrt of relative error in function (default %g)\n", p, (double)mfctx->error_rel);
238: PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_umin <umin>: see users manual (default %g)\n", p, (double)mfctx->umin);
239: PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_jorge: use Jorge More's method\n", p);
240: PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_compute_err: compute sqrt or relative error in function\n", p);
241: PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n", p);
242: PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_noise_file <file>: set file for printing noise info\n", p);
243: }
244: VecDuplicate(x, &mfctx->w);
245: PetscObjectGetComm((PetscObject)x, &comm);
246: VecGetSize(x, &n);
247: VecGetLocalSize(x, &nloc);
248: MatCreate(comm, J);
249: MatSetSizes(*J, nloc, n, n, n);
250: MatSetType(*J, MATSHELL);
251: MatShellSetContext(*J, mfctx);
252: MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))SNESMatrixFreeMult2_Private);
253: MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))SNESMatrixFreeDestroy2_Private);
254: MatShellSetOperation(*J, MATOP_VIEW, (void (*)(void))SNESMatrixFreeView2_Private);
255: MatSetUp(*J);
257: return 0;
258: }
260: /*@C
261: MatSNESMFMoreSetParameters - Sets the parameters for the approximation of
262: matrix-vector products using finite differences, see `MatCreateSNESMFMore()`
264: Input Parameters:
265: + mat - the matrix
266: . error_rel - relative error (should be set to the square root of the relative error in the function evaluations)
267: . umin - minimum allowable u-value
268: - h - differencing parameter
270: Options Database Keys:
271: + -snes_mf_err <error_rel> - see `MatCreateSNESMF()`
272: . -snes_mf_umin <umin> - see `MatCreateSNESMF()`
273: . -snes_mf_compute_err - compute the square root or relative error in function
274: . -snes_mf_freq_err <freq> - set the frequency to recompute the parameters
275: - -snes_mf_jorge - use the method of Jorge More
277: Level: advanced
279: Note:
280: If the user sets the parameter h directly, then this value will be used
281: instead of the default computation as discussed in `MatCreateSNESMFMore()`
283: .seealso: `MatCreateSNESMF()`, `MatCreateSNESMFMore()`
284: @*/
285: PetscErrorCode MatSNESMFMoreSetParameters(Mat mat, PetscReal error, PetscReal umin, PetscReal h)
286: {
287: MFCtx_Private *ctx;
289: MatShellGetContext(mat, &ctx);
290: if (ctx) {
291: if (error != PETSC_DEFAULT) ctx->error_rel = error;
292: if (umin != PETSC_DEFAULT) ctx->umin = umin;
293: if (h != PETSC_DEFAULT) {
294: ctx->h = h;
295: ctx->need_h = PETSC_FALSE;
296: }
297: }
298: return 0;
299: }
301: PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes)
302: {
303: MFCtx_Private *ctx;
304: Mat mat;
306: SNESGetJacobian(snes, &mat, NULL, NULL, NULL);
307: MatShellGetContext(mat, &ctx);
308: if (ctx) ctx->need_h = PETSC_TRUE;
309: return 0;
310: }