Actual source code: linesearchsecant.c
1: #include <petsc/private/linesearchimpl.h>
2: #include <petscsnes.h>
4: static PetscErrorCode SNESLineSearchApply_Secant(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 minlambda, maxlambda, 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_it;
21: SNESObjectiveFn *objective;
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, &minlambda, &maxlambda, NULL, &atol, <ol, &max_it));
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_it; i++) {
45: /* check whether new lambda is NaN or Inf - if so, iteratively shrink towards lambda_old */
46: while (PETSC_TRUE) {
47: PetscCall(VecWAXPY(W, -lambda_mid, Y, X));
48: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
49: if (!objective) {
50: /* compute the norm at the midpoint */
51: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
52: if (linesearch->ops->vinorm) {
53: fnrm_mid = gnorm;
54: PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid));
55: } else {
56: PetscCall(VecNorm(F, NORM_2, &fnrm_mid));
57: }
59: /* compute the norm at the new endpoint */
60: PetscCall(VecWAXPY(W, -lambda, Y, X));
61: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
62: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
63: if (linesearch->ops->vinorm) {
64: fnrm = gnorm;
65: PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm));
66: } else {
67: PetscCall(VecNorm(F, NORM_2, &fnrm));
68: }
69: fnrm_mid = fnrm_mid * fnrm_mid;
70: fnrm = fnrm * fnrm;
71: } else {
72: /* compute the objective at the midpoint */
73: PetscCall(SNESComputeObjective(snes, W, &fnrm_mid));
75: /* compute the objective at the new endpoint */
76: PetscCall(VecWAXPY(W, -lambda, Y, X));
77: PetscCall(SNESComputeObjective(snes, W, &fnrm));
78: }
80: /* if new endpoint is viable, exit */
81: if (!PetscIsInfOrNanReal(fnrm)) break;
82: if (monitor) {
83: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
84: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambda = %g is Inf or Nan, cutting lambda\n", (double)lambda));
85: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
86: }
88: /* if smallest allowable lambda gives NaN or Inf, exit line search */
89: if (lambda <= minlambda) {
90: if (monitor) {
91: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
92: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambda = %g <= lambda_min = %g is Inf or Nan, can not further cut lambda\n", (double)lambda, (double)lambda));
93: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
94: }
95: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: /* forbid the search from ever going back to the "failed" length that generates Nan or Inf */
100: maxlambda = .95 * lambda;
102: /* shrink lambda towards the previous one which was viable */
103: lambda = .5 * (lambda + lambda_old);
104: lambda_mid = .5 * (lambda + lambda_old);
105: }
107: /* monitoring output */
108: if (monitor) {
109: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
110: if (!objective) {
111: 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)));
112: } else {
113: 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));
114: }
115: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
116: }
118: /* determine change of lambda */
119: delLambda = lambda - lambda_old;
121: /* check change of lambda tolerance */
122: if (PetscAbsReal(delLambda) < ltol) {
123: if (monitor) {
124: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
125: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(delLambda) = %g < ltol = %g\n", (double)PetscAbsReal(delLambda), (double)ltol));
126: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
127: }
128: break;
129: }
131: /* compute f'() at the end points using second order one sided differencing */
132: delFnrm = (3. * fnrm - 4. * fnrm_mid + 1. * fnrm_old) / delLambda;
133: delFnrm_old = (-3. * fnrm_old + 4. * fnrm_mid - 1. * fnrm) / delLambda;
135: /* compute f''() at the midpoint using centered differencing */
136: del2Fnrm = (delFnrm - delFnrm_old) / delLambda;
138: /* check absolute tolerance */
139: if (PetscAbsReal(delFnrm) <= atol) {
140: if (monitor) {
141: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
142: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(delFnrm) = %g <= atol = %g\n", (double)PetscAbsReal(delFnrm), (double)atol));
143: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
144: }
145: break;
146: }
148: /* compute the secant (Newton) update -- always go downhill */
149: if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
150: else if (del2Fnrm < 0.) lambda_update = lambda + delFnrm / del2Fnrm;
151: else {
152: if (monitor) {
153: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
154: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: del2Fnrm = 0\n"));
155: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
156: }
157: break;
158: }
160: /* prevent secant method from stepping out of bounds */
161: if (lambda_update < minlambda) lambda_update = 0.5 * (lambda + lambda_old);
162: if (lambda_update > maxlambda) {
163: lambda_update = maxlambda;
164: break;
165: }
167: /* if lambda is NaN or Inf, do not accept update but exit with previous lambda */
168: if (PetscIsInfOrNanReal(lambda_update)) {
169: if (monitor) {
170: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
171: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambda_update is NaN or Inf\n"));
172: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
173: }
174: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
175: break;
176: }
178: /* update the endpoints and the midpoint of the bracketed secant region */
179: lambda_old = lambda;
180: lambda = lambda_update;
181: fnrm_old = fnrm;
182: lambda_mid = 0.5 * (lambda + lambda_old);
184: if ((i == max_it - 1) && monitor) {
185: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
186: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: reached maximum number of iterations!\n"));
187: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
188: }
189: }
191: /* construct the solution */
192: PetscCall(VecWAXPY(W, -lambda, Y, X));
193: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
195: /* postcheck */
196: PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
197: PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
198: if (changed_y) {
199: if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
200: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
201: }
202: PetscCall(VecCopy(W, X));
203: PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
205: PetscCall(SNESLineSearchComputeNorms(linesearch));
207: if (monitor) {
208: PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &gnorm, NULL));
209: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
210: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorm = %g\n", (double)lambda, (double)gnorm));
211: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
212: }
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: /*MC
217: SNESLINESEARCHSECANT - Secant search in the L2 norm of the function or the objective function.
219: Attempts to solve $ \min_{\lambda} f(x_k + \lambda Y_k) $ using the secant method with the initial bracketing of $ \lambda $ between [0,damping].
220: $f(x_k + \lambda Y_k)$ is either the squared L2-norm of the function $||F(x_k + \lambda Y_k)||_2^2$,
221: or the objective function $G(x_k + \lambda Y_k)$ if it is provided with `SNESSetObjective()`
222: Differences of $f()$ are used to approximate the first and second derivative of $f()$ with respect to
223: $\lambda$, $f'()$ and $f''()$.
225: When an objective function is provided $f(w)$ is the objective function otherwise $f(w) = ||F(w)||^2$.
226: $x$ is the current step and $y$ is the search direction.
228: Options Database Keys:
229: + -snes_linesearch_max_it <1> - maximum number of iterations within the line search
230: . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search
231: . -snes_linesearch_minlambda <1e\-12> - minimum allowable `lambda`
232: . -snes_linesearch_maxlambda <1.0> - maximum `lambda` (scaling of solution update) allowed
233: . -snes_linesearch_atol <1e\-15> - absolute tolerance for the secant method $ f'() < atol $
234: - -snes_linesearch_ltol <1e\-8> - minimum absolute change in `lambda` allowed
236: Level: advanced
238: .seealso: [](ch_snes), `SNESLINESEARCHBT`, `SNESLINESEARCHCP`, `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
239: M*/
240: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Secant(SNESLineSearch linesearch)
241: {
242: PetscFunctionBegin;
243: linesearch->ops->apply = SNESLineSearchApply_Secant;
244: linesearch->ops->destroy = NULL;
245: linesearch->ops->setfromoptions = NULL;
246: linesearch->ops->reset = NULL;
247: linesearch->ops->view = NULL;
248: linesearch->ops->setup = NULL;
250: linesearch->max_it = 1;
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }