1 #include <petsc/private/linesearchimpl.h>
2 #include <petsc/private/snesimpl.h>
3
SNESLineSearchApply_Bisection(SNESLineSearch linesearch)4 static PetscErrorCode SNESLineSearchApply_Bisection(SNESLineSearch linesearch)
5 {
6 PetscBool changed_y, changed_w;
7 Vec X, F, Y, W, G;
8 SNES snes;
9 PetscReal ynorm;
10 PetscReal lambda_left, lambda, lambda_right, lambda_old, fnorm;
11 PetscScalar fty_left, fty, fty_initial;
12 PetscViewer monitor;
13 PetscReal rtol, atol, ltol;
14 PetscInt it, max_it;
15
16 PetscFunctionBegin;
17 PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G));
18 PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
19 PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
20 PetscCall(SNESLineSearchGetTolerances(linesearch, NULL, NULL, &rtol, &atol, <ol, &max_it));
21 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
22 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
23
24 /* pre-check */
25 PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
26
27 /* compute ynorm to normalize search direction */
28 PetscCall(VecNorm(Y, NORM_2, &ynorm));
29
30 /* initialize interval for bisection */
31 lambda_left = 0.0;
32 lambda_right = lambda;
33
34 /* compute fty at left end of interval */
35 if (linesearch->ops->vidirderiv) {
36 PetscCall((*linesearch->ops->vidirderiv)(snes, F, X, Y, &fty_left));
37 } else {
38 PetscCall(VecDot(F, Y, &fty_left));
39 }
40 fty_initial = fty_left;
41
42 /* compute fty at right end of interval (initial lambda) */
43 PetscCall(VecWAXPY(W, -lambda, Y, X));
44 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
45 PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
46 if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
47 PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations during line search!\n"));
48 snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
49 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
50 PetscFunctionReturn(PETSC_SUCCESS);
51 }
52 if (linesearch->ops->vidirderiv) {
53 PetscCall((*linesearch->ops->vidirderiv)(snes, G, W, Y, &fty));
54 } else {
55 PetscCall(VecDot(G, Y, &fty));
56 }
57 /* check whether sign changes in interval */
58 if (!PetscIsInfOrNanScalar(fty) && (PetscRealPart(fty_left * fty) > 0.0)) {
59 /* no change of sign: accept full step */
60 if (monitor) {
61 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
62 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: sign of fty does not change in step interval, accepting full step\n"));
63 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
64 }
65 } else {
66 /* change of sign: iteratively bisect interval */
67 lambda_old = 0.0;
68 it = 0;
69
70 while (PETSC_TRUE) {
71 /* check for infinity or NaN */
72 if (PetscIsInfOrNanScalar(fty)) {
73 if (monitor) {
74 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
75 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search fty is infinity or NaN!\n"));
76 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
77 }
78 PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_CONV_FAILED, "infinity or NaN in function evaluation");
79 if (snes->functiondomainerror) {
80 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
81 snes->functiondomainerror = PETSC_FALSE;
82 } else PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
83 PetscFunctionReturn(PETSC_SUCCESS);
84 break;
85 }
86
87 /* check absolute tolerance */
88 if (PetscAbsScalar(fty) <= atol * ynorm) {
89 if (monitor) {
90 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
91 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty)/||y|| = %g <= atol = %g\n", (double)(PetscAbsScalar(fty) / ynorm), (double)atol));
92 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
93 }
94 break;
95 }
96
97 /* check relative tolerance */
98 if (PetscAbsScalar(fty) / PetscAbsScalar(fty_initial) <= rtol) {
99 if (monitor) {
100 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
101 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty/fty_initial) = %g <= rtol = %g\n", (double)(PetscAbsScalar(fty) / PetscAbsScalar(fty_initial)), (double)rtol));
102 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
103 }
104 break;
105 }
106
107 /* check maximum number of iterations */
108 if (it > max_it) {
109 if (monitor) {
110 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
111 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: maximum iterations reached\n"));
112 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
113 }
114 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
115 break;
116 }
117
118 /* check change of lambda tolerance */
119 if (PetscAbsReal(lambda - lambda_old) < ltol) {
120 if (monitor) {
121 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
122 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(dlambda) = %g < ltol = %g\n", (double)PetscAbsReal(lambda - lambda_old), (double)ltol));
123 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
124 }
125 break;
126 }
127
128 /* determine direction of bisection (not necessary for 0th iteration) */
129 if (it > 0) {
130 if (PetscRealPart(fty * fty_left) <= 0.0) {
131 lambda_right = lambda;
132 } else {
133 lambda_left = lambda;
134 /* also update fty_left for direction check in next iteration */
135 fty_left = fty;
136 }
137 }
138
139 /* bisect interval */
140 lambda_old = lambda;
141 lambda = 0.5 * (lambda_left + lambda_right);
142
143 /* compute fty at new lambda */
144 PetscCall(VecWAXPY(W, -lambda, Y, X));
145 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
146 PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
147 if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
148 PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations during line search!\n"));
149 snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
150 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
151 PetscFunctionReturn(PETSC_SUCCESS);
152 }
153 if (linesearch->ops->vidirderiv) {
154 PetscCall((*linesearch->ops->vidirderiv)(snes, G, W, Y, &fty));
155 } else {
156 PetscCall(VecDot(G, Y, &fty));
157 }
158
159 /* print iteration information */
160 if (monitor) {
161 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
162 PetscCall(PetscViewerASCIIPrintf(monitor, " %3" PetscInt_FMT " Line search: fty/||y|| = %g, lambda = %g\n", it, (double)(PetscRealPart(fty) / ynorm), (double)lambda));
163 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
164 }
165
166 /* count up iteration */
167 it++;
168 }
169 }
170
171 /* post-check */
172 PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
173 PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
174 if (changed_y) {
175 if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
176 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
177 }
178
179 /* update solution*/
180 PetscCall(VecCopy(W, X));
181 PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
182 PetscCall(SNESLineSearchComputeNorms(linesearch));
183 PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &fnorm, NULL));
184 SNESLineSearchCheckFunctionDomainError(snes, linesearch, fnorm);
185 PetscFunctionReturn(PETSC_SUCCESS);
186 }
187
188 /*MC
189 SNESLINESEARCHBISECTION - Bisection line search.
190 Similar to the critical point line search, `SNESLINESEARCHCP`, the bisection line search assumes that there exists some $G(x)$ for which the `SNESFunctionFn` $F(x) = grad G(x)$.
191 This line search seeks to find the root of the directional derivative, that is $F(x_k - \lambda Y_k) \cdot Y_k / ||Y_k|| = 0$, along the search direction $Y_k$ through bisection.
192
193 Options Database Keys:
194 + -snes_linesearch_max_it <50> - maximum number of bisection iterations for the line search
195 . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search
196 . -snes_linesearch_rtol <1e\-8> - relative tolerance for the directional derivative
197 . -snes_linesearch_atol <1e\-6> - absolute tolerance for the directional derivative
198 - -snes_linesearch_ltol <1e\-6> - minimum absolute change in `lambda` allowed (this is an alternative to setting a maximum number of iterations)
199
200 Level: intermediate
201
202 Notes:
203 `lambda` is the scaling of the search direction (vector) that is computed by this algorithm.
204 If there is no change of sign in the directional derivative from $\lambda=0$ to the initial `lambda` (the damping), then the initial `lambda` will be used.
205 Hence, this line search will always give a `lambda` in the interval $[0, damping]$.
206 This method does NOT use the objective function if it is provided with `SNESSetObjective()`.
207
208 .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHCP`
209 M*/
SNESLineSearchCreate_Bisection(SNESLineSearch linesearch)210 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Bisection(SNESLineSearch linesearch)
211 {
212 PetscFunctionBegin;
213 linesearch->ops->apply = SNESLineSearchApply_Bisection;
214 linesearch->ops->destroy = NULL;
215 linesearch->ops->setfromoptions = NULL;
216 linesearch->ops->reset = NULL;
217 linesearch->ops->view = NULL;
218 linesearch->ops->setup = NULL;
219
220 /* set default option values */
221 linesearch->max_it = 50;
222 linesearch->rtol = 1e-8;
223 linesearch->atol = 1e-6;
224 linesearch->ltol = 1e-6;
225 PetscFunctionReturn(PETSC_SUCCESS);
226 }
227