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 @*/ 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 @*/ 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 @*/ 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 @*/ 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