Actual source code: linesearchl2.c
1: #include <petsc/private/linesearchimpl.h>
2: #include <petscsnes.h>
4: static PetscErrorCode SNESLineSearchApply_L2(SNESLineSearch linesearch)
5: {
6: PetscBool changed_y, changed_w;
7: Vec X;
8: Vec F;
9: Vec Y;
10: Vec W;
11: SNES snes;
12: PetscReal gnorm;
13: PetscReal ynorm;
14: PetscReal xnorm;
15: PetscReal steptol, maxstep, rtol, atol, ltol;
16: PetscViewer monitor;
17: PetscReal lambda, lambda_old, lambda_mid, lambda_update, delLambda;
18: PetscReal fnrm, fnrm_old, fnrm_mid;
19: PetscReal delFnrm, delFnrm_old, del2Fnrm;
20: PetscInt i, max_its;
21: PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);
23: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
24: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
25: SNESLineSearchGetLambda(linesearch, &lambda);
26: SNESLineSearchGetSNES(linesearch, &snes);
27: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
28: SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);
29: SNESLineSearchGetDefaultMonitor(linesearch, &monitor);
31: SNESGetObjective(snes, &objective, NULL);
33: /* precheck */
34: SNESLineSearchPreCheck(linesearch, X, Y, &changed_y);
35: lambda_old = 0.0;
36: if (!objective) {
37: fnrm_old = gnorm * gnorm;
38: } else {
39: SNESComputeObjective(snes, X, &fnrm_old);
40: }
41: lambda_mid = 0.5 * (lambda + lambda_old);
43: for (i = 0; i < max_its; i++) {
44: while (PETSC_TRUE) {
45: VecWAXPY(W, -lambda_mid, Y, X);
46: if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
47: if (!objective) {
48: /* compute the norm at the midpoint */
49: (*linesearch->ops->snesfunc)(snes, W, F);
50: if (linesearch->ops->vinorm) {
51: fnrm_mid = gnorm;
52: (*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid);
53: } else {
54: VecNorm(F, NORM_2, &fnrm_mid);
55: }
57: /* compute the norm at the new endpoit */
58: VecWAXPY(W, -lambda, Y, X);
59: if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
60: (*linesearch->ops->snesfunc)(snes, W, F);
61: if (linesearch->ops->vinorm) {
62: fnrm = gnorm;
63: (*linesearch->ops->vinorm)(snes, F, W, &fnrm);
64: } else {
65: VecNorm(F, NORM_2, &fnrm);
66: }
67: fnrm_mid = fnrm_mid * fnrm_mid;
68: fnrm = fnrm * fnrm;
69: } else {
70: /* compute the objective at the midpoint */
71: VecWAXPY(W, -lambda_mid, Y, X);
72: SNESComputeObjective(snes, W, &fnrm_mid);
74: /* compute the objective at the new endpoint */
75: VecWAXPY(W, -lambda, Y, X);
76: SNESComputeObjective(snes, W, &fnrm);
77: }
78: if (!PetscIsInfOrNanReal(fnrm)) break;
79: if (monitor) {
80: PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel);
81: PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda);
82: PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel);
83: }
84: if (lambda <= steptol) {
85: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
86: return 0;
87: }
88: maxstep = .95 * lambda; /* forbid the search from ever going back to the "failed" length that generates Nan or Inf */
89: lambda = .5 * (lambda + lambda_old);
90: lambda_mid = .5 * (lambda + lambda_old);
91: }
93: delLambda = lambda - lambda_old;
94: /* compute f'() at the end points using second order one sided differencing */
95: delFnrm = (3. * fnrm - 4. * fnrm_mid + 1. * fnrm_old) / delLambda;
96: delFnrm_old = (-3. * fnrm_old + 4. * fnrm_mid - 1. * fnrm) / delLambda;
97: /* compute f''() at the midpoint using centered differencing */
98: del2Fnrm = (delFnrm - delFnrm_old) / delLambda;
100: if (monitor) {
101: PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel);
102: if (!objective) {
103: 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));
104: } else {
105: 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);
106: }
107: PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel);
108: }
110: /* compute the secant (Newton) update -- always go downhill */
111: if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
112: else if (del2Fnrm < 0.) lambda_update = lambda + delFnrm / del2Fnrm;
113: else break;
115: if (lambda_update < steptol) lambda_update = 0.5 * (lambda + lambda_old);
117: if (PetscIsInfOrNanReal(lambda_update)) break;
119: if (lambda_update > maxstep) break;
121: /* update the endpoints and the midpoint of the bracketed secant region */
122: lambda_old = lambda;
123: lambda = lambda_update;
124: fnrm_old = fnrm;
125: lambda_mid = 0.5 * (lambda + lambda_old);
126: }
127: /* construct the solution */
128: VecWAXPY(W, -lambda, Y, X);
129: if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
131: /* postcheck */
132: SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w);
133: if (changed_y) {
134: VecAXPY(X, -lambda, Y);
135: if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, X);
136: } else {
137: VecCopy(W, X);
138: }
139: (*linesearch->ops->snesfunc)(snes, X, F);
141: SNESLineSearchSetLambda(linesearch, lambda);
142: SNESLineSearchComputeNorms(linesearch);
143: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
145: if (monitor) {
146: PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel);
147: PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
148: PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel);
149: }
150: if (lambda <= steptol) SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
151: return 0;
152: }
154: /*MC
155: SNESLINESEARCHL2 - Secant search in the L2 norm of the function or the objective function, if it is provided with `SNESSetObjective()`.
157: 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()
158: 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.
160: 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.
162: This has no checks on whether the secant method is actually converging.
164: Options Database Keys:
165: + -snes_linesearch_max_it <maxit> - maximum number of iterations, default is 1
166: . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
167: . -snes_linesearch_damping <damping> - initial step is scaled back by this factor, default is 1.0
168: - -snes_linesearch_minlambda <minlambda> - minimum allowable lambda
170: Level: advanced
172: Developer Note:
173: A better name for this method might be `SNESLINESEARCHSECANT`, L2 is not descriptive
175: .seealso: `SNESLINESEARCHBT`, `SNESLINESEARCHCP`, `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
176: M*/
177: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
178: {
179: linesearch->ops->apply = SNESLineSearchApply_L2;
180: linesearch->ops->destroy = NULL;
181: linesearch->ops->setfromoptions = NULL;
182: linesearch->ops->reset = NULL;
183: linesearch->ops->view = NULL;
184: linesearch->ops->setup = NULL;
186: linesearch->max_its = 1;
187: return 0;
188: }