1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
2
3 /*@C
4 SNESSetObjective - Sets the objective function minimized by some of the `SNES` linesearch methods, used instead of the 2-norm of the residual in the line search
5
6 Logically Collective
7
8 Input Parameters:
9 + snes - the `SNES` context
10 . obj - objective evaluation routine; see `SNESObjectiveFn` for the calling sequence
11 - ctx - [optional] user-defined context for private data for the objective evaluation routine (may be `NULL`)
12
13 Level: intermediate
14
15 Note:
16 Some of the `SNESLineSearch` methods attempt to minimize a given objective provided by this function to determine a step length.
17
18 If not provided then this defaults to the two-norm of the function evaluation (set with `SNESSetFunction()`)
19
20 This is not used in the `SNESLINESEARCHCP` line search.
21
22 .seealso: [](ch_snes), `SNES`, `SNESLineSearch()`, `SNESGetObjective()`, `SNESComputeObjective()`, `SNESSetFunction()`, `SNESSetJacobian()`,
23 `SNESObjectiveFn`, `SNESSetObjectiveDomainError()`
24 @*/
SNESSetObjective(SNES snes,SNESObjectiveFn * obj,PetscCtx ctx)25 PetscErrorCode SNESSetObjective(SNES snes, SNESObjectiveFn *obj, PetscCtx ctx)
26 {
27 DM dm;
28
29 PetscFunctionBegin;
30 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
31 PetscCall(SNESGetDM(snes, &dm));
32 PetscCall(DMSNESSetObjective(dm, obj, ctx));
33 PetscFunctionReturn(PETSC_SUCCESS);
34 }
35
36 /*@C
37 SNESGetObjective - Returns the objective function set with `SNESSetObjective()`
38
39 Not Collective
40
41 Input Parameter:
42 . snes - the `SNES` context
43
44 Output Parameters:
45 + obj - objective evaluation routine (or `NULL`); see `SNESObjectiveFn` for the calling sequence
46 - ctx - the function context (or `NULL`)
47
48 Level: advanced
49
50 .seealso: [](ch_snes), `SNES`, `SNESSetObjective()`, `SNESGetSolution()`, `SNESObjectiveFn`
51 @*/
SNESGetObjective(SNES snes,SNESObjectiveFn ** obj,PetscCtxRt ctx)52 PetscErrorCode SNESGetObjective(SNES snes, SNESObjectiveFn **obj, PetscCtxRt ctx)
53 {
54 DM dm;
55
56 PetscFunctionBegin;
57 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
58 PetscCall(SNESGetDM(snes, &dm));
59 PetscCall(DMSNESGetObjective(dm, obj, ctx));
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
63 /*@
64 SNESComputeObjective - Computes the objective function that has been provided by `SNESSetObjective()`
65
66 Collective
67
68 Input Parameters:
69 + snes - the `SNES` context
70 - X - the state vector
71
72 Output Parameter:
73 . ob - the objective value
74
75 Level: developer
76
77 Notes:
78 `SNESComputeObjective()` is typically used within line-search routines,
79 so users would not generally call this routine themselves.
80
81 When solving for $F(x) = b$, this routine computes $objective(x) - x^T b$ where $objective(x)$ is the function provided with `SNESSetObjective()`
82
83 .seealso: [](ch_snes), `SNESLineSearch`, `SNES`, `SNESSetObjective()`, `SNESGetSolution()`
84 @*/
SNESComputeObjective(SNES snes,Vec X,PetscReal * ob)85 PetscErrorCode SNESComputeObjective(SNES snes, Vec X, PetscReal *ob)
86 {
87 DM dm;
88 DMSNES sdm;
89
90 PetscFunctionBegin;
91 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
92 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
93 PetscAssertPointer(ob, 3);
94 PetscCall(SNESGetDM(snes, &dm));
95 PetscCall(DMGetDMSNES(dm, &sdm));
96 PetscCheck(sdm->ops->computeobjective, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
97 PetscCall(PetscLogEventBegin(SNES_ObjectiveEval, snes, X, 0, 0));
98 PetscCall(sdm->ops->computeobjective(snes, X, ob, sdm->objectivectx));
99 PetscCall(PetscLogEventEnd(SNES_ObjectiveEval, snes, X, 0, 0));
100 if (snes->vec_rhs) {
101 PetscScalar dot;
102
103 PetscCall(VecDot(snes->vec_rhs, X, &dot));
104 *ob -= PetscRealPart(dot);
105 }
106 PetscFunctionReturn(PETSC_SUCCESS);
107 }
108
109 /*@C
110 SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective function
111
112 Collective
113
114 Input Parameters:
115 + snes - the `SNES` context
116 . X - the state vector
117 - ctx - the (ignored) function context
118
119 Output Parameter:
120 . F - the function value
121
122 Options Database Keys:
123 + -snes_fd_function_eps - Tolerance for including non-zero entries into the gradient, default is 1.e-6
124 - -snes_fd_function - Computes function from user provided objective function (set with `SNESSetObjective()`) with finite difference
125
126 Level: advanced
127
128 Notes:
129 This function can be used with `SNESSetFunction()` to have the nonlinear function solved for with `SNES` defined by the gradient of an objective function
130 `SNESObjectiveComputeFunctionDefaultFD()` is similar in character to `SNESComputeJacobianDefault()`.
131 Therefore, it should be used for debugging purposes only. Using it in conjunction with
132 `SNESComputeJacobianDefault()` is excessively costly and produces a Jacobian that is quite
133 noisy. This is often necessary, but should be done with care, even when debugging
134 small problems.
135
136 This uses quadratic interpolation of the objective to form each value in the function.
137
138 .seealso: [](ch_snes), `SNESSetObjective()`, `SNESSetFunction()`, `SNESComputeObjective()`, `SNESComputeJacobianDefault()`, `SNESObjectiveFn`
139 @*/
SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,PetscCtx ctx)140 PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes, Vec X, Vec F, PetscCtx ctx)
141 {
142 Vec Xh;
143 PetscInt i, N, start, end;
144 PetscReal ob, ob1, ob2, ob3, fob, dx, eps = 1e-6;
145 PetscScalar fv, xv;
146
147 PetscFunctionBegin;
148 PetscCall(VecDuplicate(X, &Xh));
149 PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing parameters", "SNES");
150 PetscCall(PetscOptionsReal("-snes_fd_function_eps", "Tolerance for nonzero entries in fd function", "None", eps, &eps, NULL));
151 PetscOptionsEnd();
152 PetscCall(VecSet(F, 0.));
153
154 PetscCall(VecNorm(X, NORM_2, &fob));
155
156 PetscCall(VecGetSize(X, &N));
157 PetscCall(VecGetOwnershipRange(X, &start, &end));
158 PetscCall(SNESComputeObjective(snes, X, &ob));
159
160 if (fob > 0.) dx = 1e-6 * fob;
161 else dx = 1e-6;
162
163 for (i = 0; i < N; i++) {
164 /* compute the 1st value */
165 PetscCall(VecCopy(X, Xh));
166 if (i >= start && i < end) {
167 xv = dx;
168 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES));
169 }
170 PetscCall(VecAssemblyBegin(Xh));
171 PetscCall(VecAssemblyEnd(Xh));
172 PetscCall(SNESComputeObjective(snes, Xh, &ob1));
173
174 /* compute the 2nd value */
175 PetscCall(VecCopy(X, Xh));
176 if (i >= start && i < end) {
177 xv = 2. * dx;
178 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES));
179 }
180 PetscCall(VecAssemblyBegin(Xh));
181 PetscCall(VecAssemblyEnd(Xh));
182 PetscCall(SNESComputeObjective(snes, Xh, &ob2));
183
184 /* compute the 3rd value */
185 PetscCall(VecCopy(X, Xh));
186 if (i >= start && i < end) {
187 xv = -dx;
188 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES));
189 }
190 PetscCall(VecAssemblyBegin(Xh));
191 PetscCall(VecAssemblyEnd(Xh));
192 PetscCall(SNESComputeObjective(snes, Xh, &ob3));
193
194 if (i >= start && i < end) {
195 /* set this entry to be the gradient of the objective */
196 fv = (-ob2 + 6. * ob1 - 3. * ob - 2. * ob3) / (6. * dx);
197 if (PetscAbsScalar(fv) > eps) {
198 PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES));
199 } else {
200 fv = 0.;
201 PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES));
202 }
203 }
204 }
205 PetscCall(VecDestroy(&Xh));
206 PetscCall(VecAssemblyBegin(F));
207 PetscCall(VecAssemblyEnd(F));
208 PetscFunctionReturn(PETSC_SUCCESS);
209 }
210