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