Actual source code: snesmfjdef.c
1: #define PETSCSNES_DLL
3: /*
4: Implements the DS PETSc approach for computing the h
5: parameter used with the finite difference based matrix-free
6: Jacobian-vector products.
8: To make your own: clone this file and modify for your needs.
10: Mandatory functions:
11: -------------------
12: MatSNESMFCompute_ - for a given point and direction computes h
14: MatSNESMFCreate_ - fills in the MatSNESMF data structure
15: for this particular implementation
17:
18: Optional functions:
19: -------------------
20: MatSNESMFView_ - prints information about the parameters being used.
21: This is called when SNESView() or -snes_view is used.
23: MatSNESMFSetFromOptions_ - checks the options database for options that
24: apply to this method.
26: MatSNESMFDestroy_ - frees any space allocated by the routines above
28: */
30: /*
31: This include file defines the data structure MatSNESMF that
32: includes information about the computation of h. It is shared by
33: all implementations that people provide
34: */
35: #include include/private/matimpl.h
36: #include src/snes/mf/snesmfj.h
38: /*
39: The method has one parameter that is used to
40: "cutoff" very small values. This is stored in a data structure
41: that is only visible to this file. If your method has no parameters
42: it can omit this, if it has several simply reorganize the data structure.
43: The data structure is "hung-off" the MatSNESMF data structure in
44: the void *hctx; field.
45: */
46: typedef struct {
47: PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */
48: } MatSNESMF_DS;
52: /*
53: MatSNESMFCompute_DS - Standard PETSc code for computing the
54: differencing paramter (h) for use with matrix-free finite differences.
56: Input Parameters:
57: + ctx - the matrix free context
58: . U - the location at which you want the Jacobian
59: - a - the direction you want the derivative
61:
62: Output Parameter:
63: . h - the scale computed
65: */
66: static PetscErrorCode MatSNESMFCompute_DS(MatSNESMFCtx ctx,Vec U,Vec a,PetscScalar *h,PetscTruth *zeroa)
67: {
68: MatSNESMF_DS *hctx = (MatSNESMF_DS*)ctx->hctx;
69: PetscReal nrm,sum,umin = hctx->umin;
70: PetscScalar dot;
71: PetscErrorCode ierr;
74: if (!(ctx->count % ctx->recomputeperiod)) {
75: /*
76: This algorithm requires 2 norms and 1 inner product. Rather than
77: use directly the VecNorm() and VecDot() routines (and thus have
78: three separate collective operations, we use the VecxxxBegin/End() routines
79: */
80: VecDotBegin(U,a,&dot);
81: VecNormBegin(a,NORM_1,&sum);
82: VecNormBegin(a,NORM_2,&nrm);
83: VecDotEnd(U,a,&dot);
84: VecNormEnd(a,NORM_1,&sum);
85: VecNormEnd(a,NORM_2,&nrm);
87: if (nrm == 0.0) {
88: *zeroa = PETSC_TRUE;
89: return(0);
90: }
91: *zeroa = PETSC_FALSE;
93: /*
94: Safeguard for step sizes that are "too small"
95: */
96: #if defined(PETSC_USE_COMPLEX)
97: if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum;
98: else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum;
99: #else
100: if (dot < umin*sum && dot >= 0.0) dot = umin*sum;
101: else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum;
102: #endif
103: *h = ctx->error_rel*dot/(nrm*nrm);
104: } else {
105: *h = ctx->currenth;
106: }
107: if (*h != *h) SETERRQ3(PETSC_ERR_PLIB,"Differencing parameter is not a number sum = %G dot = %G norm = %G",sum,PetscRealPart(dot),nrm);
108: ctx->count++;
109: return(0);
110: }
114: /*
115: MatSNESMFView_DS - Prints information about this particular
116: method for computing h. Note that this does not print the general
117: information about the matrix-free method, as such info is printed
118: by the calling routine.
120: Input Parameters:
121: + ctx - the matrix free context
122: - viewer - the PETSc viewer
123: */
124: static PetscErrorCode MatSNESMFView_DS(MatSNESMFCtx ctx,PetscViewer viewer)
125: {
126: MatSNESMF_DS *hctx = (MatSNESMF_DS *)ctx->hctx;
127: PetscErrorCode ierr;
128: PetscTruth iascii;
131: /*
132: Currently this only handles the ascii file viewers, others
133: could be added, but for this type of object other viewers
134: make less sense
135: */
136: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
137: if (iascii) {
138: PetscViewerASCIIPrintf(viewer," umin=%G (minimum iterate parameter)\n",hctx->umin);
139: } else {
140: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for this SNES matrix free matrix",((PetscObject)viewer)->type_name);
141: }
142: return(0);
143: }
147: /*
148: MatSNESMFSetFromOptions_DS - Looks in the options database for
149: any options appropriate for this method.
151: Input Parameter:
152: . ctx - the matrix free context
154: */
155: static PetscErrorCode MatSNESMFSetFromOptions_DS(MatSNESMFCtx ctx)
156: {
157: PetscErrorCode ierr;
158: MatSNESMF_DS *hctx = (MatSNESMF_DS*)ctx->hctx;
161: PetscOptionsHead("Finite difference matrix free parameters");
162: PetscOptionsReal("-snes_mf_umin","umin","MatSNESMFDSSetUmin",hctx->umin,&hctx->umin,0);
163: PetscOptionsTail();
164: return(0);
165: }
169: /*
170: MatSNESMFDestroy_DS - Frees the space allocated by
171: MatSNESMFCreate_DS().
173: Input Parameter:
174: . ctx - the matrix free context
176: Notes:
177: Does not free the ctx, that is handled by the calling routine
178: */
179: static PetscErrorCode MatSNESMFDestroy_DS(MatSNESMFCtx ctx)
180: {
184: PetscFree(ctx->hctx);
185: return(0);
186: }
191: /*
192: The following two routines use the PetscObjectCompose() and PetscObjectQuery()
193: mechanism to allow the user to change the Umin parameter used in this method.
194: */
195: PetscErrorCode MatSNESMFDSSetUmin_Private(Mat mat,PetscReal umin)
196: {
197: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
198: MatSNESMF_DS *hctx;
201: if (!ctx) {
202: SETERRQ(PETSC_ERR_ARG_WRONG,"MatSNESMFDSSetUmin() attached to non-shell matrix");
203: }
204: hctx = (MatSNESMF_DS*)ctx->hctx;
205: hctx->umin = umin;
206: return(0);
207: }
212: /*@
213: MatSNESMFDSSetUmin - Sets the "umin" parameter used by the
214: PETSc routine for computing the differencing parameter, h, which is used
215: for matrix-free Jacobian-vector products.
217: Input Parameters:
218: + A - the matrix created with MatCreateSNESMF()
219: - umin - the parameter
221: Level: advanced
223: Notes:
224: See the manual page for MatCreateSNESMF() for a complete description of the
225: algorithm used to compute h.
227: .seealso: MatSNESMFSetFunctionError(), MatCreateSNESMF()
229: @*/
230: PetscErrorCode MatSNESMFDSSetUmin(Mat A,PetscReal umin)
231: {
232: PetscErrorCode ierr,(*f)(Mat,PetscReal);
236: PetscObjectQueryFunction((PetscObject)A,"MatSNESMFDSSetUmin_C",(void (**)(void))&f);
237: if (f) {
238: (*f)(A,umin);
239: }
240: return(0);
241: }
243: /*MC
244: MATSNESMF_DS - the code for compute the "h" used in the finite difference
245: matrix-free matrix vector product. This code
246: implements the strategy in Dennis and Schnabel, "Numerical Methods for Unconstrained
247: Optimization and Nonlinear Equations".
249: Options Database Keys:
250: . -snes_mf_umin <umin> see MatSNESMFDSSetUmin()
252: Level: intermediate
254: Notes: Requires 2 norms and 1 inner product, but they are computed together
255: so only one parallel collective operation is needed. See MATSNESMF_WP for a method
256: (with GMRES) that requires NO collective operations.
258: Formula used:
259: F'(u)*a = [F(u+h*a) - F(u)]/h where
260: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
261: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
262: where
263: error_rel = square root of relative error in function evaluation
264: umin = minimum iterate parameter
266: .seealso: MATMFFD, MatCreateMF(), MatCreateSNESMF(), MATSNESMF_WP, MatSNESMFDSSetUmin()
268: M*/
272: PetscErrorCode MatSNESMFCreate_DS(MatSNESMFCtx ctx)
273: {
274: MatSNESMF_DS *hctx;
275: PetscErrorCode ierr;
279: /* allocate my own private data structure */
280: PetscNew(MatSNESMF_DS,&hctx);
281: ctx->hctx = (void*)hctx;
282: /* set a default for my parameter */
283: hctx->umin = 1.e-6;
285: /* set the functions I am providing */
286: ctx->ops->compute = MatSNESMFCompute_DS;
287: ctx->ops->destroy = MatSNESMFDestroy_DS;
288: ctx->ops->view = MatSNESMFView_DS;
289: ctx->ops->setfromoptions = MatSNESMFSetFromOptions_DS;
291: PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFDSSetUmin_C",
292: "MatSNESMFDSSetUmin_Private",
293: MatSNESMFDSSetUmin_Private);
294: return(0);
295: }