Actual source code: linesearchl2.c
petsc-3.5.4 2015-05-23
1: #include <petsc-private/linesearchimpl.h>
2: #include <petscsnes.h>
6: static PetscErrorCode SNESLineSearchApply_L2(SNESLineSearch linesearch)
7: {
9: PetscBool changed_y, changed_w;
11: Vec X;
12: Vec F;
13: Vec Y;
14: Vec W;
15: SNES snes;
16: PetscReal gnorm;
17: PetscReal ynorm;
18: PetscReal xnorm;
19: PetscReal steptol, maxstep, rtol, atol, ltol;
21: PetscViewer monitor;
22: PetscBool domainerror;
23: PetscReal lambda, lambda_old, lambda_mid, lambda_update, delLambda;
24: PetscReal fnrm, fnrm_old, fnrm_mid;
25: PetscReal delFnrm, delFnrm_old, del2Fnrm;
26: PetscInt i, max_its;
28: PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
31: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
32: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
33: SNESLineSearchGetLambda(linesearch, &lambda);
34: SNESLineSearchGetSNES(linesearch, &snes);
35: SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);
36: SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);
37: SNESLineSearchGetMonitor(linesearch, &monitor);
39: SNESGetObjective(snes,&objective,NULL);
41: /* precheck */
42: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
43: lambda_old = 0.0;
44: if (!objective) {
45: fnrm_old = gnorm*gnorm;
46: } else {
47: SNESComputeObjective(snes,X,&fnrm_old);
48: }
49: lambda_mid = 0.5*(lambda + lambda_old);
51: for (i = 0; i < max_its; i++) {
53: /* compute the norm at the midpoint */
54: VecCopy(X, W);
55: VecAXPY(W, -lambda_mid, Y);
56: if (linesearch->ops->viproject) {
57: (*linesearch->ops->viproject)(snes, W);
58: }
59: if (!objective) {
60: (*linesearch->ops->snesfunc)(snes, W, F);
61: if (linesearch->ops->vinorm) {
62: fnrm_mid = gnorm;
63: (*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid);
64: } else {
65: VecNorm(F,NORM_2,&fnrm_mid);
66: }
68: /* compute the norm at lambda */
69: VecCopy(X, W);
70: VecAXPY(W, -lambda, Y);
71: if (linesearch->ops->viproject) {
72: (*linesearch->ops->viproject)(snes, W);
73: }
74: (*linesearch->ops->snesfunc)(snes, W, F);
75: if (linesearch->ops->vinorm) {
76: fnrm = gnorm;
77: (*linesearch->ops->vinorm)(snes, F, W, &fnrm);
78: } else {
79: VecNorm(F,NORM_2,&fnrm);
80: }
81: fnrm_mid = fnrm_mid*fnrm_mid;
82: fnrm = fnrm*fnrm;
83: } else {
84: /* compute the objective at the midpoint */
85: VecCopy(X, W);
86: VecAXPY(W, -lambda_mid, Y);
87: SNESComputeObjective(snes,W,&fnrm_mid);
89: /* compute the objective at the midpoint */
90: VecCopy(X, W);
91: VecAXPY(W, -lambda, Y);
92: SNESComputeObjective(snes,W,&fnrm);
93: }
94: /* this gives us the derivatives at the endpoints -- compute them from here
96: a = x - a0
98: p_0(x) = (x / dA - 1)(2x / dA - 1)
99: p_1(x) = 4(x / dA)(1 - x / dA)
100: p_2(x) = (x / dA)(2x / dA - 1)
102: dp_0[0] / dx = 3 / dA
103: dp_1[0] / dx = -4 / dA
104: dp_2[0] / dx = 1 / dA
106: dp_0[dA] / dx = -1 / dA
107: dp_1[dA] / dx = 4 / dA
108: dp_2[dA] / dx = -3 / dA
110: d^2p_0[0] / dx^2 = 4 / dA^2
111: d^2p_1[0] / dx^2 = -8 / dA^2
112: d^2p_2[0] / dx^2 = 4 / dA^2
113: */
115: delLambda = lambda - lambda_old;
116: delFnrm = (3.*fnrm - 4.*fnrm_mid + 1.*fnrm_old) / delLambda;
117: delFnrm_old = (-3.*fnrm_old + 4.*fnrm_mid -1.*fnrm) / delLambda;
118: del2Fnrm = (delFnrm - delFnrm_old) / delLambda;
120: if (monitor) {
121: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
122: if (!objective) {
123: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g, %g], fnorms = [%g, %g, %g]\n",
124: (double)lambda, (double)lambda_mid, (double)lambda_old, (double)PetscSqrtReal(fnrm), (double)PetscSqrtReal(fnrm_mid), (double)PetscSqrtReal(fnrm_old));
125: } else {
126: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g, %g], obj = [%g, %g, %g]\n",
127: (double)lambda, (double)lambda_mid, (double)lambda_old, (double)fnrm, (double)fnrm_mid, (double)fnrm_old);
129: }
130: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
131: }
133: /* compute the search direction -- always go downhill */
134: if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
135: else lambda_update = lambda + delFnrm / del2Fnrm;
137: if (lambda_update < steptol) lambda_update = 0.5*(lambda + lambda_old);
139: if (PetscIsInfOrNanScalar(lambda_update)) break;
141: if (lambda_update > maxstep) break;
143: /* compute the new state of the line search */
144: lambda_old = lambda;
145: lambda = lambda_update;
146: fnrm_old = fnrm;
147: lambda_mid = 0.5*(lambda + lambda_old);
148: }
149: /* construct the solution */
150: VecCopy(X, W);
151: VecAXPY(W, -lambda, Y);
152: if (linesearch->ops->viproject) {
153: (*linesearch->ops->viproject)(snes, W);
154: }
156: /* postcheck */
157: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
158: if (changed_y) {
159: VecAXPY(X, -lambda, Y);
160: if (linesearch->ops->viproject) {
161: (*linesearch->ops->viproject)(snes, X);
162: }
163: } else {
164: VecCopy(W, X);
165: }
166: (*linesearch->ops->snesfunc)(snes,X,F);
167: SNESGetFunctionDomainError(snes, &domainerror);
168: if (domainerror) {
169: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
170: return(0);
171: }
173: SNESLineSearchSetLambda(linesearch, lambda);
174: SNESLineSearchComputeNorms(linesearch);
175: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
177: if (monitor) {
178: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
179: PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
180: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
181: }
182: if (lambda <= steptol) {
183: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
184: }
185: return(0);
186: }
190: /*MC
191: SNESLINESEARCHL2 - Secant search in the L2 norm of the function or the objective function if it is provided with SNESSetObjective().
193: The function norm is evaluated at points in [0, damping] to construct
194: a polynomial fitting. This fitting is used to construct a new lambda
195: based upon secant descent. The process is repeated on the new
196: interval, [lambda, lambda_old], max_it - 1 times.
198: Options Database Keys:
199: + -snes_linesearch_max_it <maxit> - maximum number of iterations, default is 1
200: . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
201: . -snes_linesearch_damping <damping> - initial step is scaled back by this factor, default is 1.0
202: - -snes_linesearch_minlambda <minlambda> - minimum allowable lambda
204: Level: advanced
206: .keywords: SNES, nonlinear, line search, norm, secant
208: .seealso: SNESLineSearchBT, SNESLineSearchCP, SNESLineSearch
209: M*/
210: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
211: {
213: linesearch->ops->apply = SNESLineSearchApply_L2;
214: linesearch->ops->destroy = NULL;
215: linesearch->ops->setfromoptions = NULL;
216: linesearch->ops->reset = NULL;
217: linesearch->ops->view = NULL;
218: linesearch->ops->setup = NULL;
220: linesearch->max_its = 1;
221: return(0);
222: }