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