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 /* Quadratically interpolates the change in the objective based upon a change in a particular direction */ 129 130 Vec Xh; 131 PetscErrorCode ierr; 132 PetscInt i,N,start,end; 133 PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6; 134 PetscScalar fv,xv; 135 136 PetscFunctionBegin; 137 ierr = VecDuplicate(X,&Xh);CHKERRQ(ierr); 138 ierr = PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,PETSC_NULL);CHKERRQ(ierr); 139 ierr = VecSet(F,0.);CHKERRQ(ierr); 140 141 ierr = VecNorm(X,NORM_2,&fob);CHKERRQ(ierr); 142 143 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 144 ierr = VecGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 145 ierr = SNESComputeObjective(snes,X,&ob);CHKERRQ(ierr); 146 147 if (fob > 0.) { 148 dx =1e-6*fob; 149 } else { 150 dx = 1e-6; 151 } 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