Actual source code: snesob.c
petsc-3.4.5 2014-06-29
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()
20: M*/
25: /*@C
26: SNESSetObjective - Sets the objective function minimized by the SNES methods.
28: Logically Collective on SNES
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)
36: Level: beginner
38: Note: If not provided then this defaults to the two norm of the function evaluation (set with SNESSetFunction())
40: .keywords: SNES, nonlinear, set, objective
42: .seealso: SNESGetObjective(), SNESComputeObjective(), SNESSetFunction(), SNESSetJacobian(), SNESObjectiveFunction
43: @*/
44: PetscErrorCode SNESSetObjective(SNES snes,PetscErrorCode (*SNESObjectiveFunction)(SNES,Vec,PetscReal*,void*),void *ctx)
45: {
47: DM dm;
51: SNESGetDM(snes,&dm);
52: DMSNESSetObjective(dm,SNESObjectiveFunction,ctx);
53: return(0);
54: }
58: /*@C
59: SNESGetObjective - Returns the objective function.
61: Not Collective
63: Input Parameter:
64: . snes - the SNES context
66: Output Parameter:
67: + SNESObjectiveFunction - objective evaluation routine (or NULL)
68: - ctx - the function context (or NULL)
70: Level: advanced
72: .keywords: SNES, nonlinear, get, objective
74: .seealso: SNESSetObjective(), SNESGetSolution()
75: @*/
76: PetscErrorCode SNESGetObjective(SNES snes,PetscErrorCode (**SNESObjectiveFunction)(SNES,Vec,PetscReal*,void*),void **ctx)
77: {
79: DM dm;
83: SNESGetDM(snes,&dm);
84: DMSNESGetObjective(dm,SNESObjectiveFunction,ctx);
85: return(0);
86: }
90: /*@C
91: SNESComputeObjective - Computes the objective.
93: Collective on SNES
95: Input Parameter:
96: + snes - the SNES context
97: - X - the state vector
99: Output Parameter:
100: . ob - the objective value
102: Level: advanced
104: .keywords: SNES, nonlinear, compute, objective
106: .seealso: SNESSetObjective(), SNESGetSolution()
107: @*/
108: PetscErrorCode SNESComputeObjective(SNES snes,Vec X,PetscReal *ob)
109: {
111: DM dm;
112: DMSNES sdm;
118: SNESGetDM(snes,&dm);
119: DMGetDMSNES(dm,&sdm);
120: if (sdm->ops->computeobjective) {
121: (sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx);
122: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
123: return(0);
124: }
129: /*@C
130: SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective
132: Collective on SNES
134: Input Parameter:
135: + snes - the SNES context
136: . X - the state vector
137: - ctx - the (ignored) function context
139: Output Parameter:
140: . F - the function value
142: Options Database Key:
143: + -snes_fd_function_eps - The differencing parameter
144: - -snes_fd_function - Compute function from user provided objective with finite difference
146: Notes:
147: SNESObjectiveComputeFunctionDefaultFD is similar in character to SNESComputeJacobianDefault.
148: Therefore, it should be used for debugging purposes only. Using it in conjunction with
149: SNESComputeJacobianDefault is excessively costly and produces a Jacobian that is quite
150: noisy. This is often necessary, but should be done with a grain of salt, even when debugging
151: small problems.
153: Note that this uses quadratic interpolation of the objective to form each value in the function.
155: Level: advanced
157: .keywords: SNES, objective, debugging, finite differences, function
159: .seealso: SNESSetFunction(), SNESComputeObjective(), SNESComputeJacobianDefault()
160: @*/
161: PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,void *ctx)
162: {
163: Vec Xh;
165: PetscInt i,N,start,end;
166: PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6;
167: PetscScalar fv,xv;
170: VecDuplicate(X,&Xh);
171: PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL);
172: VecSet(F,0.);
174: VecNorm(X,NORM_2,&fob);
176: VecGetSize(X,&N);
177: VecGetOwnershipRange(X,&start,&end);
178: SNESComputeObjective(snes,X,&ob);
180: if (fob > 0.) dx =1e-6*fob;
181: else dx = 1e-6;
183: for (i=0; i<N; i++) {
184: /* compute the 1st value */
185: VecCopy(X,Xh);
186: if (i>= start && i<end) {
187: xv = dx;
188: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
189: }
190: VecAssemblyBegin(Xh);
191: VecAssemblyEnd(Xh);
192: SNESComputeObjective(snes,Xh,&ob1);
194: /* compute the 2nd value */
195: VecCopy(X,Xh);
196: if (i>= start && i<end) {
197: xv = 2.*dx;
198: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
199: }
200: VecAssemblyBegin(Xh);
201: VecAssemblyEnd(Xh);
202: SNESComputeObjective(snes,Xh,&ob2);
204: /* compute the 3rd value */
205: VecCopy(X,Xh);
206: if (i>= start && i<end) {
207: xv = -dx;
208: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
209: }
210: VecAssemblyBegin(Xh);
211: VecAssemblyEnd(Xh);
212: SNESComputeObjective(snes,Xh,&ob3);
214: if (i >= start && i<end) {
215: /* set this entry to be the gradient of the objective */
216: fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx);
217: if (PetscAbsScalar(fv) > eps) {
218: VecSetValues(F,1,&i,&fv,INSERT_VALUES);
219: } else {
220: fv = 0.;
221: VecSetValues(F,1,&i,&fv,INSERT_VALUES);
222: }
223: }
224: }
226: VecDestroy(&Xh);
228: VecAssemblyBegin(F);
229: VecAssemblyEnd(F);
230: return(0);
231: }