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 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetObjectiveFunction(), SNESGetObjectiveFunction() 18 M*/ 19 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "SNESSetObjective" 23 /*@C 24 SNESSetObjective - Sets the objective function minimized by the SNES methods. 25 26 Logically Collective on SNES 27 28 Input Parameters: 29 + snes - the SNES context 30 . SNESObjectiveFunction - objective evaluation routine 31 - ctx - [optional] user-defined context for private data for the 32 function evaluation routine (may be PETSC_NULL) 33 34 Level: beginner 35 36 Note: If not provided then this defaults to the two norm of the function evaluation (set with SNESSetFunction()) 37 38 .keywords: SNES, nonlinear, set, objective 39 40 .seealso: SNESGetObjective(), SNESComputeObjective(), SNESSetFunction(), SNESSetJacobian(), SNESObjectiveFunction 41 @*/ 42 PetscErrorCode SNESSetObjective(SNES snes,PetscErrorCode (*SNESObjectiveFunction)(SNES,Vec,PetscReal *,void*),void *ctx) 43 { 44 PetscErrorCode ierr; 45 DM dm; 46 47 PetscFunctionBegin; 48 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 49 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 50 ierr = DMSNESSetObjective(dm,SNESObjectiveFunction,ctx);CHKERRQ(ierr); 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNCT__ 55 #define __FUNCT__ "SNESGetObjective" 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 + func - the function (or PETSC_NULL) 66 - ctx - the function context (or PETSC_NULL) 67 68 Level: advanced 69 70 .keywords: SNES, nonlinear, get, objective 71 72 .seealso: SNESSetObjective(), SNESGetSolution() 73 @*/ 74 PetscErrorCode SNESGetObjective(SNES snes,PetscErrorCode (**SNESObjectiveFunction)(SNES,Vec,PetscReal *,void*),void **ctx) 75 { 76 PetscErrorCode ierr; 77 DM dm; 78 PetscFunctionBegin; 79 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 80 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 81 ierr = DMSNESGetObjective(dm,SNESObjectiveFunction,ctx);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 #undef __FUNCT__ 86 #define __FUNCT__ "SNESComputeObjective" 87 /*@C 88 SNESComputeObjective - Computes the objective. 89 90 Collective on SNES 91 92 Input Parameter: 93 + snes - the SNES context 94 - X - the state vector 95 96 Output Parameter: 97 . ob - the objective value 98 99 Level: advanced 100 101 .keywords: SNES, nonlinear, compute, objective 102 103 .seealso: SNESSetObjective(), SNESGetSolution() 104 @*/ 105 PetscErrorCode SNESComputeObjective(SNES snes,Vec X,PetscReal *ob) 106 { 107 PetscErrorCode ierr; 108 DM dm; 109 DMSNES sdm; 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 = (sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx);CHKERRQ(ierr); 118 } else { 119 SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective()."); 120 } 121 PetscFunctionReturn(0); 122 } 123 124 125 #undef __FUNCT__ 126 #define __FUNCT__ "SNESDefaultObjectiveComputeFunctionFD" 127 PetscErrorCode SNESDefaultObjectiveComputeFunctionFD(SNES snes,Vec X,Vec F,void *ctx) 128 { 129 /* Quadratically interpolates the change in the objective based upon a change in a particular direction */ 130 131 Vec Xh; 132 PetscErrorCode ierr; 133 PetscInt i,N,start,end; 134 PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6; 135 PetscScalar fv,xv; 136 137 PetscFunctionBegin; 138 ierr = VecDuplicate(X,&Xh);CHKERRQ(ierr); 139 ierr = PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,PETSC_NULL);CHKERRQ(ierr); 140 ierr = VecSet(F,0.);CHKERRQ(ierr); 141 142 ierr = VecNorm(X,NORM_2,&fob);CHKERRQ(ierr); 143 144 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 145 ierr = VecGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 146 ierr = SNESComputeObjective(snes,X,&ob);CHKERRQ(ierr); 147 148 if (fob > 0.) { 149 dx =1e-6*fob; 150 } else { 151 dx = 1e-6; 152 } 153 154 for (i=0; i<N; i++) { 155 /* compute the 1st value */ 156 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 157 if (i>= start && i<end) { 158 xv = dx; 159 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 160 } 161 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 162 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 163 ierr = SNESComputeObjective(snes,Xh,&ob1);CHKERRQ(ierr); 164 165 /* compute the 2nd value */ 166 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 167 if (i>= start && i<end) { 168 xv = 2.*dx; 169 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 170 } 171 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 172 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 173 ierr = SNESComputeObjective(snes,Xh,&ob2);CHKERRQ(ierr); 174 175 /* compute the 3rd value */ 176 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 177 if (i>= start && i<end) { 178 xv = -dx; 179 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 180 } 181 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 182 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 183 ierr = SNESComputeObjective(snes,Xh,&ob3);CHKERRQ(ierr); 184 185 if (i >= start && i<end) { 186 /* set this entry to be the gradient of the objective */ 187 fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx); 188 if (PetscAbsScalar(fv) > eps) { 189 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr); 190 } else { 191 fv = 0.; 192 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr); 193 } 194 } 195 } 196 197 ierr = VecDestroy(&Xh);CHKERRQ(ierr); 198 199 ierr = VecAssemblyBegin(F);CHKERRQ(ierr); 200 ierr = VecAssemblyEnd(F);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203