Actual source code: mffddef.c
2: /*
3: Implements the DS PETSc approach for computing the h
4: parameter used with the finite difference based matrix-free
5: Jacobian-vector products.
7: To make your own: clone this file and modify for your needs.
9: Mandatory functions:
10: -------------------
11: MatMFFDCompute_ - for a given point and direction computes h
13: MatCreateMFFD _ - fills in the MatMFFD data structure
14: for this particular implementation
16: Optional functions:
17: -------------------
18: MatMFFDView_ - prints information about the parameters being used.
19: This is called when SNESView() or -snes_view is used.
21: MatMFFDSetFromOptions_ - checks the options database for options that
22: apply to this method.
24: MatMFFDDestroy_ - frees any space allocated by the routines above
26: */
28: /*
29: This include file defines the data structure MatMFFD that
30: includes information about the computation of h. It is shared by
31: all implementations that people provide
32: */
33: #include <petsc/private/matimpl.h>
34: #include <../src/mat/impls/mffd/mffdimpl.h>
36: /*
37: The method has one parameter that is used to
38: "cutoff" very small values. This is stored in a data structure
39: that is only visible to this file. If your method has no parameters
40: it can omit this, if it has several simply reorganize the data structure.
41: The data structure is "hung-off" the MatMFFD data structure in
42: the void *hctx; field.
43: */
44: typedef struct {
45: PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */
46: } MatMFFD_DS;
48: /*
49: MatMFFDCompute_DS - Standard PETSc code for computing the
50: differencing parameter (h) for use with matrix-free finite differences.
52: Input Parameters:
53: + ctx - the matrix free context
54: . U - the location at which you want the Jacobian
55: - a - the direction you want the derivative
57: Output Parameter:
58: . h - the scale computed
60: */
61: static PetscErrorCode MatMFFDCompute_DS(MatMFFD ctx,Vec U,Vec a,PetscScalar *h,PetscBool *zeroa)
62: {
63: MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx;
64: PetscReal nrm,sum,umin = hctx->umin;
65: PetscScalar dot;
67: if (!(ctx->count % ctx->recomputeperiod)) {
68: /*
69: This algorithm requires 2 norms and 1 inner product. Rather than
70: use directly the VecNorm() and VecDot() routines (and thus have
71: three separate collective operations, we use the VecxxxBegin/End() routines
72: */
73: VecDotBegin(U,a,&dot);
74: VecNormBegin(a,NORM_1,&sum);
75: VecNormBegin(a,NORM_2,&nrm);
76: VecDotEnd(U,a,&dot);
77: VecNormEnd(a,NORM_1,&sum);
78: VecNormEnd(a,NORM_2,&nrm);
80: if (nrm == 0.0) {
81: *zeroa = PETSC_TRUE;
82: return 0;
83: }
84: *zeroa = PETSC_FALSE;
86: /*
87: Safeguard for step sizes that are "too small"
88: */
89: if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum;
90: else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum;
91: *h = ctx->error_rel*dot/(nrm*nrm);
93: } else {
94: *h = ctx->currenth;
95: }
96: ctx->count++;
97: return 0;
98: }
100: /*
101: MatMFFDView_DS - Prints information about this particular
102: method for computing h. Note that this does not print the general
103: information about the matrix-free method, as such info is printed
104: by the calling routine.
106: Input Parameters:
107: + ctx - the matrix free context
108: - viewer - the PETSc viewer
109: */
110: static PetscErrorCode MatMFFDView_DS(MatMFFD ctx,PetscViewer viewer)
111: {
112: MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx;
113: PetscBool iascii;
115: /*
116: Currently this only handles the ascii file viewers, others
117: could be added, but for this type of object other viewers
118: make less sense
119: */
120: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
121: if (iascii) {
122: PetscViewerASCIIPrintf(viewer," umin=%g (minimum iterate parameter)\n",(double)hctx->umin);
123: }
124: return 0;
125: }
127: /*
128: MatMFFDSetFromOptions_DS - Looks in the options database for
129: any options appropriate for this method.
131: Input Parameter:
132: . ctx - the matrix free context
134: */
135: static PetscErrorCode MatMFFDSetFromOptions_DS(PetscOptionItems *PetscOptionsObject,MatMFFD ctx)
136: {
137: MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx;
139: PetscOptionsHead(PetscOptionsObject,"Finite difference matrix free parameters");
140: PetscOptionsReal("-mat_mffd_umin","umin","MatMFFDDSSetUmin",hctx->umin,&hctx->umin,NULL);
141: PetscOptionsTail();
142: return 0;
143: }
145: /*
146: MatMFFDDestroy_DS - Frees the space allocated by
147: MatCreateMFFD_DS().
149: Input Parameter:
150: . ctx - the matrix free context
152: Notes:
153: Does not free the ctx, that is handled by the calling routine
154: */
155: static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx)
156: {
157: PetscFree(ctx->hctx);
158: return 0;
159: }
161: /*
162: The following two routines use the PetscObjectCompose() and PetscObjectQuery()
163: mechanism to allow the user to change the Umin parameter used in this method.
164: */
165: PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat,PetscReal umin)
166: {
167: MatMFFD ctx=NULL;
168: MatMFFD_DS *hctx;
170: MatShellGetContext(mat,&ctx);
172: hctx = (MatMFFD_DS*)ctx->hctx;
173: hctx->umin = umin;
174: return 0;
175: }
177: /*@
178: MatMFFDDSSetUmin - Sets the "umin" parameter used by the
179: PETSc routine for computing the differencing parameter, h, which is used
180: for matrix-free Jacobian-vector products.
182: Input Parameters:
183: + A - the matrix created with MatCreateSNESMF()
184: - umin - the parameter
186: Level: advanced
188: Notes:
189: See the manual page for MatCreateSNESMF() for a complete description of the
190: algorithm used to compute h.
192: .seealso: MatMFFDSetFunctionError(), MatCreateSNESMF()
194: @*/
195: PetscErrorCode MatMFFDDSSetUmin(Mat A,PetscReal umin)
196: {
198: PetscTryMethod(A,"MatMFFDDSSetUmin_C",(Mat,PetscReal),(A,umin));
199: return 0;
200: }
202: /*MC
203: MATMFFD_DS - the code for compute the "h" used in the finite difference
204: matrix-free matrix vector product. This code
205: implements the strategy in Dennis and Schnabel, "Numerical Methods for Unconstrained
206: Optimization and Nonlinear Equations".
208: Options Database Keys:
209: . -mat_mffd_umin <umin> - see MatMFFDDSSetUmin()
211: Level: intermediate
213: Notes:
214: Requires 2 norms and 1 inner product, but they are computed together
215: so only one parallel collective operation is needed. See MATMFFD_WP for a method
216: (with GMRES) that requires NO collective operations.
218: Formula used:
219: F'(u)*a = [F(u+h*a) - F(u)]/h where
220: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
221: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
222: where
223: error_rel = square root of relative error in function evaluation
224: umin = minimum iterate parameter
226: .seealso: MATMFFD, MatCreateMFFD(), MatCreateSNESMF(), MATMFFD_WP, MatMFFDDSSetUmin()
228: M*/
229: PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx)
230: {
231: MatMFFD_DS *hctx;
233: /* allocate my own private data structure */
234: PetscNewLog(ctx,&hctx);
235: ctx->hctx = (void*)hctx;
236: /* set a default for my parameter */
237: hctx->umin = 1.e-6;
239: /* set the functions I am providing */
240: ctx->ops->compute = MatMFFDCompute_DS;
241: ctx->ops->destroy = MatMFFDDestroy_DS;
242: ctx->ops->view = MatMFFDView_DS;
243: ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS;
245: PetscObjectComposeFunction((PetscObject)ctx->mat,"MatMFFDDSSetUmin_C",MatMFFDDSSetUmin_DS);
246: return 0;
247: }