Actual source code: snesob.c
1: #include <petsc/private/snesimpl.h>
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
6: Logically Collective
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`)
13: Level: intermediate
15: Note:
16: Some of the `SNESLineSearch` methods attempt to minimize a given objective provided by this function to determine a step length.
18: If not provided then this defaults to the two-norm of the function evaluation (set with `SNESSetFunction()`)
20: This is not used in the `SNESLINESEARCHCP` line search.
22: .seealso: [](ch_snes), `SNES`, `SNESLineSearch()`, `SNESGetObjective()`, `SNESComputeObjective()`, `SNESSetFunction()`, `SNESSetJacobian()`,
23: `SNESObjectiveFn`
24: @*/
25: PetscErrorCode SNESSetObjective(SNES snes, SNESObjectiveFn *obj, void *ctx)
26: {
27: DM dm;
29: PetscFunctionBegin;
31: PetscCall(SNESGetDM(snes, &dm));
32: PetscCall(DMSNESSetObjective(dm, obj, ctx));
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: /*@C
37: SNESGetObjective - Returns the objective function set with `SNESSetObjective()`
39: Not Collective
41: Input Parameter:
42: . snes - the `SNES` context
44: Output Parameters:
45: + obj - objective evaluation routine (or `NULL`); see `SNESObjectiveFn` for the calling sequence
46: - ctx - the function context (or `NULL`)
48: Level: advanced
50: .seealso: [](ch_snes), `SNES`, `SNESSetObjective()`, `SNESGetSolution()`, `SNESObjectiveFn`
51: @*/
52: PetscErrorCode SNESGetObjective(SNES snes, SNESObjectiveFn **obj, void **ctx)
53: {
54: DM dm;
56: PetscFunctionBegin;
58: PetscCall(SNESGetDM(snes, &dm));
59: PetscCall(DMSNESGetObjective(dm, obj, ctx));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: /*@
64: SNESComputeObjective - Computes the objective function that has been provided by `SNESSetObjective()`
66: Collective
68: Input Parameters:
69: + snes - the `SNES` context
70: - X - the state vector
72: Output Parameter:
73: . ob - the objective value
75: Level: developer
77: Notes:
78: `SNESComputeObjective()` is typically used within line-search routines,
79: so users would not generally call this routine themselves.
81: When solving for $F(x) = b$, this routine computes $objective(x) - x^T b$ where $objective(x)$ is the function provided with `SNESSetObjective()`
83: .seealso: [](ch_snes), `SNESLineSearch`, `SNES`, `SNESSetObjective()`, `SNESGetSolution()`
84: @*/
85: PetscErrorCode SNESComputeObjective(SNES snes, Vec X, PetscReal *ob)
86: {
87: DM dm;
88: DMSNES sdm;
90: PetscFunctionBegin;
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;
103: PetscCall(VecDot(snes->vec_rhs, X, &dot));
104: *ob -= PetscRealPart(dot);
105: }
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: /*@C
110: SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective function
112: Collective
114: Input Parameters:
115: + snes - the `SNES` context
116: . X - the state vector
117: - ctx - the (ignored) function context
119: Output Parameter:
120: . F - the function value
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
126: Level: advanced
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.
136: This uses quadratic interpolation of the objective to form each value in the function.
138: .seealso: [](ch_snes), `SNESSetObjective()`, `SNESSetFunction()`, `SNESComputeObjective()`, `SNESComputeJacobianDefault()`, `SNESObjectiveFn`
139: @*/
140: PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes, Vec X, Vec F, void *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;
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.));
154: PetscCall(VecNorm(X, NORM_2, &fob));
156: PetscCall(VecGetSize(X, &N));
157: PetscCall(VecGetOwnershipRange(X, &start, &end));
158: PetscCall(SNESComputeObjective(snes, X, &ob));
160: if (fob > 0.) dx = 1e-6 * fob;
161: else dx = 1e-6;
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));
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));
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));
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: }