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: 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> /*I "petscmat.h" I*/
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;
51: /*
52: MatMFFDCompute_DS - Standard PETSc code for computing the
53: differencing paramter (h) for use with matrix-free finite differences.
55: Input Parameters:
56: + ctx - the matrix free context
57: . U - the location at which you want the Jacobian
58: - a - the direction you want the derivative
60: 61: Output Parameter:
62: . h - the scale computed
64: */
65: static PetscErrorCode MatMFFDCompute_DS(MatMFFD ctx,Vec U,Vec a,PetscScalar *h,PetscBool *zeroa) 66: {
67: MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx;
68: PetscReal nrm,sum,umin = hctx->umin;
69: PetscScalar dot;
70: PetscErrorCode ierr;
73: if (!(ctx->count % ctx->recomputeperiod)) {
74: /*
75: This algorithm requires 2 norms and 1 inner product. Rather than
76: use directly the VecNorm() and VecDot() routines (and thus have
77: three separate collective operations, we use the VecxxxBegin/End() routines
78: */
79: VecDotBegin(U,a,&dot);
80: VecNormBegin(a,NORM_1,&sum);
81: VecNormBegin(a,NORM_2,&nrm);
82: VecDotEnd(U,a,&dot);
83: VecNormEnd(a,NORM_1,&sum);
84: VecNormEnd(a,NORM_2,&nrm);
86: if (nrm == 0.0) {
87: *zeroa = PETSC_TRUE;
88: return(0);
89: }
90: *zeroa = PETSC_FALSE;
92: /*
93: Safeguard for step sizes that are "too small"
94: */
95: #if defined(PETSC_USE_COMPLEX)
96: if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum;
97: else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum;
98: #else
99: if (dot < umin*sum && dot >= 0.0) dot = umin*sum;
100: else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum;
101: #endif
102: *h = ctx->error_rel*dot/(nrm*nrm);
103: } else {
104: *h = ctx->currenth;
105: }
106: if (*h != *h) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Differencing parameter is not a number sum = %G dot = %G norm = %G",sum,PetscRealPart(dot),nrm);
107: ctx->count++;
108: return(0);
109: }
113: /*
114: MatMFFDView_DS - Prints information about this particular
115: method for computing h. Note that this does not print the general
116: information about the matrix-free method, as such info is printed
117: by the calling routine.
119: Input Parameters:
120: + ctx - the matrix free context
121: - viewer - the PETSc viewer
122: */
123: static PetscErrorCode MatMFFDView_DS(MatMFFD ctx,PetscViewer viewer)124: {
125: MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
126: PetscErrorCode ierr;
127: PetscBool iascii;
130: /*
131: Currently this only handles the ascii file viewers, others
132: could be added, but for this type of object other viewers
133: make less sense
134: */
135: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
136: if (iascii) {
137: PetscViewerASCIIPrintf(viewer," umin=%G (minimum iterate parameter)\n",hctx->umin);
138: } else {
139: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for this SNES matrix free matrix",((PetscObject)viewer)->type_name);
140: }
141: return(0);
142: }
146: /*
147: MatMFFDSetFromOptions_DS - Looks in the options database for
148: any options appropriate for this method.
150: Input Parameter:
151: . ctx - the matrix free context
153: */
154: static PetscErrorCode MatMFFDSetFromOptions_DS(MatMFFD ctx)155: {
156: PetscErrorCode ierr;
157: MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx;
160: PetscOptionsHead("Finite difference matrix free parameters");
161: PetscOptionsReal("-mat_mffd_umin","umin","MatMFFDDSSetUmin",hctx->umin,&hctx->umin,0);
162: PetscOptionsTail();
163: return(0);
164: }
168: /*
169: MatMFFDDestroy_DS - Frees the space allocated by
170: MatCreateMFFD_DS().
172: Input Parameter:
173: . ctx - the matrix free context
175: Notes:
176: Does not free the ctx, that is handled by the calling routine
177: */
178: static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx)179: {
183: PetscFree(ctx->hctx);
184: return(0);
185: }
187: EXTERN_C_BEGIN
190: /*
191: The following two routines use the PetscObjectCompose() and PetscObjectQuery()
192: mechanism to allow the user to change the Umin parameter used in this method.
193: */
194: PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat,PetscReal umin)195: {
196: MatMFFD ctx = (MatMFFD)mat->data;
197: MatMFFD_DS *hctx;
200: if (!ctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatMFFDDSSetUmin() attached to non-shell matrix");
201: hctx = (MatMFFD_DS*)ctx->hctx;
202: hctx->umin = umin;
203: return(0);
204: }
205: EXTERN_C_END
209: /*@
210: MatMFFDDSSetUmin - Sets the "umin" parameter used by the
211: PETSc routine for computing the differencing parameter, h, which is used
212: for matrix-free Jacobian-vector products.
214: Input Parameters:
215: + A - the matrix created with MatCreateSNESMF()
216: - umin - the parameter
218: Level: advanced
220: Notes:
221: See the manual page for MatCreateSNESMF() for a complete description of the
222: algorithm used to compute h.
224: .seealso: MatMFFDSetFunctionError(), MatCreateSNESMF()
226: @*/
227: PetscErrorCodeMatMFFDDSSetUmin(Mat A,PetscReal umin)228: {
233: PetscTryMethod(A,"MatMFFDDSSetUmin_C",(Mat,PetscReal),(A,umin));
234: return(0);
235: }
237: /*MC
238: MATMFFD_DS - the code for compute the "h" used in the finite difference
239: matrix-free matrix vector product. This code
240: implements the strategy in Dennis and Schnabel, "Numerical Methods for Unconstrained
241: Optimization and Nonlinear Equations".
243: Options Database Keys:
244: . -mat_mffd_umin <umin> see MatMFFDDSSetUmin()
246: Level: intermediate
248: Notes: Requires 2 norms and 1 inner product, but they are computed together
249: so only one parallel collective operation is needed. See MATMFFD_WP for a method
250: (with GMRES) that requires NO collective operations.
252: Formula used:
253: F'(u)*a = [F(u+h*a) - F(u)]/h where
254: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
255: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
256: where
257: error_rel = square root of relative error in function evaluation
258: umin = minimum iterate parameter
260: .seealso: MATMFFD, MatCreateMFFD(), MatCreateSNESMF(), MATMFFD_WP, MatMFFDDSSetUmin()
262: M*/
263: EXTERN_C_BEGIN
266: PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx)267: {
268: MatMFFD_DS *hctx;
269: PetscErrorCode ierr;
273: /* allocate my own private data structure */
274: PetscNewLog(ctx,MatMFFD_DS,&hctx);
275: ctx->hctx = (void*)hctx;
276: /* set a default for my parameter */
277: hctx->umin = 1.e-6;
279: /* set the functions I am providing */
280: ctx->ops->compute = MatMFFDCompute_DS;
281: ctx->ops->destroy = MatMFFDDestroy_DS;
282: ctx->ops->view = MatMFFDView_DS;
283: ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS;
285: PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatMFFDDSSetUmin_C",
286: "MatMFFDDSSetUmin_DS",
287: MatMFFDDSSetUmin_DS);
288: return(0);
289: }
290: EXTERN_C_END