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.
5
6 To make your own: clone this file and modify for your needs.
7
8 Mandatory functions:
9 -------------------
10 MatMFFDCompute_ - for a given point and direction computes h
11
12 MatCreateMFFD _ - fills in the MatMFFD data structure
13 for this particular implementation
14
15 Optional functions:
16 -------------------
17 MatMFFDView_ - prints information about the parameters being used.
18 This is called when SNESView() or -snes_view is used.
19
20 MatMFFDSetFromOptions_ - checks the options database for options that
21 apply to this method.
22
23 MatMFFDDestroy_ - frees any space allocated by the routines above
24
25 */
26
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> /*I "petscmat.h" I*/
34
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;
46
MatMFFDCompute_DS(MatMFFD ctx,Vec U,Vec a,PetscScalar * h,PetscBool * zeroa)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;
52
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));
66
67 if (nrm == 0.0) {
68 *zeroa = PETSC_TRUE;
69 PetscFunctionReturn(PETSC_SUCCESS);
70 }
71 *zeroa = PETSC_FALSE;
72
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 = dot * ctx->error_rel / (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 }
86
MatMFFDView_DS(MatMFFD ctx,PetscViewer viewer)87 static PetscErrorCode MatMFFDView_DS(MatMFFD ctx, PetscViewer viewer)
88 {
89 MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
90 PetscBool isascii;
91
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, &isascii));
99 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, " umin=%g (minimum iterate parameter)\n", (double)hctx->umin));
100 PetscFunctionReturn(PETSC_SUCCESS);
101 }
102
MatMFFDSetFromOptions_DS(MatMFFD ctx,PetscOptionItems PetscOptionsObject)103 static PetscErrorCode MatMFFDSetFromOptions_DS(MatMFFD ctx, PetscOptionItems PetscOptionsObject)
104 {
105 MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
106
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 }
113
MatMFFDDestroy_DS(MatMFFD ctx)114 static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx)
115 {
116 PetscFunctionBegin;
117 PetscCall(PetscFree(ctx->hctx));
118 PetscFunctionReturn(PETSC_SUCCESS);
119 }
120
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 */
MatMFFDDSSetUmin_DS(Mat mat,PetscReal umin)125 static PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat, PetscReal umin)
126 {
127 MatMFFD ctx = NULL;
128 MatMFFD_DS *hctx;
129
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 }
137
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.
142
143 Input Parameters:
144 + A - the `MATMFFD` matrix
145 - umin - the parameter
146
147 Level: advanced
148
149 Note:
150 See the manual page for `MatCreateSNESMF()` for a complete description of the
151 algorithm used to compute h.
152
153 .seealso: `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
154 @*/
MatMFFDDSSetUmin(Mat A,PetscReal umin)155 PetscErrorCode MatMFFDDSSetUmin(Mat A, PetscReal umin)
156 {
157 PetscFunctionBegin;
158 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
159 PetscTryMethod(A, "MatMFFDDSSetUmin_C", (Mat, PetscReal), (A, umin));
160 PetscFunctionReturn(PETSC_SUCCESS);
161 }
162
163 /*MC
164 MATMFFD_DS - algorithm for compute the "h" used in the finite difference matrix-free matrix vector product, `MatMult()`.
165
166 Options Database Keys:
167 . -mat_mffd_umin <umin> - see `MatMFFDDSSetUmin()`
168
169 Level: intermediate
170
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.
175
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
183
184 Method taken from {cite}`dennis:83`
185
186 .seealso: `MATMFFD`, `MATMFFD_WP`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_WP`, `MatMFFDDSSetUmin()`
187 M*/
MatCreateMFFD_DS(MatMFFD ctx)188 PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx)
189 {
190 MatMFFD_DS *hctx;
191
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;
198
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;
204
205 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDDSSetUmin_C", MatMFFDDSSetUmin_DS));
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208