Actual source code: snesob.c
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: . obj - real to hold the objective value
14: - ctx - optional user-defined objective context
16: Level: advanced
18: .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetObjective(), SNESGetObjective(), SNESJacobianFunction, SNESFunction
19: M*/
22: /*@C
23: SNESSetObjective - Sets the objective function minimized by some of the SNES linesearch methods.
25: Logically Collective on SNES
27: Input Parameters:
28: + snes - the SNES context
29: . obj - objective evaluation routine; see SNESObjectiveFunction for details
30: - ctx - [optional] user-defined context for private data for the
31: function evaluation routine (may be NULL)
33: Level: intermediate
35: Note: This is not used in the SNESLINESEARCHCP line search.
37: If not provided then this defaults to the two norm of the function evaluation (set with SNESSetFunction())
39: .seealso: SNESGetObjective(), SNESComputeObjective(), SNESSetFunction(), SNESSetJacobian(), SNESObjectiveFunction
40: @*/
41: PetscErrorCode SNESSetObjective(SNES snes,PetscErrorCode (*obj)(SNES,Vec,PetscReal*,void*),void *ctx)
42: {
44: DM dm;
48: SNESGetDM(snes,&dm);
49: DMSNESSetObjective(dm,obj,ctx);
50: return(0);
51: }
53: /*@C
54: SNESGetObjective - Returns the objective function.
56: Not Collective
58: Input Parameter:
59: . snes - the SNES context
61: Output Parameter:
62: + obj - objective evaluation routine (or NULL); see SNESObjectFunction for details
63: - ctx - the function context (or NULL)
65: Level: advanced
67: .seealso: SNESSetObjective(), SNESGetSolution()
68: @*/
69: PetscErrorCode SNESGetObjective(SNES snes,PetscErrorCode (**obj)(SNES,Vec,PetscReal*,void*),void **ctx)
70: {
72: DM dm;
76: SNESGetDM(snes,&dm);
77: DMSNESGetObjective(dm,obj,ctx);
78: return(0);
79: }
81: /*@C
82: SNESComputeObjective - Computes the objective.
84: Collective on SNES
86: Input Parameter:
87: + snes - the SNES context
88: - X - the state vector
90: Output Parameter:
91: . ob - the objective value
93: Level: advanced
95: .seealso: SNESSetObjective(), SNESGetSolution()
96: @*/
97: PetscErrorCode SNESComputeObjective(SNES snes,Vec X,PetscReal *ob)
98: {
100: DM dm;
101: DMSNES sdm;
107: SNESGetDM(snes,&dm);
108: DMGetDMSNES(dm,&sdm);
109: if (sdm->ops->computeobjective) {
110: PetscLogEventBegin(SNES_ObjectiveEval,snes,X,0,0);
111: (sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx);
112: PetscLogEventEnd(SNES_ObjectiveEval,snes,X,0,0);
113: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
114: return(0);
115: }
118: /*@C
119: SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective
121: Collective on SNES
123: Input Parameter:
124: + snes - the SNES context
125: . X - the state vector
126: - ctx - the (ignored) function context
128: Output Parameter:
129: . F - the function value
131: Options Database Key:
132: + -snes_fd_function_eps - The differencing parameter
133: - -snes_fd_function - Compute function from user provided objective with finite difference
135: Notes:
136: SNESObjectiveComputeFunctionDefaultFD is similar in character to SNESComputeJacobianDefault.
137: Therefore, it should be used for debugging purposes only. Using it in conjunction with
138: SNESComputeJacobianDefault is excessively costly and produces a Jacobian that is quite
139: noisy. This is often necessary, but should be done with a grain of salt, even when debugging
140: small problems.
142: Note that this uses quadratic interpolation of the objective to form each value in the function.
144: Level: advanced
146: .seealso: SNESSetFunction(), SNESComputeObjective(), SNESComputeJacobianDefault()
147: @*/
148: PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,void *ctx)
149: {
150: Vec Xh;
152: PetscInt i,N,start,end;
153: PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6;
154: PetscScalar fv,xv;
157: VecDuplicate(X,&Xh);
158: PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing parameters","SNES");
159: PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL);
160: PetscOptionsEnd();
161: VecSet(F,0.);
163: VecNorm(X,NORM_2,&fob);
165: VecGetSize(X,&N);
166: VecGetOwnershipRange(X,&start,&end);
167: SNESComputeObjective(snes,X,&ob);
169: if (fob > 0.) dx =1e-6*fob;
170: else dx = 1e-6;
172: for (i=0; i<N; i++) {
173: /* compute the 1st value */
174: VecCopy(X,Xh);
175: if (i>= start && i<end) {
176: xv = dx;
177: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
178: }
179: VecAssemblyBegin(Xh);
180: VecAssemblyEnd(Xh);
181: SNESComputeObjective(snes,Xh,&ob1);
183: /* compute the 2nd value */
184: VecCopy(X,Xh);
185: if (i>= start && i<end) {
186: xv = 2.*dx;
187: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
188: }
189: VecAssemblyBegin(Xh);
190: VecAssemblyEnd(Xh);
191: SNESComputeObjective(snes,Xh,&ob2);
193: /* compute the 3rd value */
194: VecCopy(X,Xh);
195: if (i>= start && i<end) {
196: xv = -dx;
197: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
198: }
199: VecAssemblyBegin(Xh);
200: VecAssemblyEnd(Xh);
201: SNESComputeObjective(snes,Xh,&ob3);
203: if (i >= start && i<end) {
204: /* set this entry to be the gradient of the objective */
205: fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx);
206: if (PetscAbsScalar(fv) > eps) {
207: VecSetValues(F,1,&i,&fv,INSERT_VALUES);
208: } else {
209: fv = 0.;
210: VecSetValues(F,1,&i,&fv,INSERT_VALUES);
211: }
212: }
213: }
215: VecDestroy(&Xh);
217: VecAssemblyBegin(F);
218: VecAssemblyEnd(F);
219: return(0);
220: }