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