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: intermediately 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 = PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL);CHKERRQ(ierr); 174 ierr = VecSet(F,0.);CHKERRQ(ierr); 175 176 ierr = VecNorm(X,NORM_2,&fob);CHKERRQ(ierr); 177 178 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 179 ierr = VecGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 180 ierr = SNESComputeObjective(snes,X,&ob);CHKERRQ(ierr); 181 182 if (fob > 0.) dx =1e-6*fob; 183 else dx = 1e-6; 184 185 for (i=0; i<N; i++) { 186 /* compute the 1st value */ 187 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 188 if (i>= start && i<end) { 189 xv = dx; 190 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 191 } 192 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 193 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 194 ierr = SNESComputeObjective(snes,Xh,&ob1);CHKERRQ(ierr); 195 196 /* compute the 2nd value */ 197 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 198 if (i>= start && i<end) { 199 xv = 2.*dx; 200 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 201 } 202 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 203 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 204 ierr = SNESComputeObjective(snes,Xh,&ob2);CHKERRQ(ierr); 205 206 /* compute the 3rd value */ 207 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 208 if (i>= start && i<end) { 209 xv = -dx; 210 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 211 } 212 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 213 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 214 ierr = SNESComputeObjective(snes,Xh,&ob3);CHKERRQ(ierr); 215 216 if (i >= start && i<end) { 217 /* set this entry to be the gradient of the objective */ 218 fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx); 219 if (PetscAbsScalar(fv) > eps) { 220 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr); 221 } else { 222 fv = 0.; 223 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr); 224 } 225 } 226 } 227 228 ierr = VecDestroy(&Xh);CHKERRQ(ierr); 229 230 ierr = VecAssemblyBegin(F);CHKERRQ(ierr); 231 ierr = VecAssemblyEnd(F);CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234