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