Actual source code: snesob.c
petsc-3.9.4 2018-09-11
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*/
23: /*@C
24: SNESSetObjective - Sets the objective function minimized by some of the SNES linesearch methods.
26: Logically Collective on SNES
28: Input Parameters:
29: + snes - the SNES context
30: . obj - objective evaluation routine; see SNESObjectiveFunction for details
31: - ctx - [optional] user-defined context for private data for the
32: function evaluation routine (may be NULL)
34: Level: intermediate
36: Note: This is not used in the SNESLINESEARCHCP line search.
38: 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 (*obj)(SNES,Vec,PetscReal*,void*),void *ctx)
45: {
47: DM dm;
51: SNESGetDM(snes,&dm);
52: DMSNESSetObjective(dm,obj,ctx);
53: return(0);
54: }
56: /*@C
57: SNESGetObjective - Returns the objective function.
59: Not Collective
61: Input Parameter:
62: . snes - the SNES context
64: Output Parameter:
65: + obj - objective evaluation routine (or NULL); see SNESObjectFunction for details
66: - ctx - the function context (or NULL)
68: Level: advanced
70: .keywords: SNES, nonlinear, get, objective
72: .seealso: SNESSetObjective(), SNESGetSolution()
73: @*/
74: PetscErrorCode SNESGetObjective(SNES snes,PetscErrorCode (**obj)(SNES,Vec,PetscReal*,void*),void **ctx)
75: {
77: DM dm;
81: SNESGetDM(snes,&dm);
82: DMSNESGetObjective(dm,obj,ctx);
83: return(0);
84: }
86: /*@C
87: SNESComputeObjective - Computes the objective.
89: Collective on SNES
91: Input Parameter:
92: + snes - the SNES context
93: - X - the state vector
95: Output Parameter:
96: . ob - the objective value
98: Level: advanced
100: .keywords: SNES, nonlinear, compute, objective
102: .seealso: SNESSetObjective(), SNESGetSolution()
103: @*/
104: PetscErrorCode SNESComputeObjective(SNES snes,Vec X,PetscReal *ob)
105: {
107: DM dm;
108: DMSNES sdm;
114: SNESGetDM(snes,&dm);
115: DMGetDMSNES(dm,&sdm);
116: if (sdm->ops->computeobjective) {
117: PetscLogEventBegin(SNES_ObjectiveEval,snes,X,0,0);
118: (sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx);
119: PetscLogEventEnd(SNES_ObjectiveEval,snes,X,0,0);
120: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
121: return(0);
122: }
125: /*@C
126: SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective
128: Collective on SNES
130: Input Parameter:
131: + snes - the SNES context
132: . X - the state vector
133: - ctx - the (ignored) function context
135: Output Parameter:
136: . F - the function value
138: Options Database Key:
139: + -snes_fd_function_eps - The differencing parameter
140: - -snes_fd_function - Compute function from user provided objective with finite difference
142: Notes:
143: SNESObjectiveComputeFunctionDefaultFD is similar in character to SNESComputeJacobianDefault.
144: Therefore, it should be used for debugging purposes only. Using it in conjunction with
145: SNESComputeJacobianDefault is excessively costly and produces a Jacobian that is quite
146: noisy. This is often necessary, but should be done with a grain of salt, even when debugging
147: small problems.
149: Note that this uses quadratic interpolation of the objective to form each value in the function.
151: Level: advanced
153: .keywords: SNES, objective, debugging, finite differences, function
155: .seealso: SNESSetFunction(), SNESComputeObjective(), SNESComputeJacobianDefault()
156: @*/
157: PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,void *ctx)
158: {
159: Vec Xh;
161: PetscInt i,N,start,end;
162: PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6;
163: PetscScalar fv,xv;
166: VecDuplicate(X,&Xh);
167: PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing parameters","SNES");
168: PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL);
169: PetscOptionsEnd();
170: VecSet(F,0.);
172: VecNorm(X,NORM_2,&fob);
174: VecGetSize(X,&N);
175: VecGetOwnershipRange(X,&start,&end);
176: SNESComputeObjective(snes,X,&ob);
178: if (fob > 0.) dx =1e-6*fob;
179: else dx = 1e-6;
181: for (i=0; i<N; i++) {
182: /* compute the 1st value */
183: VecCopy(X,Xh);
184: if (i>= start && i<end) {
185: xv = dx;
186: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
187: }
188: VecAssemblyBegin(Xh);
189: VecAssemblyEnd(Xh);
190: SNESComputeObjective(snes,Xh,&ob1);
192: /* compute the 2nd value */
193: VecCopy(X,Xh);
194: if (i>= start && i<end) {
195: xv = 2.*dx;
196: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
197: }
198: VecAssemblyBegin(Xh);
199: VecAssemblyEnd(Xh);
200: SNESComputeObjective(snes,Xh,&ob2);
202: /* compute the 3rd value */
203: VecCopy(X,Xh);
204: if (i>= start && i<end) {
205: xv = -dx;
206: VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
207: }
208: VecAssemblyBegin(Xh);
209: VecAssemblyEnd(Xh);
210: SNESComputeObjective(snes,Xh,&ob3);
212: if (i >= start && i<end) {
213: /* set this entry to be the gradient of the objective */
214: fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx);
215: if (PetscAbsScalar(fv) > eps) {
216: VecSetValues(F,1,&i,&fv,INSERT_VALUES);
217: } else {
218: fv = 0.;
219: VecSetValues(F,1,&i,&fv,INSERT_VALUES);
220: }
221: }
222: }
224: VecDestroy(&Xh);
226: VecAssemblyBegin(F);
227: VecAssemblyEnd(F);
228: return(0);
229: }