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