Actual source code: snesob.c

petsc-3.5.3 2015-01-31
Report Typos and Errors
  1: #include <petsc-private/snesimpl.h>

  3: /*MC
  4:     SNESObjectiveFunction - functional form used to convey the objective function to the nonlinear solver

  6:      Synopsis:
  7:      #include <petscsnes.h>
  8:        SNESObjectiveFunction(SNES snes,Vec x,PetscReal *obj,void *ctx);

 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

 17:    Level: advanced

 19: .seealso:   SNESSetFunction(), SNESGetFunction(), SNESSetObjective(), SNESGetObjective(), SNESJacobianFunction, SNESFunction
 20: M*/


 25: /*@C
 26:    SNESSetObjective - Sets the objective function minimized by some of the SNES linesearch methods.

 28:    Logically Collective on SNES

 30:    Input Parameters:
 31: +  snes - the SNES context
 32: .  obj - objective evaluation routine; see SNESObjectiveFunction for details
 33: -  ctx - [optional] user-defined context for private data for the
 34:          function evaluation routine (may be NULL)

 36:    Level: intermediately

 38:    Note: This is not used in the SNESLINESEARCHCP line search.

 40:          If not provided then this defaults to the two norm of the function evaluation (set with SNESSetFunction())

 42: .keywords: SNES, nonlinear, set, objective

 44: .seealso: SNESGetObjective(), SNESComputeObjective(), SNESSetFunction(), SNESSetJacobian(), SNESObjectiveFunction
 45: @*/
 46: PetscErrorCode  SNESSetObjective(SNES snes,PetscErrorCode (*obj)(SNES,Vec,PetscReal*,void*),void *ctx)
 47: {
 49:   DM             dm;

 53:   SNESGetDM(snes,&dm);
 54:   DMSNESSetObjective(dm,obj,ctx);
 55:   return(0);
 56: }

 60: /*@C
 61:    SNESGetObjective - Returns the objective function.

 63:    Not Collective

 65:    Input Parameter:
 66: .  snes - the SNES context

 68:    Output Parameter:
 69: +  obj - objective evaluation routine (or NULL); see SNESObjectFunction for details
 70: -  ctx - the function context (or NULL)

 72:    Level: advanced

 74: .keywords: SNES, nonlinear, get, objective

 76: .seealso: SNESSetObjective(), SNESGetSolution()
 77: @*/
 78: PetscErrorCode SNESGetObjective(SNES snes,PetscErrorCode (**obj)(SNES,Vec,PetscReal*,void*),void **ctx)
 79: {
 81:   DM             dm;

 85:   SNESGetDM(snes,&dm);
 86:   DMSNESGetObjective(dm,obj,ctx);
 87:   return(0);
 88: }

 92: /*@C
 93:    SNESComputeObjective - Computes the objective.

 95:    Collective on SNES

 97:    Input Parameter:
 98: +  snes - the SNES context
 99: -  X    - the state vector

101:    Output Parameter:
102: .  ob   - the objective value

104:    Level: advanced

106: .keywords: SNES, nonlinear, compute, objective

108: .seealso: SNESSetObjective(), SNESGetSolution()
109: @*/
110: PetscErrorCode SNESComputeObjective(SNES snes,Vec X,PetscReal *ob)
111: {
113:   DM             dm;
114:   DMSNES         sdm;

120:   SNESGetDM(snes,&dm);
121:   DMGetDMSNES(dm,&sdm);
122:   if (sdm->ops->computeobjective) {
123:     (sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx);
124:   } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
125:   return(0);
126: }


131: /*@C
132:    SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective

134:    Collective on SNES

136:    Input Parameter:
137: +  snes - the SNES context
138: .  X    - the state vector
139: -  ctx  - the (ignored) function context

141:    Output Parameter:
142: .  F   - the function value

144:    Options Database Key:
145: +  -snes_fd_function_eps - The differencing parameter
146: -  -snes_fd_function - Compute function from user provided objective with finite difference

148:    Notes:
149:    SNESObjectiveComputeFunctionDefaultFD is similar in character to SNESComputeJacobianDefault.
150:    Therefore, it should be used for debugging purposes only.  Using it in conjunction with
151:    SNESComputeJacobianDefault is excessively costly and produces a Jacobian that is quite
152:    noisy.  This is often necessary, but should be done with a grain of salt, even when debugging
153:    small problems.

155:    Note that this uses quadratic interpolation of the objective to form each value in the function.

157:    Level: advanced

159: .keywords: SNES, objective, debugging, finite differences, function

161: .seealso: SNESSetFunction(), SNESComputeObjective(), SNESComputeJacobianDefault()
162: @*/
163: PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,void *ctx)
164: {
165:   Vec            Xh;
167:   PetscInt       i,N,start,end;
168:   PetscReal      ob,ob1,ob2,ob3,fob,dx,eps=1e-6;
169:   PetscScalar    fv,xv;

172:   VecDuplicate(X,&Xh);
173:   PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL);
174:   VecSet(F,0.);

176:   VecNorm(X,NORM_2,&fob);

178:   VecGetSize(X,&N);
179:   VecGetOwnershipRange(X,&start,&end);
180:   SNESComputeObjective(snes,X,&ob);

182:   if (fob > 0.) dx =1e-6*fob;
183:   else dx = 1e-6;

185:   for (i=0; i<N; i++) {
186:     /* compute the 1st value */
187:     VecCopy(X,Xh);
188:     if (i>= start && i<end) {
189:       xv   = dx;
190:       VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
191:     }
192:     VecAssemblyBegin(Xh);
193:     VecAssemblyEnd(Xh);
194:     SNESComputeObjective(snes,Xh,&ob1);

196:     /* compute the 2nd value */
197:     VecCopy(X,Xh);
198:     if (i>= start && i<end) {
199:       xv   = 2.*dx;
200:       VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
201:     }
202:     VecAssemblyBegin(Xh);
203:     VecAssemblyEnd(Xh);
204:     SNESComputeObjective(snes,Xh,&ob2);

206:     /* compute the 3rd value */
207:     VecCopy(X,Xh);
208:     if (i>= start && i<end) {
209:       xv   = -dx;
210:       VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
211:     }
212:     VecAssemblyBegin(Xh);
213:     VecAssemblyEnd(Xh);
214:     SNESComputeObjective(snes,Xh,&ob3);

216:     if (i >= start && i<end) {
217:       /* set this entry to be the gradient of the objective */
218:       fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx);
219:       if (PetscAbsScalar(fv) > eps) {
220:         VecSetValues(F,1,&i,&fv,INSERT_VALUES);
221:       } else {
222:         fv   = 0.;
223:         VecSetValues(F,1,&i,&fv,INSERT_VALUES);
224:       }
225:     }
226:   }

228:   VecDestroy(&Xh);

230:   VecAssemblyBegin(F);
231:   VecAssemblyEnd(F);
232:   return(0);
233: }