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 CHKERRQ(SNESGetDM(snes,&dm)); 47 CHKERRQ(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 CHKERRQ(SNESGetDM(snes,&dm)); 74 CHKERRQ(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 CHKERRQ(SNESGetDM(snes,&dm)); 104 CHKERRQ(DMGetDMSNES(dm,&sdm)); 105 if (sdm->ops->computeobjective) { 106 CHKERRQ(PetscLogEventBegin(SNES_ObjectiveEval,snes,X,0,0)); 107 CHKERRQ((sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx)); 108 CHKERRQ(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 CHKERRQ(VecDuplicate(X,&Xh)); 153 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing parameters","SNES");CHKERRQ(ierr); 154 CHKERRQ(PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL)); 155 ierr = PetscOptionsEnd();CHKERRQ(ierr); 156 CHKERRQ(VecSet(F,0.)); 157 158 CHKERRQ(VecNorm(X,NORM_2,&fob)); 159 160 CHKERRQ(VecGetSize(X,&N)); 161 CHKERRQ(VecGetOwnershipRange(X,&start,&end)); 162 CHKERRQ(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 CHKERRQ(VecCopy(X,Xh)); 170 if (i>= start && i<end) { 171 xv = dx; 172 CHKERRQ(VecSetValues(Xh,1,&i,&xv,ADD_VALUES)); 173 } 174 CHKERRQ(VecAssemblyBegin(Xh)); 175 CHKERRQ(VecAssemblyEnd(Xh)); 176 CHKERRQ(SNESComputeObjective(snes,Xh,&ob1)); 177 178 /* compute the 2nd value */ 179 CHKERRQ(VecCopy(X,Xh)); 180 if (i>= start && i<end) { 181 xv = 2.*dx; 182 CHKERRQ(VecSetValues(Xh,1,&i,&xv,ADD_VALUES)); 183 } 184 CHKERRQ(VecAssemblyBegin(Xh)); 185 CHKERRQ(VecAssemblyEnd(Xh)); 186 CHKERRQ(SNESComputeObjective(snes,Xh,&ob2)); 187 188 /* compute the 3rd value */ 189 CHKERRQ(VecCopy(X,Xh)); 190 if (i>= start && i<end) { 191 xv = -dx; 192 CHKERRQ(VecSetValues(Xh,1,&i,&xv,ADD_VALUES)); 193 } 194 CHKERRQ(VecAssemblyBegin(Xh)); 195 CHKERRQ(VecAssemblyEnd(Xh)); 196 CHKERRQ(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 CHKERRQ(VecSetValues(F,1,&i,&fv,INSERT_VALUES)); 203 } else { 204 fv = 0.; 205 CHKERRQ(VecSetValues(F,1,&i,&fv,INSERT_VALUES)); 206 } 207 } 208 } 209 210 CHKERRQ(VecDestroy(&Xh)); 211 212 CHKERRQ(VecAssemblyBegin(F)); 213 CHKERRQ(VecAssemblyEnd(F)); 214 PetscFunctionReturn(0); 215 } 216