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: PetscFunctionBegin;
24: PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL));
25: PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
26: PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
27: PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
28: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
29: PetscCall(SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its));
30: PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
32: PetscCall(SNESGetObjective(snes, &objective, NULL));
34: /* precheck */
35: PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
36: lambda_old = 0.0;
37: if (!objective) {
38: fnrm_old = gnorm * gnorm;
39: } else {
40: PetscCall(SNESComputeObjective(snes, X, &fnrm_old));
41: }
42: lambda_mid = 0.5 * (lambda + lambda_old);
44: for (i = 0; i < max_its; i++) {
45: while (PETSC_TRUE) {
46: PetscCall(VecWAXPY(W, -lambda_mid, Y, X));
47: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
48: if (!objective) {
49: /* compute the norm at the midpoint */
50: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
51: if (linesearch->ops->vinorm) {
52: fnrm_mid = gnorm;
53: PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid));
54: } else {
55: PetscCall(VecNorm(F, NORM_2, &fnrm_mid));
56: }
58: /* compute the norm at the new endpoint */
59: PetscCall(VecWAXPY(W, -lambda, Y, X));
60: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
61: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
62: if (linesearch->ops->vinorm) {
63: fnrm = gnorm;
64: PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm));
65: } else {
66: PetscCall(VecNorm(F, NORM_2, &fnrm));
67: }
68: fnrm_mid = fnrm_mid * fnrm_mid;
69: fnrm = fnrm * fnrm;
70: } else {
71: /* compute the objective at the midpoint */
72: PetscCall(SNESComputeObjective(snes, W, &fnrm_mid));
74: /* compute the objective at the new endpoint */
75: PetscCall(VecWAXPY(W, -lambda, Y, X));
76: PetscCall(SNESComputeObjective(snes, W, &fnrm));
77: }
78: if (!PetscIsInfOrNanReal(fnrm)) break;
79: if (monitor) {
80: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
81: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda));
82: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
83: }
84: if (lambda <= steptol) {
85: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
86: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
102: if (!objective) {
103: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(VecWAXPY(W, -lambda, Y, X));
129: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
131: /* postcheck */
132: PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
133: PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
134: if (changed_y) {
135: if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
136: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
137: }
138: PetscCall(VecCopy(W, X));
139: PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
141: PetscCall(SNESLineSearchComputeNorms(linesearch));
143: if (monitor) {
144: PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &gnorm, NULL));
145: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
146: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorm = %g\n", (double)lambda, (double)gnorm));
147: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
148: }
149: if (lambda <= steptol) PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: /*MC
154: SNESLINESEARCHL2 - Secant search in the L2 norm of the function or the objective function, if it is provided with `SNESSetObjective()`.
156: Attempts to solve $ \min_{\lambda} f(x + \lambda y) $ using the secant method with the initial bracketing of $ \lambda $ between [0,damping].
157: Differences of $f()$ are used to approximate the first and second derivative of $f()$ with respect to
158: $\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$.
161: $x$ is the current step and $y$ is the search direction.
163: This has no checks on whether the secant method is actually converging.
165: Options Database Keys:
166: + -snes_linesearch_max_it <maxit> - maximum number of iterations, default is 1
167: . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
168: . -snes_linesearch_damping <damping> - initial step is scaled back by this factor, default is 1.0
169: - -snes_linesearch_minlambda <minlambda> - minimum allowable lambda
171: Level: advanced
173: Developer Note:
174: A better name for this method might be `SNESLINESEARCHSECANT`, L2 is not descriptive
176: .seealso: [](ch_snes), `SNESLINESEARCHBT`, `SNESLINESEARCHCP`, `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
177: M*/
178: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
179: {
180: PetscFunctionBegin;
181: linesearch->ops->apply = SNESLineSearchApply_L2;
182: linesearch->ops->destroy = NULL;
183: linesearch->ops->setfromoptions = NULL;
184: linesearch->ops->reset = NULL;
185: linesearch->ops->view = NULL;
186: linesearch->ops->setup = NULL;
188: linesearch->max_its = 1;
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }