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 = (sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx);CHKERRQ(ierr); 124 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective()."); 125 PetscFunctionReturn(0); 126 } 127 128 129 #undef __FUNCT__ 130 #define __FUNCT__ "SNESObjectiveComputeFunctionDefaultFD" 131 /*@C 132 SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective 133 134 Collective on SNES 135 136 Input Parameter: 137 + snes - the SNES context 138 . X - the state vector 139 - ctx - the (ignored) function context 140 141 Output Parameter: 142 . F - the function value 143 144 Options Database Key: 145 + -snes_fd_function_eps - The differencing parameter 146 - -snes_fd_function - Compute function from user provided objective with finite difference 147 148 Notes: 149 SNESObjectiveComputeFunctionDefaultFD is similar in character to SNESComputeJacobianDefault. 150 Therefore, it should be used for debugging purposes only. Using it in conjunction with 151 SNESComputeJacobianDefault is excessively costly and produces a Jacobian that is quite 152 noisy. This is often necessary, but should be done with a grain of salt, even when debugging 153 small problems. 154 155 Note that this uses quadratic interpolation of the objective to form each value in the function. 156 157 Level: advanced 158 159 .keywords: SNES, objective, debugging, finite differences, function 160 161 .seealso: SNESSetFunction(), SNESComputeObjective(), SNESComputeJacobianDefault() 162 @*/ 163 PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,void *ctx) 164 { 165 Vec Xh; 166 PetscErrorCode ierr; 167 PetscInt i,N,start,end; 168 PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6; 169 PetscScalar fv,xv; 170 171 PetscFunctionBegin; 172 ierr = VecDuplicate(X,&Xh);CHKERRQ(ierr); 173 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing parameters","SNES");CHKERRQ(ierr); 174 ierr = PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL);CHKERRQ(ierr); 175 ierr = PetscOptionsEnd();CHKERRQ(ierr); 176 ierr = VecSet(F,0.);CHKERRQ(ierr); 177 178 ierr = VecNorm(X,NORM_2,&fob);CHKERRQ(ierr); 179 180 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 181 ierr = VecGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 182 ierr = SNESComputeObjective(snes,X,&ob);CHKERRQ(ierr); 183 184 if (fob > 0.) dx =1e-6*fob; 185 else dx = 1e-6; 186 187 for (i=0; i<N; i++) { 188 /* compute the 1st value */ 189 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 190 if (i>= start && i<end) { 191 xv = dx; 192 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 193 } 194 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 195 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 196 ierr = SNESComputeObjective(snes,Xh,&ob1);CHKERRQ(ierr); 197 198 /* compute the 2nd value */ 199 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 200 if (i>= start && i<end) { 201 xv = 2.*dx; 202 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 203 } 204 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 205 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 206 ierr = SNESComputeObjective(snes,Xh,&ob2);CHKERRQ(ierr); 207 208 /* compute the 3rd value */ 209 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 210 if (i>= start && i<end) { 211 xv = -dx; 212 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 213 } 214 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 215 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 216 ierr = SNESComputeObjective(snes,Xh,&ob3);CHKERRQ(ierr); 217 218 if (i >= start && i<end) { 219 /* set this entry to be the gradient of the objective */ 220 fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx); 221 if (PetscAbsScalar(fv) > eps) { 222 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr); 223 } else { 224 fv = 0.; 225 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr); 226 } 227 } 228 } 229 230 ierr = VecDestroy(&Xh);CHKERRQ(ierr); 231 232 ierr = VecAssemblyBegin(F);CHKERRQ(ierr); 233 ierr = VecAssemblyEnd(F);CHKERRQ(ierr); 234 PetscFunctionReturn(0); 235 } 236