1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
2 /* matimpl.h is needed only for logging of matrix operation */
3 #include <petsc/private/matimpl.h>
4
5 PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES);
6
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 *);
10
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;
26
SNESMatrixFreeDestroy2_Private(Mat mat)27 static PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat)
28 {
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(PETSC_SUCCESS);
38 }
39
40 /*
41 SNESMatrixFreeView2_Private - Views matrix-free parameters.
42 */
SNESMatrixFreeView2_Private(Mat J,PetscViewer viewer)43 static PetscErrorCode SNESMatrixFreeView2_Private(Mat J, PetscViewer viewer)
44 {
45 MFCtx_Private *ctx;
46 PetscBool isascii;
47
48 PetscFunctionBegin;
49 PetscCall(MatShellGetContext(J, &ctx));
50 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
51 if (isascii) {
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 }
60
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 */
SNESMatrixFreeMult2_Private(Mat mat,Vec a,Vec y)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);
79
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));
86
87 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
88 PetscCall(MatShellGetContext(mat, &ctx));
89 snes = ctx->snes;
90 w = ctx->w;
91 umin = ctx->umin;
92
93 PetscCall(SNESGetSolution(snes, &U));
94 eval_fct = SNESComputeFunction;
95 PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
96
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));
102
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));
110
111 ctx->error_rel = PetscSqrtReal(noise);
112
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));
114
115 ctx->compute_err_iter = iter;
116 ctx->need_err = PETSC_FALSE;
117 }
118
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));
125
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;
135
136 if (!ctx->jorge || !ctx->need_h) PetscCall(PetscInfo(snes, "h = %g\n", (double)h));
137
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));
145
146 PetscCall(PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0));
147 PetscFunctionReturn(PETSC_SUCCESS);
148 }
149
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()`.
154
155 Input Parameters:
156 + snes - the `SNES` context
157 - x - vector where `SNES` solution is to be stored.
158
159 Output Parameter:
160 . J - the matrix-free matrix
161
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
168
169 Level: advanced
170
171 Notes:
172 This is an experimental approach, use `MatCreateSNESMF()`.
173
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
177
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
190
191 The user can set these parameters via `MatMFFDSetFunctionError()`.
192
193 The user should call `MatDestroy()` when finished with the matrix-free
194 matrix context.
195
196 .seealso: [](ch_snes), `SNESCreateMF()`, `MatCreateMFFD()`, `MatDestroy()`, `MatMFFDSetFunctionError()`
197 @*/
MatCreateSNESMFMore(SNES snes,Vec x,Mat * J)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];
205
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;
220
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) PetscCall(SNESDiffParameterCreate_More(snes, x, &mfctx->data));
232 else mfctx->data = NULL;
233
234 PetscCall(PetscOptionsHasHelp(((PetscObject)snes)->options, &flg));
235 PetscCall(PetscStrncpy(p, "-", sizeof(p)));
236 if (((PetscObject)snes)->prefix) PetscCall(PetscStrlcat(p, ((PetscObject)snes)->prefix, sizeof(p)));
237 if (flg) {
238 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " Matrix-free Options (via SNES):\n"));
239 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_err <err>: set sqrt of relative error in function (default %g)\n", p, (double)mfctx->error_rel));
240 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_umin <umin>: see users manual (default %g)\n", p, (double)mfctx->umin));
241 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_jorge: use Jorge More's method\n", p));
242 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_compute_err: compute sqrt or relative error in function\n", p));
243 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n", p));
244 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " %ssnes_mf_noise_file <file>: set file for printing noise info\n", p));
245 }
246 PetscCall(VecDuplicate(x, &mfctx->w));
247 PetscCall(PetscObjectGetComm((PetscObject)x, &comm));
248 PetscCall(VecGetSize(x, &n));
249 PetscCall(VecGetLocalSize(x, &nloc));
250 PetscCall(MatCreate(comm, J));
251 PetscCall(MatSetSizes(*J, nloc, n, n, n));
252 PetscCall(MatSetType(*J, MATSHELL));
253 PetscCall(MatShellSetContext(*J, mfctx));
254 PetscCall(MatShellSetOperation(*J, MATOP_MULT, (PetscErrorCodeFn *)SNESMatrixFreeMult2_Private));
255 PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (PetscErrorCodeFn *)SNESMatrixFreeDestroy2_Private));
256 PetscCall(MatShellSetOperation(*J, MATOP_VIEW, (PetscErrorCodeFn *)SNESMatrixFreeView2_Private));
257 PetscCall(MatSetUp(*J));
258 PetscFunctionReturn(PETSC_SUCCESS);
259 }
260
261 /*@
262 MatSNESMFMoreSetParameters - Sets the parameters for the approximation of
263 matrix-vector products using finite differences, see `MatCreateSNESMFMore()`
264
265 Input Parameters:
266 + mat - the matrix
267 . error - relative error (should be set to the square root of the relative error in the function evaluations)
268 . umin - minimum allowable u-value
269 - h - differencing parameter
270
271 Options Database Keys:
272 + -snes_mf_err <error_rel> - see `MatCreateSNESMF()`
273 . -snes_mf_umin <umin> - see `MatCreateSNESMF()`
274 . -snes_mf_compute_err - compute the square root or relative error in function
275 . -snes_mf_freq_err <freq> - set the frequency to recompute the parameters
276 - -snes_mf_jorge - use the method of Jorge More
277
278 Level: advanced
279
280 Note:
281 If the user sets the parameter `h` directly, then this value will be used
282 instead of the default computation as discussed in `MatCreateSNESMFMore()`
283
284 .seealso: [](ch_snes), `SNES`, `MatCreateSNESMF()`, `MatCreateSNESMFMore()`
285 @*/
MatSNESMFMoreSetParameters(Mat mat,PetscReal error,PetscReal umin,PetscReal h)286 PetscErrorCode MatSNESMFMoreSetParameters(Mat mat, PetscReal error, PetscReal umin, PetscReal h)
287 {
288 MFCtx_Private *ctx;
289
290 PetscFunctionBegin;
291 PetscCall(MatShellGetContext(mat, &ctx));
292 if (ctx) {
293 if (error != (PetscReal)PETSC_DEFAULT) ctx->error_rel = error;
294 if (umin != (PetscReal)PETSC_DEFAULT) ctx->umin = umin;
295 if (h != (PetscReal)PETSC_DEFAULT) {
296 ctx->h = h;
297 ctx->need_h = PETSC_FALSE;
298 }
299 }
300 PetscFunctionReturn(PETSC_SUCCESS);
301 }
302
SNESUnSetMatrixFreeParameter(SNES snes)303 PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes)
304 {
305 MFCtx_Private *ctx;
306 Mat mat;
307
308 PetscFunctionBegin;
309 PetscCall(SNESGetJacobian(snes, &mat, NULL, NULL, NULL));
310 PetscCall(MatShellGetContext(mat, &ctx));
311 if (ctx) ctx->need_h = PETSC_TRUE;
312 PetscFunctionReturn(PETSC_SUCCESS);
313 }
314