Actual source code: mffddef.c
1: /*
2: Implements the DS PETSc approach for computing the h
3: parameter used with the finite difference based matrix-free
4: Jacobian-vector products.
6: To make your own: clone this file and modify for your needs.
8: Mandatory functions:
9: -------------------
10: MatMFFDCompute_ - for a given point and direction computes h
12: MatCreateMFFD _ - fills in the MatMFFD data structure
13: for this particular implementation
15: Optional functions:
16: -------------------
17: MatMFFDView_ - prints information about the parameters being used.
18: This is called when SNESView() or -snes_view is used.
20: MatMFFDSetFromOptions_ - checks the options database for options that
21: apply to this method.
23: MatMFFDDestroy_ - frees any space allocated by the routines above
25: */
27: /*
28: This include file defines the data structure MatMFFD that
29: includes information about the computation of h. It is shared by
30: all implementations that people provide
31: */
32: #include <petsc/private/matimpl.h>
33: #include <../src/mat/impls/mffd/mffdimpl.h>
35: /*
36: The method has one parameter that is used to
37: "cutoff" very small values. This is stored in a data structure
38: that is only visible to this file. If your method has no parameters
39: it can omit this, if it has several simply reorganize the data structure.
40: The data structure is "hung-off" the MatMFFD data structure in
41: the void *hctx; field.
42: */
43: typedef struct {
44: PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */
45: } MatMFFD_DS;
47: static PetscErrorCode MatMFFDCompute_DS(MatMFFD ctx, Vec U, Vec a, PetscScalar *h, PetscBool *zeroa)
48: {
49: MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
50: PetscReal nrm, sum, umin = hctx->umin;
51: PetscScalar dot;
53: PetscFunctionBegin;
54: if (!(ctx->count % ctx->recomputeperiod)) {
55: /*
56: This algorithm requires 2 norms and 1 inner product. Rather than
57: use directly the VecNorm() and VecDot() routines (and thus have
58: three separate collective operations, we use the VecxxxBegin/End() routines
59: */
60: PetscCall(VecDotBegin(U, a, &dot));
61: PetscCall(VecNormBegin(a, NORM_1, &sum));
62: PetscCall(VecNormBegin(a, NORM_2, &nrm));
63: PetscCall(VecDotEnd(U, a, &dot));
64: PetscCall(VecNormEnd(a, NORM_1, &sum));
65: PetscCall(VecNormEnd(a, NORM_2, &nrm));
67: if (nrm == 0.0) {
68: *zeroa = PETSC_TRUE;
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
71: *zeroa = PETSC_FALSE;
73: /*
74: Safeguard for step sizes that are "too small"
75: */
76: if (PetscAbsScalar(dot) < umin * sum && PetscRealPart(dot) >= 0.0) dot = umin * sum;
77: else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin * sum) dot = -umin * sum;
78: *h = ctx->error_rel * dot / (nrm * nrm);
79: PetscCheck(!PetscIsInfOrNanScalar(*h), 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);
80: } else {
81: *h = ctx->currenth;
82: }
83: ctx->count++;
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode MatMFFDView_DS(MatMFFD ctx, PetscViewer viewer)
88: {
89: MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
90: PetscBool iascii;
92: PetscFunctionBegin;
93: /*
94: Currently this only handles the ascii file viewers, others
95: could be added, but for this type of object other viewers
96: make less sense
97: */
98: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
99: if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " umin=%g (minimum iterate parameter)\n", (double)hctx->umin));
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: static PetscErrorCode MatMFFDSetFromOptions_DS(MatMFFD ctx, PetscOptionItems *PetscOptionsObject)
104: {
105: MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
107: PetscFunctionBegin;
108: PetscOptionsHeadBegin(PetscOptionsObject, "Finite difference matrix-free parameters");
109: PetscCall(PetscOptionsReal("-mat_mffd_umin", "umin", "MatMFFDDSSetUmin", hctx->umin, &hctx->umin, NULL));
110: PetscOptionsHeadEnd();
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx)
115: {
116: PetscFunctionBegin;
117: PetscCall(PetscFree(ctx->hctx));
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: /*
122: The following two routines use the PetscObjectCompose() and PetscObjectQuery()
123: mechanism to allow the user to change the Umin parameter used in this method.
124: */
125: static PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat, PetscReal umin)
126: {
127: MatMFFD ctx = NULL;
128: MatMFFD_DS *hctx;
130: PetscFunctionBegin;
131: PetscCall(MatShellGetContext(mat, &ctx));
132: PetscCheck(ctx, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatMFFDDSSetUmin() attached to non-shell matrix");
133: hctx = (MatMFFD_DS *)ctx->hctx;
134: hctx->umin = umin;
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: /*@
139: MatMFFDDSSetUmin - Sets the "umin" parameter used by the
140: PETSc routine for computing the differencing parameter, h, which is used
141: for matrix-free Jacobian-vector products for a `MATMFFD` matrix.
143: Input Parameters:
144: + A - the `MATMFFD` matrix
145: - umin - the parameter
147: Level: advanced
149: Note:
150: See the manual page for `MatCreateSNESMF()` for a complete description of the
151: algorithm used to compute h.
153: .seealso: `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
154: @*/
155: PetscErrorCode MatMFFDDSSetUmin(Mat A, PetscReal umin)
156: {
157: PetscFunctionBegin;
159: PetscTryMethod(A, "MatMFFDDSSetUmin_C", (Mat, PetscReal), (A, umin));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: /*MC
164: MATMFFD_DS - algorithm for compute the "h" used in the finite difference matrix-free matrix vector product, `MatMult()`.
166: Options Database Keys:
167: . -mat_mffd_umin <umin> - see `MatMFFDDSSetUmin()`
169: Level: intermediate
171: Notes:
172: Requires 2 norms and 1 inner product, but they are computed together
173: so only one parallel collective operation is needed. See `MATMFFD_WP` for a method
174: (with `KSPGMRES`) that requires NO collective operations.
176: Formula used:
177: F'(u)*a = [F(u+h*a) - F(u)]/h where
178: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
179: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
180: where
181: error_rel = square root of relative error in function evaluation
182: umin = minimum iterate parameter
184: Method taken from {cite}`dennis:83`
186: .seealso: `MATMFFD`, `MATMFFD_WP`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_WP`, `MatMFFDDSSetUmin()`
187: M*/
188: PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx)
189: {
190: MatMFFD_DS *hctx;
192: PetscFunctionBegin;
193: /* allocate my own private data structure */
194: PetscCall(PetscNew(&hctx));
195: ctx->hctx = (void *)hctx;
196: /* set a default for my parameter */
197: hctx->umin = 1.e-6;
199: /* set the functions I am providing */
200: ctx->ops->compute = MatMFFDCompute_DS;
201: ctx->ops->destroy = MatMFFDDestroy_DS;
202: ctx->ops->view = MatMFFDView_DS;
203: ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS;
205: PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDDSSetUmin_C", MatMFFDDSSetUmin_DS));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }