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: static PetscErrorCode MatMFFDCompute_DS(MatMFFD ctx, Vec U, Vec a, PetscScalar *h, PetscBool *zeroa)
49: {
50: MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
51: PetscReal nrm, sum, umin = hctx->umin;
52: PetscScalar dot;
54: PetscFunctionBegin;
55: if (!(ctx->count % ctx->recomputeperiod)) {
56: /*
57: This algorithm requires 2 norms and 1 inner product. Rather than
58: use directly the VecNorm() and VecDot() routines (and thus have
59: three separate collective operations, we use the VecxxxBegin/End() routines
60: */
61: PetscCall(VecDotBegin(U, a, &dot));
62: PetscCall(VecNormBegin(a, NORM_1, &sum));
63: PetscCall(VecNormBegin(a, NORM_2, &nrm));
64: PetscCall(VecDotEnd(U, a, &dot));
65: PetscCall(VecNormEnd(a, NORM_1, &sum));
66: PetscCall(VecNormEnd(a, NORM_2, &nrm));
68: if (nrm == 0.0) {
69: *zeroa = PETSC_TRUE;
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
72: *zeroa = PETSC_FALSE;
74: /*
75: Safeguard for step sizes that are "too small"
76: */
77: if (PetscAbsScalar(dot) < umin * sum && PetscRealPart(dot) >= 0.0) dot = umin * sum;
78: else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin * sum) dot = -umin * sum;
79: *h = ctx->error_rel * dot / (nrm * nrm);
80: 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);
81: } else {
82: *h = ctx->currenth;
83: }
84: ctx->count++;
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: static PetscErrorCode MatMFFDView_DS(MatMFFD ctx, PetscViewer viewer)
89: {
90: MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
91: PetscBool iascii;
93: PetscFunctionBegin;
94: /*
95: Currently this only handles the ascii file viewers, others
96: could be added, but for this type of object other viewers
97: make less sense
98: */
99: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
100: if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " umin=%g (minimum iterate parameter)\n", (double)hctx->umin));
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: static PetscErrorCode MatMFFDSetFromOptions_DS(MatMFFD ctx, PetscOptionItems *PetscOptionsObject)
105: {
106: MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
108: PetscFunctionBegin;
109: PetscOptionsHeadBegin(PetscOptionsObject, "Finite difference matrix free parameters");
110: PetscCall(PetscOptionsReal("-mat_mffd_umin", "umin", "MatMFFDDSSetUmin", hctx->umin, &hctx->umin, NULL));
111: PetscOptionsHeadEnd();
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx)
116: {
117: PetscFunctionBegin;
118: PetscCall(PetscFree(ctx->hctx));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: /*
123: The following two routines use the PetscObjectCompose() and PetscObjectQuery()
124: mechanism to allow the user to change the Umin parameter used in this method.
125: */
126: PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat, PetscReal umin)
127: {
128: MatMFFD ctx = NULL;
129: MatMFFD_DS *hctx;
131: PetscFunctionBegin;
132: PetscCall(MatShellGetContext(mat, &ctx));
133: PetscCheck(ctx, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatMFFDDSSetUmin() attached to non-shell matrix");
134: hctx = (MatMFFD_DS *)ctx->hctx;
135: hctx->umin = umin;
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: /*@
140: MatMFFDDSSetUmin - Sets the "umin" parameter used by the
141: PETSc routine for computing the differencing parameter, h, which is used
142: for matrix-free Jacobian-vector products for a `MATMFFD` matrix.
144: Input Parameters:
145: + A - the `MATMFFD` matrix
146: - umin - the parameter
148: Level: advanced
150: Note:
151: See the manual page for `MatCreateSNESMF()` for a complete description of the
152: algorithm used to compute h.
154: .seealso: `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
155: @*/
156: PetscErrorCode MatMFFDDSSetUmin(Mat A, PetscReal umin)
157: {
158: PetscFunctionBegin;
160: PetscTryMethod(A, "MatMFFDDSSetUmin_C", (Mat, PetscReal), (A, umin));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*MC
165: MATMFFD_DS - algorithm for compute the "h" used in the finite difference matrix-free matrix vector product, `MatMult()`.
167: Options Database Keys:
168: . -mat_mffd_umin <umin> - see `MatMFFDDSSetUmin()`
170: Level: intermediate
172: Notes:
173: Requires 2 norms and 1 inner product, but they are computed together
174: so only one parallel collective operation is needed. See `MATMFFD_WP` for a method
175: (with `KSPGMRES`) that requires NO collective operations.
177: Formula used:
178: F'(u)*a = [F(u+h*a) - F(u)]/h where
179: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
180: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
181: where
182: error_rel = square root of relative error in function evaluation
183: umin = minimum iterate parameter
185: References:
186: . * - Dennis and Schnabel, "Numerical Methods for Unconstrained Optimization and Nonlinear Equations"
188: .seealso: `MATMFFD`, `MATMFFD_WP`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_WP`, `MatMFFDDSSetUmin()`
189: M*/
190: PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx)
191: {
192: MatMFFD_DS *hctx;
194: PetscFunctionBegin;
195: /* allocate my own private data structure */
196: PetscCall(PetscNew(&hctx));
197: ctx->hctx = (void *)hctx;
198: /* set a default for my parameter */
199: hctx->umin = 1.e-6;
201: /* set the functions I am providing */
202: ctx->ops->compute = MatMFFDCompute_DS;
203: ctx->ops->destroy = MatMFFDDestroy_DS;
204: ctx->ops->view = MatMFFDView_DS;
205: ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS;
207: PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDDSSetUmin_C", MatMFFDDSSetUmin_DS));
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }