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