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 PetscErrorCode ierr; 147 PetscInt i,N,start,end; 148 PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6; 149 PetscScalar fv,xv; 150 151 PetscFunctionBegin; 152 PetscCall(VecDuplicate(X,&Xh)); 153 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing parameters","SNES");PetscCall(ierr); 154 PetscCall(PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL)); 155 ierr = PetscOptionsEnd();PetscCall(ierr); 156 PetscCall(VecSet(F,0.)); 157 158 PetscCall(VecNorm(X,NORM_2,&fob)); 159 160 PetscCall(VecGetSize(X,&N)); 161 PetscCall(VecGetOwnershipRange(X,&start,&end)); 162 PetscCall(SNESComputeObjective(snes,X,&ob)); 163 164 if (fob > 0.) dx =1e-6*fob; 165 else dx = 1e-6; 166 167 for (i=0; i<N; i++) { 168 /* compute the 1st value */ 169 PetscCall(VecCopy(X,Xh)); 170 if (i>= start && i<end) { 171 xv = dx; 172 PetscCall(VecSetValues(Xh,1,&i,&xv,ADD_VALUES)); 173 } 174 PetscCall(VecAssemblyBegin(Xh)); 175 PetscCall(VecAssemblyEnd(Xh)); 176 PetscCall(SNESComputeObjective(snes,Xh,&ob1)); 177 178 /* compute the 2nd value */ 179 PetscCall(VecCopy(X,Xh)); 180 if (i>= start && i<end) { 181 xv = 2.*dx; 182 PetscCall(VecSetValues(Xh,1,&i,&xv,ADD_VALUES)); 183 } 184 PetscCall(VecAssemblyBegin(Xh)); 185 PetscCall(VecAssemblyEnd(Xh)); 186 PetscCall(SNESComputeObjective(snes,Xh,&ob2)); 187 188 /* compute the 3rd value */ 189 PetscCall(VecCopy(X,Xh)); 190 if (i>= start && i<end) { 191 xv = -dx; 192 PetscCall(VecSetValues(Xh,1,&i,&xv,ADD_VALUES)); 193 } 194 PetscCall(VecAssemblyBegin(Xh)); 195 PetscCall(VecAssemblyEnd(Xh)); 196 PetscCall(SNESComputeObjective(snes,Xh,&ob3)); 197 198 if (i >= start && i<end) { 199 /* set this entry to be the gradient of the objective */ 200 fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx); 201 if (PetscAbsScalar(fv) > eps) { 202 PetscCall(VecSetValues(F,1,&i,&fv,INSERT_VALUES)); 203 } else { 204 fv = 0.; 205 PetscCall(VecSetValues(F,1,&i,&fv,INSERT_VALUES)); 206 } 207 } 208 } 209 210 PetscCall(VecDestroy(&Xh)); 211 212 PetscCall(VecAssemblyBegin(F)); 213 PetscCall(VecAssemblyEnd(F)); 214 PetscFunctionReturn(0); 215 } 216