Actual source code: mffddef.c
petsc-3.9.4 2018-09-11
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
17: Optional functions:
18: -------------------
19: MatMFFDView_ - prints information about the parameters being used.
20: This is called when SNESView() or -snes_view is used.
22: MatMFFDSetFromOptions_ - checks the options database for options that
23: apply to this method.
25: MatMFFDDestroy_ - frees any space allocated by the routines above
27: */
29: /*
30: This include file defines the data structure MatMFFD that
31: includes information about the computation of h. It is shared by
32: all implementations that people provide
33: */
34: #include <petsc/private/matimpl.h>
35: #include <../src/mat/impls/mffd/mffdimpl.h>
37: /*
38: The method has one parameter that is used to
39: "cutoff" very small values. This is stored in a data structure
40: that is only visible to this file. If your method has no parameters
41: it can omit this, if it has several simply reorganize the data structure.
42: The data structure is "hung-off" the MatMFFD data structure in
43: the void *hctx; field.
44: */
45: typedef struct {
46: PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */
47: } MatMFFD_DS;
49: /*
50: MatMFFDCompute_DS - Standard PETSc code for computing the
51: differencing paramter (h) for use with matrix-free finite differences.
53: Input Parameters:
54: + ctx - the matrix free context
55: . U - the location at which you want the Jacobian
56: - a - the direction you want the derivative
59: Output Parameter:
60: . h - the scale computed
62: */
63: static PetscErrorCode MatMFFDCompute_DS(MatMFFD ctx,Vec U,Vec a,PetscScalar *h,PetscBool *zeroa)
64: {
65: MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx;
66: PetscReal nrm,sum,umin = hctx->umin;
67: PetscScalar dot;
71: if (!(ctx->count % ctx->recomputeperiod)) {
72: /*
73: This algorithm requires 2 norms and 1 inner product. Rather than
74: use directly the VecNorm() and VecDot() routines (and thus have
75: three separate collective operations, we use the VecxxxBegin/End() routines
76: */
77: VecDotBegin(U,a,&dot);
78: VecNormBegin(a,NORM_1,&sum);
79: VecNormBegin(a,NORM_2,&nrm);
80: VecDotEnd(U,a,&dot);
81: VecNormEnd(a,NORM_1,&sum);
82: VecNormEnd(a,NORM_2,&nrm);
84: if (nrm == 0.0) {
85: *zeroa = PETSC_TRUE;
86: return(0);
87: }
88: *zeroa = PETSC_FALSE;
90: /*
91: Safeguard for step sizes that are "too small"
92: */
93: if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum;
94: else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum;
95: *h = ctx->error_rel*dot/(nrm*nrm);
96: if (PetscIsInfOrNanScalar(*h)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Differencing parameter is not a number sum = %g dot = %g norm = %g",(double)sum,(double)PetscRealPart(dot),(double)nrm);
97: } else {
98: *h = ctx->currenth;
99: }
100: ctx->count++;
101: return(0);
102: }
104: /*
105: MatMFFDView_DS - Prints information about this particular
106: method for computing h. Note that this does not print the general
107: information about the matrix-free method, as such info is printed
108: by the calling routine.
110: Input Parameters:
111: + ctx - the matrix free context
112: - viewer - the PETSc viewer
113: */
114: static PetscErrorCode MatMFFDView_DS(MatMFFD ctx,PetscViewer viewer)
115: {
116: MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx;
118: PetscBool iascii;
121: /*
122: Currently this only handles the ascii file viewers, others
123: could be added, but for this type of object other viewers
124: make less sense
125: */
126: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
127: if (iascii) {
128: PetscViewerASCIIPrintf(viewer," umin=%g (minimum iterate parameter)\n",(double)hctx->umin);
129: }
130: return(0);
131: }
133: /*
134: MatMFFDSetFromOptions_DS - Looks in the options database for
135: any options appropriate for this method.
137: Input Parameter:
138: . ctx - the matrix free context
140: */
141: static PetscErrorCode MatMFFDSetFromOptions_DS(PetscOptionItems *PetscOptionsObject,MatMFFD ctx)
142: {
144: MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx;
147: PetscOptionsHead(PetscOptionsObject,"Finite difference matrix free parameters");
148: PetscOptionsReal("-mat_mffd_umin","umin","MatMFFDDSSetUmin",hctx->umin,&hctx->umin,0);
149: PetscOptionsTail();
150: return(0);
151: }
153: /*
154: MatMFFDDestroy_DS - Frees the space allocated by
155: MatCreateMFFD_DS().
157: Input Parameter:
158: . ctx - the matrix free context
160: Notes:
161: Does not free the ctx, that is handled by the calling routine
162: */
163: static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx)
164: {
168: PetscFree(ctx->hctx);
169: return(0);
170: }
172: /*
173: The following two routines use the PetscObjectCompose() and PetscObjectQuery()
174: mechanism to allow the user to change the Umin parameter used in this method.
175: */
176: PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat,PetscReal umin)
177: {
178: MatMFFD ctx = (MatMFFD)mat->data;
179: MatMFFD_DS *hctx;
182: if (!ctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatMFFDDSSetUmin() attached to non-shell matrix");
183: hctx = (MatMFFD_DS*)ctx->hctx;
184: hctx->umin = umin;
185: return(0);
186: }
188: /*@
189: MatMFFDDSSetUmin - Sets the "umin" parameter used by the
190: PETSc routine for computing the differencing parameter, h, which is used
191: for matrix-free Jacobian-vector products.
193: Input Parameters:
194: + A - the matrix created with MatCreateSNESMF()
195: - umin - the parameter
197: Level: advanced
199: Notes:
200: See the manual page for MatCreateSNESMF() for a complete description of the
201: algorithm used to compute h.
203: .seealso: MatMFFDSetFunctionError(), MatCreateSNESMF()
205: @*/
206: PetscErrorCode MatMFFDDSSetUmin(Mat A,PetscReal umin)
207: {
212: PetscTryMethod(A,"MatMFFDDSSetUmin_C",(Mat,PetscReal),(A,umin));
213: return(0);
214: }
216: /*MC
217: MATMFFD_DS - the code for compute the "h" used in the finite difference
218: matrix-free matrix vector product. This code
219: implements the strategy in Dennis and Schnabel, "Numerical Methods for Unconstrained
220: Optimization and Nonlinear Equations".
222: Options Database Keys:
223: . -mat_mffd_umin <umin> see MatMFFDDSSetUmin()
225: Level: intermediate
227: Notes: Requires 2 norms and 1 inner product, but they are computed together
228: so only one parallel collective operation is needed. See MATMFFD_WP for a method
229: (with GMRES) that requires NO collective operations.
231: Formula used:
232: F'(u)*a = [F(u+h*a) - F(u)]/h where
233: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
234: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
235: where
236: error_rel = square root of relative error in function evaluation
237: umin = minimum iterate parameter
239: .seealso: MATMFFD, MatCreateMFFD(), MatCreateSNESMF(), MATMFFD_WP, MatMFFDDSSetUmin()
241: M*/
242: PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx)
243: {
244: MatMFFD_DS *hctx;
248: /* allocate my own private data structure */
249: PetscNewLog(ctx,&hctx);
250: ctx->hctx = (void*)hctx;
251: /* set a default for my parameter */
252: hctx->umin = 1.e-6;
254: /* set the functions I am providing */
255: ctx->ops->compute = MatMFFDCompute_DS;
256: ctx->ops->destroy = MatMFFDDestroy_DS;
257: ctx->ops->view = MatMFFDView_DS;
258: ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS;
260: PetscObjectComposeFunction((PetscObject)ctx->mat,"MatMFFDDSSetUmin_C",MatMFFDDSSetUmin_DS);
261: return(0);
262: }