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);
7: PETSC_INTERN PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*);
8: PETSC_INTERN PetscErrorCode SNESDefaultMatrixFreeSetParameters2(Mat,PetscReal,PetscReal,PetscReal);
10: PETSC_INTERN PetscErrorCode SNESDiffParameterCreate_More(SNES,Vec,void**);
11: PETSC_INTERN PetscErrorCode SNESDiffParameterCompute_More(SNES,void*,Vec,Vec,PetscReal*,PetscReal*);
12: PETSC_INTERN PetscErrorCode SNESDiffParameterDestroy_More(void*);
14: typedef struct { /* default context for matrix-free SNES */
15: SNES snes; /* SNES context */
16: Vec w; /* work vector */
17: MatNullSpace sp; /* null space context */
18: PetscReal error_rel; /* square root of relative error in computing function */
19: PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */
20: PetscBool jorge; /* flag indicating use of Jorge's method for determining the differencing parameter */
21: PetscReal h; /* differencing parameter */
22: PetscBool need_h; /* flag indicating whether we must compute h */
23: PetscBool need_err; /* flag indicating whether we must currently compute error_rel */
24: PetscBool compute_err; /* flag indicating whether we must ever compute error_rel */
25: PetscInt compute_err_iter; /* last iter where we've computer error_rel */
26: PetscInt compute_err_freq; /* frequency of computing error_rel */
27: void *data; /* implementation-specific data */
28: } MFCtx_Private;
30: PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat)
31: {
32: MFCtx_Private *ctx;
34: MatShellGetContext(mat,&ctx);
35: VecDestroy(&ctx->w);
36: MatNullSpaceDestroy(&ctx->sp);
37: if (ctx->jorge || ctx->compute_err) SNESDiffParameterDestroy_More(ctx->data);
38: PetscFree(ctx);
39: return 0;
40: }
42: /*
43: SNESMatrixFreeView2_Private - Views matrix-free parameters.
44: */
45: PetscErrorCode SNESMatrixFreeView2_Private(Mat J,PetscViewer viewer)
46: {
47: MFCtx_Private *ctx;
48: PetscBool iascii;
50: MatShellGetContext(J,&ctx);
51: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
52: if (iascii) {
53: PetscViewerASCIIPrintf(viewer," SNES matrix-free approximation:\n");
54: if (ctx->jorge) {
55: PetscViewerASCIIPrintf(viewer," using Jorge's method of determining differencing parameter\n");
56: }
57: PetscViewerASCIIPrintf(viewer," err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
58: PetscViewerASCIIPrintf(viewer," umin=%g (minimum iterate parameter)\n",(double)ctx->umin);
59: if (ctx->compute_err) {
60: PetscViewerASCIIPrintf(viewer," freq_err=%D (frequency for computing err)\n",ctx->compute_err_freq);
61: }
62: }
63: return 0;
64: }
66: /*
67: SNESMatrixFreeMult2_Private - Default matrix-free form for Jacobian-vector
68: product, y = F'(u)*a:
69: y = (F(u + ha) - F(u)) /h,
70: where F = nonlinear function, as set by SNESSetFunction()
71: u = current iterate
72: h = difference interval
73: */
74: PetscErrorCode SNESMatrixFreeMult2_Private(Mat mat,Vec a,Vec y)
75: {
76: MFCtx_Private *ctx;
77: SNES snes;
78: PetscReal h,norm,sum,umin,noise;
79: PetscScalar hs,dot;
80: Vec w,U,F;
81: MPI_Comm comm;
82: PetscInt iter;
83: PetscErrorCode (*eval_fct)(SNES,Vec,Vec);
85: /* We log matrix-free matrix-vector products separately, so that we can
86: separate the performance monitoring from the cases that use conventional
87: storage. We may eventually modify event logging to associate events
88: with particular objects, hence alleviating the more general problem. */
89: PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);
91: PetscObjectGetComm((PetscObject)mat,&comm);
92: MatShellGetContext(mat,&ctx);
93: snes = ctx->snes;
94: w = ctx->w;
95: umin = ctx->umin;
97: SNESGetSolution(snes,&U);
98: eval_fct = SNESComputeFunction;
99: SNESGetFunction(snes,&F,NULL,NULL);
101: /* Determine a "good" step size, h */
102: if (ctx->need_h) {
104: /* Use Jorge's method to compute h */
105: if (ctx->jorge) {
106: SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);
108: /* Use the Brown/Saad method to compute h */
109: } else {
110: /* Compute error if desired */
111: SNESGetIterationNumber(snes,&iter);
112: if ((ctx->need_err) || ((ctx->compute_err_freq) && (ctx->compute_err_iter != iter) && (!((iter-1)%ctx->compute_err_freq)))) {
113: /* Use Jorge's method to compute noise */
114: SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);
116: ctx->error_rel = PetscSqrtReal(noise);
118: PetscInfo(snes,"Using Jorge's noise: noise=%g, sqrt(noise)=%g, h_more=%g\n",(double)noise,(double)ctx->error_rel,(double)h);
120: ctx->compute_err_iter = iter;
121: ctx->need_err = PETSC_FALSE;
122: }
124: VecDotBegin(U,a,&dot);
125: VecNormBegin(a,NORM_1,&sum);
126: VecNormBegin(a,NORM_2,&norm);
127: VecDotEnd(U,a,&dot);
128: VecNormEnd(a,NORM_1,&sum);
129: VecNormEnd(a,NORM_2,&norm);
131: /* Safeguard for step sizes too small */
132: if (sum == 0.0) {
133: dot = 1.0;
134: norm = 1.0;
135: } else if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum;
136: else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum;
137: h = PetscRealPart(ctx->error_rel*dot/(norm*norm));
138: }
139: } else h = ctx->h;
141: if (!ctx->jorge || !ctx->need_h) PetscInfo(snes,"h = %g\n",(double)h);
143: /* Evaluate function at F(u + ha) */
144: hs = h;
145: VecWAXPY(w,hs,a,U);
146: eval_fct(snes,w,y);
147: VecAXPY(y,-1.0,F);
148: VecScale(y,1.0/hs);
149: if (mat->nullsp) MatNullSpaceRemove(mat->nullsp,y);
151: PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
152: return 0;
153: }
155: /*@C
156: SNESMatrixFreeCreate2 - Creates a matrix-free matrix
157: context for use with a SNES solver. This matrix can be used as
158: the Jacobian argument for the routine SNESSetJacobian().
160: Input Parameters:
161: + snes - the SNES context
162: - x - vector where SNES solution is to be stored.
164: Output Parameter:
165: . J - the matrix-free matrix
167: Level: advanced
169: Notes:
170: The matrix-free matrix context merely contains the function pointers
171: and work space for performing finite difference approximations of
172: Jacobian-vector products, J(u)*a, via
174: $ J(u)*a = [J(u+h*a) - J(u)]/h,
175: $ where by default
176: $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
177: $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
178: $ where
179: $ error_rel = square root of relative error in
180: $ function evaluation
181: $ umin = minimum iterate parameter
182: $ Alternatively, the differencing parameter, h, can be set using
183: $ Jorge's nifty new strategy if one specifies the option
184: $ -snes_mf_jorge
186: The user can set these parameters via MatMFFDSetFunctionError().
187: See Users-Manual: ch_snes for details
189: The user should call MatDestroy() when finished with the matrix-free
190: matrix context.
192: Options Database Keys:
193: $ -snes_mf_err <error_rel>
194: $ -snes_mf_unim <umin>
195: $ -snes_mf_compute_err
196: $ -snes_mf_freq_err <freq>
197: $ -snes_mf_jorge
199: .seealso: MatDestroy(), MatMFFDSetFunctionError()
200: @*/
201: PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES snes,Vec x,Mat *J)
202: {
203: MPI_Comm comm;
204: MFCtx_Private *mfctx;
205: PetscInt n,nloc;
206: PetscBool flg;
207: char p[64];
209: PetscNewLog(snes,&mfctx);
210: mfctx->sp = NULL;
211: mfctx->snes = snes;
212: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
213: mfctx->umin = 1.e-6;
214: mfctx->h = 0.0;
215: mfctx->need_h = PETSC_TRUE;
216: mfctx->need_err = PETSC_FALSE;
217: mfctx->compute_err = PETSC_FALSE;
218: mfctx->compute_err_freq = 0;
219: mfctx->compute_err_iter = -1;
220: mfctx->compute_err = PETSC_FALSE;
221: mfctx->jorge = PETSC_FALSE;
223: PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_err",&mfctx->error_rel,NULL);
224: PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_umin",&mfctx->umin,NULL);
225: PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_jorge",&mfctx->jorge,NULL);
226: PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_compute_err",&mfctx->compute_err,NULL);
227: PetscOptionsGetInt(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_freq_err",&mfctx->compute_err_freq,&flg);
228: if (flg) {
229: if (mfctx->compute_err_freq < 0) mfctx->compute_err_freq = 0;
230: mfctx->compute_err = PETSC_TRUE;
231: }
232: if (mfctx->compute_err) mfctx->need_err = PETSC_TRUE;
233: if (mfctx->jorge || mfctx->compute_err) {
234: SNESDiffParameterCreate_More(snes,x,&mfctx->data);
235: } else mfctx->data = NULL;
237: PetscOptionsHasHelp(((PetscObject)snes)->options,&flg);
238: PetscStrncpy(p,"-",sizeof(p));
239: if (((PetscObject)snes)->prefix) PetscStrlcat(p,((PetscObject)snes)->prefix,sizeof(p));
240: if (flg) {
241: PetscPrintf(PetscObjectComm((PetscObject)snes)," Matrix-free Options (via SNES):\n");
242: PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_err <err>: set sqrt of relative error in function (default %g)\n",p,(double)mfctx->error_rel);
243: PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_umin <umin>: see users manual (default %g)\n",p,(double)mfctx->umin);
244: PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_jorge: use Jorge More's method\n",p);
245: PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p);
246: PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p);
247: PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_noise_file <file>: set file for printing noise info\n",p);
248: }
249: VecDuplicate(x,&mfctx->w);
250: PetscObjectGetComm((PetscObject)x,&comm);
251: VecGetSize(x,&n);
252: VecGetLocalSize(x,&nloc);
253: MatCreate(comm,J);
254: MatSetSizes(*J,nloc,n,n,n);
255: MatSetType(*J,MATSHELL);
256: MatShellSetContext(*J,mfctx);
257: MatShellSetOperation(*J,MATOP_MULT,(void (*)(void))SNESMatrixFreeMult2_Private);
258: MatShellSetOperation(*J,MATOP_DESTROY,(void (*)(void))SNESMatrixFreeDestroy2_Private);
259: MatShellSetOperation(*J,MATOP_VIEW,(void (*)(void))SNESMatrixFreeView2_Private);
260: MatSetUp(*J);
262: PetscLogObjectParent((PetscObject)*J,(PetscObject)mfctx->w);
263: PetscLogObjectParent((PetscObject)snes,(PetscObject)*J);
264: return 0;
265: }
267: /*@C
268: SNESDefaultMatrixFreeSetParameters2 - Sets the parameters for the approximation of
269: matrix-vector products using finite differences.
271: $ J(u)*a = [J(u+h*a) - J(u)]/h where
273: either the user sets h directly here, or this parameter is computed via
275: $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
276: $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
277: $
279: Input Parameters:
280: + mat - the matrix
281: . error_rel - relative error (should be set to the square root of
282: the relative error in the function evaluations)
283: . umin - minimum allowable u-value
284: - h - differencing parameter
286: Level: advanced
288: Notes:
289: If the user sets the parameter h directly, then this value will be used
290: instead of the default computation indicated above.
292: .seealso: MatCreateSNESMF()
293: @*/
294: PetscErrorCode SNESDefaultMatrixFreeSetParameters2(Mat mat,PetscReal error,PetscReal umin,PetscReal h)
295: {
296: MFCtx_Private *ctx;
298: MatShellGetContext(mat,&ctx);
299: if (ctx) {
300: if (error != PETSC_DEFAULT) ctx->error_rel = error;
301: if (umin != PETSC_DEFAULT) ctx->umin = umin;
302: if (h != PETSC_DEFAULT) {
303: ctx->h = h;
304: ctx->need_h = PETSC_FALSE;
305: }
306: }
307: return 0;
308: }
310: PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes)
311: {
312: MFCtx_Private *ctx;
313: Mat mat;
315: SNESGetJacobian(snes,&mat,NULL,NULL,NULL);
316: MatShellGetContext(mat,&ctx);
317: if (ctx) ctx->need_h = PETSC_TRUE;
318: return 0;
319: }