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