Actual source code: linesearchl2.c
petsc-3.9.4 2018-09-11
1: #include <petsc/private/linesearchimpl.h>
2: #include <petscsnes.h>
4: static PetscErrorCode SNESLineSearchApply_L2(SNESLineSearch linesearch)
5: {
7: PetscBool changed_y, changed_w;
9: Vec X;
10: Vec F;
11: Vec Y;
12: Vec W;
13: SNES snes;
14: PetscReal gnorm;
15: PetscReal ynorm;
16: PetscReal xnorm;
17: PetscReal steptol, maxstep, rtol, atol, ltol;
19: PetscViewer monitor;
20: PetscReal lambda, lambda_old, lambda_mid, lambda_update, delLambda;
21: PetscReal fnrm, fnrm_old, fnrm_mid;
22: PetscReal delFnrm, delFnrm_old, del2Fnrm;
23: PetscInt i, max_its;
25: PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
28: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
29: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
30: SNESLineSearchGetLambda(linesearch, &lambda);
31: SNESLineSearchGetSNES(linesearch, &snes);
32: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
33: SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);
34: SNESLineSearchGetDefaultMonitor(linesearch, &monitor);
36: SNESGetObjective(snes,&objective,NULL);
38: /* precheck */
39: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
40: lambda_old = 0.0;
41: if (!objective) {
42: fnrm_old = gnorm*gnorm;
43: } else {
44: SNESComputeObjective(snes,X,&fnrm_old);
45: }
46: lambda_mid = 0.5*(lambda + lambda_old);
48: for (i = 0; i < max_its; i++) {
50: while (PETSC_TRUE) {
51: VecCopy(X, W);
52: VecAXPY(W, -lambda_mid, Y);
53: if (linesearch->ops->viproject) {
54: (*linesearch->ops->viproject)(snes, W);
55: }
56: if (!objective) {
57: /* compute the norm at the midpoint */
58: (*linesearch->ops->snesfunc)(snes, W, F);
59: if (linesearch->ops->vinorm) {
60: fnrm_mid = gnorm;
61: (*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid);
62: } else {
63: VecNorm(F,NORM_2,&fnrm_mid);
64: }
66: /* compute the norm at the new endpoit */
67: VecCopy(X, W);
68: VecAXPY(W, -lambda, Y);
69: if (linesearch->ops->viproject) {
70: (*linesearch->ops->viproject)(snes, W);
71: }
72: (*linesearch->ops->snesfunc)(snes, W, F);
73: if (linesearch->ops->vinorm) {
74: fnrm = gnorm;
75: (*linesearch->ops->vinorm)(snes, F, W, &fnrm);
76: } else {
77: VecNorm(F,NORM_2,&fnrm);
78: }
79: fnrm_mid = fnrm_mid*fnrm_mid;
80: fnrm = fnrm*fnrm;
81: } else {
82: /* compute the objective at the midpoint */
83: VecCopy(X, W);
84: VecAXPY(W, -lambda_mid, Y);
85: SNESComputeObjective(snes,W,&fnrm_mid);
87: /* compute the objective at the new endpoint */
88: VecCopy(X, W);
89: VecAXPY(W, -lambda, Y);
90: SNESComputeObjective(snes,W,&fnrm);
91: }
92: if (!PetscIsInfOrNanReal(fnrm)) break;
93: if (monitor) {
94: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
95: PetscViewerASCIIPrintf(monitor," Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda);
96: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
97: }
98: if (lambda <= steptol) {
99: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
100: return(0);
101: }
102: maxstep = .95*lambda; /* forbid the search from ever going back to the "failed" length that generates Nan or Inf */
103: lambda = .5*(lambda + lambda_old);
104: lambda_mid = .5*(lambda + lambda_old);
105: }
107: delLambda = lambda - lambda_old;
108: /* compute f'() at the end points using second order one sided differencing */
109: delFnrm = (3.*fnrm - 4.*fnrm_mid + 1.*fnrm_old) / delLambda;
110: delFnrm_old = (-3.*fnrm_old + 4.*fnrm_mid -1.*fnrm) / delLambda;
111: /* compute f''() at the midpoint using centered differencing */
112: del2Fnrm = (delFnrm - delFnrm_old) / delLambda;
114: if (monitor) {
115: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
116: if (!objective) {
117: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g, %g], fnorms = [%g, %g, %g]\n",(double)lambda, (double)lambda_mid, (double)lambda_old, (double)PetscSqrtReal(fnrm), (double)PetscSqrtReal(fnrm_mid), (double)PetscSqrtReal(fnrm_old));
118: } else {
119: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g, %g], obj = [%g, %g, %g]\n",(double)lambda, (double)lambda_mid, (double)lambda_old, (double)fnrm, (double)fnrm_mid, (double)fnrm_old);
120: }
121: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
122: }
124: /* compute the secant (Newton) update -- always go downhill */
125: if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
126: else if (del2Fnrm < 0.) lambda_update = lambda + delFnrm / del2Fnrm;
127: else break;
129: if (lambda_update < steptol) lambda_update = 0.5*(lambda + lambda_old);
131: if (PetscIsInfOrNanReal(lambda_update)) break;
133: if (lambda_update > maxstep) break;
135: /* update the endpoints and the midpoint of the bracketed secant region */
136: lambda_old = lambda;
137: lambda = lambda_update;
138: fnrm_old = fnrm;
139: lambda_mid = 0.5*(lambda + lambda_old);
140: }
141: /* construct the solution */
142: VecCopy(X, W);
143: VecAXPY(W, -lambda, Y);
144: if (linesearch->ops->viproject) {
145: (*linesearch->ops->viproject)(snes, W);
146: }
148: /* postcheck */
149: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
150: if (changed_y) {
151: VecAXPY(X, -lambda, Y);
152: if (linesearch->ops->viproject) {
153: (*linesearch->ops->viproject)(snes, X);
154: }
155: } else {
156: VecCopy(W, X);
157: }
158: (*linesearch->ops->snesfunc)(snes,X,F);
160: SNESLineSearchSetLambda(linesearch, lambda);
161: SNESLineSearchComputeNorms(linesearch);
162: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
164: if (monitor) {
165: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
166: PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
167: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
168: }
169: if (lambda <= steptol) {
170: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
171: }
172: return(0);
173: }
175: /*MC
176: SNESLINESEARCHL2 - Secant search in the L2 norm of the function or the objective function, if it is provided with SNESSetObjective().
178: Attempts to solve min_lambda f(x + lambda y) using the secant method with the initial bracketing of lambda between [0,damping]. Differences of f()
179: are used to approximate the first and second derivative of f() with respect to lambda, f'() and f''(). The secant method is run for maxit iterations.
181: When an objective function is provided f(w) is the objective function otherwise f(w) = ||F(w)||^2. x is the current step and y is the search direction.
183: This has no checks on whether the secant method is actually converging.
185: Options Database Keys:
186: + -snes_linesearch_max_it <maxit> - maximum number of iterations, default is 1
187: . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
188: . -snes_linesearch_damping <damping> - initial step is scaled back by this factor, default is 1.0
189: - -snes_linesearch_minlambda <minlambda> - minimum allowable lambda
191: Level: advanced
193: Developer Notes: A better name for this method might be SNESLINESEARCHSECANT, L2 is not descriptive
195: .keywords: SNES, nonlinear, line search, norm, secant
197: .seealso: SNESLINESEARCHBT, SNESLINESEARCHCP, SNESLineSearch, SNESLineSearchCreate(), SNESLineSearchSetType()
198: M*/
199: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
200: {
202: linesearch->ops->apply = SNESLineSearchApply_L2;
203: linesearch->ops->destroy = NULL;
204: linesearch->ops->setfromoptions = NULL;
205: linesearch->ops->reset = NULL;
206: linesearch->ops->view = NULL;
207: linesearch->ops->setup = NULL;
209: linesearch->max_its = 1;
210: return(0);
211: }