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