Actual source code: linesearchsecant.c
1: #include <petsc/private/linesearchimpl.h>
2: #include <petsc/private/snesimpl.h>
3: #include <petscsnes.h>
5: static PetscErrorCode SNESLineSearchApply_Secant(SNESLineSearch linesearch)
6: {
7: PetscBool changed_y, changed_w;
8: Vec X;
9: Vec F;
10: Vec Y;
11: Vec W;
12: SNES snes;
13: PetscReal gnorm;
14: PetscReal ynorm;
15: PetscReal xnorm;
16: PetscReal minlambda, maxlambda, atol, ltol;
17: PetscViewer monitor;
18: PetscReal lambda, lambda_old, lambda_mid, lambda_update, delLambda;
19: PetscReal fnrm, fnrm_old, fnrm_mid;
20: PetscReal delFnrm, delFnrm_old, del2Fnrm;
21: PetscInt max_it;
22: SNESObjectiveFn *objective;
24: PetscFunctionBegin;
25: PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL));
26: PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
27: PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
28: PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
29: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
30: PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxlambda, NULL, &atol, <ol, &max_it));
31: PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
33: PetscCall(SNESGetObjective(snes, &objective, NULL));
35: /* precheck */
36: PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
37: lambda_old = 0.0;
38: if (!objective) {
39: fnrm_old = gnorm * gnorm;
40: } else {
41: PetscCall(SNESComputeObjective(snes, X, &fnrm_old));
42: SNESLineSearchCheckObjectiveDomainError(snes, fnrm_old);
43: }
44: lambda_mid = 0.5 * (lambda + lambda_old);
46: for (PetscInt i = 0; i < max_it; i++) {
47: /* repeatedly cut lambda until the norm or objective function is not infinity or NaN or lambda is too small */
48: while (PETSC_TRUE) {
49: PetscCall(VecWAXPY(W, -lambda_mid, Y, X));
50: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
51: if (!objective) {
52: /* compute the norm at the midpoint */
53: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
54: if (linesearch->ops->vinorm) {
55: fnrm_mid = gnorm;
56: PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid));
57: } else {
58: PetscCall(VecNorm(F, NORM_2, &fnrm_mid));
59: }
61: /* compute the norm at the new endpoint */
62: PetscCall(VecWAXPY(W, -lambda, Y, X));
63: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
64: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
65: if (linesearch->ops->vinorm) {
66: fnrm = gnorm;
67: PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm));
68: } else {
69: PetscCall(VecNorm(F, NORM_2, &fnrm));
70: }
71: fnrm_mid = fnrm_mid * fnrm_mid;
72: fnrm = fnrm * fnrm;
73: } else {
74: /* compute the objective at the midpoint */
75: PetscCall(SNESComputeObjective(snes, W, &fnrm_mid));
77: /* compute the objective at the new endpoint */
78: PetscCall(VecWAXPY(W, -lambda, Y, X));
79: PetscCall(SNESComputeObjective(snes, W, &fnrm));
80: }
82: /* if new endpoint is viable, exit */
83: if (!PetscIsInfOrNanReal(fnrm)) break;
84: if (monitor) {
85: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
86: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambda = %g is infinity or NaN, cutting lambda\n", (double)lambda));
87: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
88: }
90: /* if smallest allowable lambda gives NaN or Inf, exit line search */
91: if (lambda <= minlambda) {
92: if (monitor) {
93: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
94: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambda = %g <= lambda_min = %g is infinity or NaN, can not further cut lambda\n", (double)lambda, (double)lambda));
95: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
96: }
97: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: /* forbid the search from ever going back to the "failed" length that generates infinity or NaN */
102: maxlambda = .95 * lambda;
104: /* shrink lambda towards the previous one which was viable */
105: lambda = .5 * (lambda + lambda_old);
106: lambda_mid = .5 * (lambda + lambda_old);
107: }
109: /* monitoring output */
110: if (monitor) {
111: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
112: if (!objective) {
113: 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)));
114: } else {
115: 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));
116: }
117: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
118: }
120: /* determine change of lambda */
121: delLambda = lambda - lambda_old;
123: /* check change of lambda tolerance */
124: if (PetscAbsReal(delLambda) < ltol) {
125: if (monitor) {
126: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
127: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(delLambda) = %g < ltol = %g\n", (double)PetscAbsReal(delLambda), (double)ltol));
128: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
129: }
130: break;
131: }
133: /* compute f'() at the end points using second order one sided differencing */
134: delFnrm = (3. * fnrm - 4. * fnrm_mid + 1. * fnrm_old) / delLambda;
135: delFnrm_old = (-3. * fnrm_old + 4. * fnrm_mid - 1. * fnrm) / delLambda;
137: /* compute f''() at the midpoint using centered differencing */
138: del2Fnrm = (delFnrm - delFnrm_old) / delLambda;
140: /* check absolute tolerance */
141: if (PetscAbsReal(delFnrm) <= atol) {
142: if (monitor) {
143: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
144: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(delFnrm) = %g <= atol = %g\n", (double)PetscAbsReal(delFnrm), (double)atol));
145: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
146: }
147: break;
148: }
150: /* compute the secant (Newton) update -- always go downhill */
151: if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
152: else if (del2Fnrm < 0.) lambda_update = lambda + delFnrm / del2Fnrm;
153: else {
154: if (monitor) {
155: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
156: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: del2Fnrm = 0\n"));
157: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
158: }
159: break;
160: }
162: /* prevent secant method from stepping out of bounds */
163: if (lambda_update < minlambda) lambda_update = 0.5 * (lambda + lambda_old);
164: if (lambda_update > maxlambda) break;
166: /* if lambda is NaN or Inf, do not accept update but exit with previous lambda */
167: if (PetscIsInfOrNanReal(lambda_update)) {
168: if (monitor) {
169: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
170: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambda_update is NaN or Inf\n"));
171: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
172: }
173: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
174: break;
175: }
177: /* update the endpoints and the midpoint of the bracketed secant region */
178: lambda_old = lambda;
179: lambda = lambda_update;
180: fnrm_old = fnrm;
181: lambda_mid = 0.5 * (lambda + lambda_old);
183: if ((i == max_it - 1) && monitor) {
184: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
185: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: reached maximum number of iterations!\n"));
186: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
187: }
188: }
190: /* construct the solution */
191: PetscCall(VecWAXPY(W, -lambda, Y, X));
192: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
194: /* postcheck */
195: PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
196: PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
197: if (changed_y) {
198: if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
199: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
200: }
201: PetscCall(VecCopy(W, X));
202: PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
204: PetscCall(SNESLineSearchComputeNorms(linesearch));
205: PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &gnorm, NULL));
206: SNESLineSearchCheckFunctionDomainError(snes, linesearch, gnorm);
207: if (monitor) {
208: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
209: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorm = %g\n", (double)lambda, (double)gnorm));
210: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
211: }
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /*MC
216: SNESLINESEARCHSECANT - Secant search in the L2 norm of the function or the objective function.
218: Attempts to solve $ \min_{\lambda} f(x_k + \lambda Y_k) $ using the secant method with the initial bracketing of $ \lambda $ between [0,damping].
219: $f(x_k + \lambda Y_k)$ is either the squared L2-norm of the function $||F(x_k + \lambda Y_k)||_2^2$,
220: or the objective function $G(x_k + \lambda Y_k)$ if it is provided with `SNESSetObjective()`
221: Differences of $f()$ are used to approximate the first and second derivative of $f()$ with respect to
222: $\lambda$, $f'()$ and $f''()$.
224: When an objective function is provided $f(w)$ is the objective function otherwise $f(w) = ||F(w)||^2$.
225: $x$ is the current step and $y$ is the search direction.
227: Options Database Keys:
228: + -snes_linesearch_max_it 1 - maximum number of iterations within the line search
229: . -snes_linesearch_damping 1.0 - initial `lambda` on entry to the line search
230: . -snes_linesearch_minlambda 1e\-12 - minimum allowable `lambda`
231: . -snes_linesearch_maxlambda 1.0 - maximum `lambda` (scaling of solution update) allowed
232: . -snes_linesearch_atol 1e\-15 - absolute tolerance for the secant method $ f'() < atol $
233: - -snes_linesearch_ltol 1e\-8 - minimum absolute change in `lambda` allowed
235: Level: advanced
237: .seealso: [](ch_snes), `SNESLINESEARCHBT`, `SNESLINESEARCHCP`, `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
238: M*/
239: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Secant(SNESLineSearch linesearch)
240: {
241: PetscFunctionBegin;
242: linesearch->ops->apply = SNESLineSearchApply_Secant;
243: linesearch->ops->destroy = NULL;
244: linesearch->ops->setfromoptions = NULL;
245: linesearch->ops->reset = NULL;
246: linesearch->ops->view = NULL;
247: linesearch->ops->setup = NULL;
249: linesearch->max_it = 1;
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }