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(), SNESSetObjectiveFunction(), SNESGetObjectiveFunction() 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 + func - the function (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 PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,void *ctx) 130 { 131 /* Quadratically interpolates the change in the objective based upon a change in a particular direction */ 132 133 Vec Xh; 134 PetscErrorCode ierr; 135 PetscInt i,N,start,end; 136 PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6; 137 PetscScalar fv,xv; 138 139 PetscFunctionBegin; 140 ierr = VecDuplicate(X,&Xh);CHKERRQ(ierr); 141 ierr = PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL);CHKERRQ(ierr); 142 ierr = VecSet(F,0.);CHKERRQ(ierr); 143 144 ierr = VecNorm(X,NORM_2,&fob);CHKERRQ(ierr); 145 146 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 147 ierr = VecGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 148 ierr = SNESComputeObjective(snes,X,&ob);CHKERRQ(ierr); 149 150 if (fob > 0.) dx =1e-6*fob; 151 else dx = 1e-6; 152 153 for (i=0; i<N; i++) { 154 /* compute the 1st value */ 155 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 156 if (i>= start && i<end) { 157 xv = dx; 158 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 159 } 160 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 161 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 162 ierr = SNESComputeObjective(snes,Xh,&ob1);CHKERRQ(ierr); 163 164 /* compute the 2nd value */ 165 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 166 if (i>= start && i<end) { 167 xv = 2.*dx; 168 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 169 } 170 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 171 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 172 ierr = SNESComputeObjective(snes,Xh,&ob2);CHKERRQ(ierr); 173 174 /* compute the 3rd value */ 175 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 176 if (i>= start && i<end) { 177 xv = -dx; 178 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 179 } 180 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 181 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 182 ierr = SNESComputeObjective(snes,Xh,&ob3);CHKERRQ(ierr); 183 184 if (i >= start && i<end) { 185 /* set this entry to be the gradient of the objective */ 186 fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx); 187 if (PetscAbsScalar(fv) > eps) { 188 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr); 189 } else { 190 fv = 0.; 191 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr); 192 } 193 } 194 } 195 196 ierr = VecDestroy(&Xh);CHKERRQ(ierr); 197 198 ierr = VecAssemblyBegin(F);CHKERRQ(ierr); 199 ierr = VecAssemblyEnd(F);CHKERRQ(ierr); 200 PetscFunctionReturn(0); 201 } 202