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