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