Actual source code: snesob.c
petsc-3.7.5 2017-01-01
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: intermediate
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: PetscLogEventBegin(SNES_ObjectiveEval,snes,X,0,0);
124: (sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx);
125: PetscLogEventEnd(SNES_ObjectiveEval,snes,X,0,0);
126: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
127: return(0);
128: }
133: /*@C
134: SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective
136: Collective on SNES
138: Input Parameter:
139: + snes - the SNES context
140: . X - the state vector
141: - ctx - the (ignored) function context
143: Output Parameter:
144: . F - the function value
146: Options Database Key:
147: + -snes_fd_function_eps - The differencing parameter
148: - -snes_fd_function - Compute function from user provided objective with finite difference
150: Notes:
151: SNESObjectiveComputeFunctionDefaultFD is similar in character to SNESComputeJacobianDefault.
152: Therefore, it should be used for debugging purposes only. Using it in conjunction with
153: SNESComputeJacobianDefault is excessively costly and produces a Jacobian that is quite
154: noisy. This is often necessary, but should be done with a grain of salt, even when debugging
155: small problems.
157: Note that this uses quadratic interpolation of the objective to form each value in the function.
159: Level: advanced
161: .keywords: SNES, objective, debugging, finite differences, function
163: .seealso: SNESSetFunction(), SNESComputeObjective(), SNESComputeJacobianDefault()
164: @*/
165: PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,void *ctx)
166: {
167: Vec Xh;
169: PetscInt i,N,start,end;
170: PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6;
171: PetscScalar fv,xv;
174: VecDuplicate(X,&Xh);
175: PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing parameters","SNES");
176: PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL);
177: PetscOptionsEnd();
178: VecSet(F,0.);
180: VecNorm(X,NORM_2,&fob);
182: VecGetSize(X,&N);
183: VecGetOwnershipRange(X,&start,&end);
184: SNESComputeObjective(snes,X,&ob);
186: if (fob > 0.) dx =1e-6*fob;
187: else dx = 1e-6;
189: for (i=0; i<N; i++) {
190: /* compute the 1st value */
191: VecCopy(X,Xh);
192: if (i>= start && i<end) {
193: xv = dx;
194: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
195: }
196: VecAssemblyBegin(Xh);
197: VecAssemblyEnd(Xh);
198: SNESComputeObjective(snes,Xh,&ob1);
200: /* compute the 2nd value */
201: VecCopy(X,Xh);
202: if (i>= start && i<end) {
203: xv = 2.*dx;
204: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
205: }
206: VecAssemblyBegin(Xh);
207: VecAssemblyEnd(Xh);
208: SNESComputeObjective(snes,Xh,&ob2);
210: /* compute the 3rd value */
211: VecCopy(X,Xh);
212: if (i>= start && i<end) {
213: xv = -dx;
214: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
215: }
216: VecAssemblyBegin(Xh);
217: VecAssemblyEnd(Xh);
218: SNESComputeObjective(snes,Xh,&ob3);
220: if (i >= start && i<end) {
221: /* set this entry to be the gradient of the objective */
222: fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx);
223: if (PetscAbsScalar(fv) > eps) {
224: VecSetValues(F,1,&i,&fv,INSERT_VALUES);
225: } else {
226: fv = 0.;
227: VecSetValues(F,1,&i,&fv,INSERT_VALUES);
228: }
229: }
230: }
232: VecDestroy(&Xh);
234: VecAssemblyBegin(F);
235: VecAssemblyEnd(F);
236: return(0);
237: }