1 #include <petsc/private/linesearchimpl.h>
2 #include <petsc/private/snesimpl.h>
3 #include <petscsnes.h>
4
SNESLineSearchApply_Secant(SNESLineSearch linesearch)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 i, max_it;
22 SNESObjectiveFn *objective;
23
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));
32
33 PetscCall(SNESGetObjective(snes, &objective, NULL));
34
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);
45
46 for (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 }
60
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));
76
77 /* compute the objective at the new endpoint */
78 PetscCall(VecWAXPY(W, -lambda, Y, X));
79 PetscCall(SNESComputeObjective(snes, W, &fnrm));
80 }
81
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 }
89
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 }
100
101 /* forbid the search from ever going back to the "failed" length that generates infinity or NaN */
102 maxlambda = .95 * lambda;
103
104 /* shrink lambda towards the previous one which was viable */
105 lambda = .5 * (lambda + lambda_old);
106 lambda_mid = .5 * (lambda + lambda_old);
107 }
108
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 }
119
120 /* determine change of lambda */
121 delLambda = lambda - lambda_old;
122
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 }
132
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;
136
137 /* compute f''() at the midpoint using centered differencing */
138 del2Fnrm = (delFnrm - delFnrm_old) / delLambda;
139
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 }
149
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 }
161
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) {
165 lambda_update = maxlambda;
166 break;
167 }
168
169 /* if lambda is NaN or Inf, do not accept update but exit with previous lambda */
170 if (PetscIsInfOrNanReal(lambda_update)) {
171 if (monitor) {
172 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
173 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambda_update is NaN or Inf\n"));
174 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
175 }
176 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
177 break;
178 }
179
180 /* update the endpoints and the midpoint of the bracketed secant region */
181 lambda_old = lambda;
182 lambda = lambda_update;
183 fnrm_old = fnrm;
184 lambda_mid = 0.5 * (lambda + lambda_old);
185
186 if ((i == max_it - 1) && monitor) {
187 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
188 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: reached maximum number of iterations!\n"));
189 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
190 }
191 }
192
193 /* construct the solution */
194 PetscCall(VecWAXPY(W, -lambda, Y, X));
195 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
196
197 /* postcheck */
198 PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
199 PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
200 if (changed_y) {
201 if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
202 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
203 }
204 PetscCall(VecCopy(W, X));
205 PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
206
207 PetscCall(SNESLineSearchComputeNorms(linesearch));
208 PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &gnorm, NULL));
209 SNESLineSearchCheckFunctionDomainError(snes, linesearch, gnorm);
210 if (monitor) {
211 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
212 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorm = %g\n", (double)lambda, (double)gnorm));
213 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
214 }
215 PetscFunctionReturn(PETSC_SUCCESS);
216 }
217
218 /*MC
219 SNESLINESEARCHSECANT - Secant search in the L2 norm of the function or the objective function.
220
221 Attempts to solve $ \min_{\lambda} f(x_k + \lambda Y_k) $ using the secant method with the initial bracketing of $ \lambda $ between [0,damping].
222 $f(x_k + \lambda Y_k)$ is either the squared L2-norm of the function $||F(x_k + \lambda Y_k)||_2^2$,
223 or the objective function $G(x_k + \lambda Y_k)$ if it is provided with `SNESSetObjective()`
224 Differences of $f()$ are used to approximate the first and second derivative of $f()$ with respect to
225 $\lambda$, $f'()$ and $f''()$.
226
227 When an objective function is provided $f(w)$ is the objective function otherwise $f(w) = ||F(w)||^2$.
228 $x$ is the current step and $y$ is the search direction.
229
230 Options Database Keys:
231 + -snes_linesearch_max_it <1> - maximum number of iterations within the line search
232 . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search
233 . -snes_linesearch_minlambda <1e\-12> - minimum allowable `lambda`
234 . -snes_linesearch_maxlambda <1.0> - maximum `lambda` (scaling of solution update) allowed
235 . -snes_linesearch_atol <1e\-15> - absolute tolerance for the secant method $ f'() < atol $
236 - -snes_linesearch_ltol <1e\-8> - minimum absolute change in `lambda` allowed
237
238 Level: advanced
239
240 .seealso: [](ch_snes), `SNESLINESEARCHBT`, `SNESLINESEARCHCP`, `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
241 M*/
SNESLineSearchCreate_Secant(SNESLineSearch linesearch)242 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Secant(SNESLineSearch linesearch)
243 {
244 PetscFunctionBegin;
245 linesearch->ops->apply = SNESLineSearchApply_Secant;
246 linesearch->ops->destroy = NULL;
247 linesearch->ops->setfromoptions = NULL;
248 linesearch->ops->reset = NULL;
249 linesearch->ops->view = NULL;
250 linesearch->ops->setup = NULL;
251
252 linesearch->max_it = 1;
253 PetscFunctionReturn(PETSC_SUCCESS);
254 }
255