1 #include <petsc/private/linesearchimpl.h>
2 #include <petscsnes.h>
3
SNESLineSearchApply_CP(SNESLineSearch linesearch)4 static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch)
5 {
6 PetscBool changed_y, changed_w;
7 Vec X, Y, F, W;
8 SNES snes;
9 PetscReal xnorm, ynorm, gnorm, minlambda, maxlambda, rtol, atol, ltol;
10 PetscReal lambda, lambda_old, lambda_update, delLambda;
11 PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s;
12 PetscInt i, max_it;
13 PetscViewer monitor;
14
15 PetscFunctionBegin;
16 PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL));
17 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
18 PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
19 PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
20 PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxlambda, &rtol, &atol, <ol, &max_it));
21 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
22 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
23
24 /* precheck */
25 PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
26 lambda_old = 0.0;
27
28 /* evaluate initial point */
29 if (linesearch->ops->vidirderiv) {
30 PetscCall((*linesearch->ops->vidirderiv)(snes, F, X, Y, &fty_old));
31 } else {
32 PetscCall(VecDot(F, Y, &fty_old));
33 }
34 if (PetscAbsScalar(fty_old) < atol * ynorm) {
35 if (monitor) {
36 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
37 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated at initial point because dot(F,Y) = %g < atol*||y|| = %g\n", (double)PetscAbsScalar(fty_old), (double)(atol * ynorm)));
38 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
39 }
40 PetscCall(SNESSetConvergedReason(linesearch->snes, SNES_CONVERGED_FNORM_ABS));
41 PetscFunctionReturn(PETSC_SUCCESS);
42 }
43 fty_init = fty_old;
44
45 for (i = 0; i < max_it; i++) {
46 /* compute the norm at lambda */
47 PetscCall(VecWAXPY(W, -lambda, Y, X));
48 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
49 PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
50 if (linesearch->ops->vidirderiv) {
51 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty));
52 } else {
53 PetscCall(VecDot(F, Y, &fty));
54 }
55
56 /* compute change of lambda */
57 delLambda = lambda - lambda_old;
58
59 /* check change of lambda tolerance */
60 if (PetscAbsReal(delLambda) < ltol) {
61 if (monitor) {
62 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
63 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(dlambda) = %g < ltol = %g\n", (double)PetscAbsReal(delLambda), (double)ltol));
64 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
65 }
66 break;
67 }
68
69 /* check relative tolerance */
70 if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) {
71 if (monitor) {
72 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
73 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty/fty_init) = %g <= rtol = %g\n", (double)(PetscAbsScalar(fty) / PetscAbsScalar(fty_init)), (double)rtol));
74 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
75 }
76 break;
77 }
78
79 /* check absolute tolerance */
80 if (PetscAbsScalar(fty) < atol * ynorm && i > 0) {
81 if (monitor) {
82 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
83 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty)/||y|| = %g <= atol = %g\n", (double)(PetscAbsScalar(fty) / ynorm), (double)atol));
84 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
85 }
86 break;
87 }
88
89 /* print iteration information */
90 if (monitor) {
91 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
92 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambdas = [%g, %g], ftys = [%g, %g]\n", (double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old)));
93 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
94 }
95
96 /* approximate the second derivative */
97 if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) {
98 /* first order finite difference approximation */
99 s = (fty - fty_old) / delLambda;
100 } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
101 /* evaluate function at midpoint 0.5 * (lambda + lambda_old) */
102 PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X));
103 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
104 PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
105 if (linesearch->ops->vidirderiv) {
106 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1));
107 } else {
108 PetscCall(VecDot(F, Y, &fty_mid1));
109 }
110
111 /* second order finite difference approximation */
112 s = (3. * fty - 4. * fty_mid1 + fty_old) / delLambda;
113
114 } else {
115 /* evaluate function at midpoint 0.5 * (lambda + lambda_old) */
116 PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X));
117 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
118 PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
119 if (linesearch->ops->vidirderiv) {
120 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1));
121 } else {
122 PetscCall(VecDot(F, Y, &fty_mid1));
123 }
124
125 /* evaluate function at lambda + 0.5 * (lambda - lambda_old) */
126 PetscCall(VecWAXPY(W, -(lambda + 0.5 * (lambda - lambda_old)), Y, X));
127 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
128 PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
129 if (linesearch->ops->vidirderiv) {
130 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid2));
131 } else {
132 PetscCall(VecDot(F, Y, &fty_mid2));
133 }
134
135 /* third order finite difference approximation */
136 s = (2. * fty_mid2 + 3. * fty - 6. * fty_mid1 + fty_old) / (3. * delLambda);
137 }
138
139 /* compute secant update (modifying the search direction if necessary) */
140 if (PetscRealPart(s) > 0.) s = -s;
141 if (s == 0.0) {
142 if (monitor) {
143 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
144 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: del2Fnrm = 0\n"));
145 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
146 }
147 break;
148 }
149 lambda_update = lambda - PetscRealPart(fty / s);
150
151 /* if step is too small, go the opposite direction */
152 if (lambda_update < minlambda) lambda_update = lambda + PetscRealPart(fty / s);
153 /* if secant method would step out of bounds, exit with the respective bound */
154 if (lambda_update > maxlambda) {
155 lambda = maxlambda;
156 break;
157 }
158
159 /* if lambda is infinity or NaN, do not accept update but exit with previous lambda */
160 if (PetscIsInfOrNanReal(lambda_update)) {
161 if (monitor) {
162 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
163 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambda_update is infinity or NaN\n"));
164 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
165 }
166 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
167 break;
168 }
169
170 /* compute the new state of the line search */
171 lambda_old = lambda;
172 lambda = lambda_update;
173 fty_old = fty;
174 }
175
176 /* construct the solution */
177 PetscCall(VecWAXPY(W, -lambda, Y, X));
178 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
179 /* postcheck */
180 PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
181 PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
182 if (changed_y) {
183 if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
184 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
185 }
186 PetscCall(VecCopy(W, X));
187 PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
188
189 PetscCall(SNESLineSearchComputeNorms(linesearch));
190 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
191
192 if (monitor) {
193 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
194 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm));
195 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
196 }
197 PetscFunctionReturn(PETSC_SUCCESS);
198 }
199
200 /*MC
201 SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some
202 artificial $G(x)$ for which the `SNESFunctionFn` $F(x) = grad G(x)$. Therefore, this line search seeks
203 to find roots of the directional derivative via a secant method, that is $F(x_k - \lambda Y_k) \cdot Y_k / ||Y_k|| = 0$.
204
205 Options Database Keys:
206 + -snes_linesearch_minlambda <1e\-12> - the minimum acceptable `lambda` (scaling of solution update)
207 . -snes_linesearch_maxlambda <1.0> - the algorithm ensures that `lambda` is never larger than this value
208 . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search
209 . -snes_linesearch_order <1> - order of the approximation in the secant method, must be 1, 2, or 3
210 . -snes_linesearch_max_it <1> - the maximum number of secant iterations performed
211 . -snes_linesearch_rtol <1e\-8> - relative tolerance for the directional derivative
212 . -snes_linesearch_atol <1e\-15> - absolute tolerance for the directional derivative
213 - -snes_linesearch_ltol <1e\-8> - minimum absolute change in `lambda` allowed
214
215 Level: advanced
216
217 Notes:
218 This method does NOT use the objective function if it is provided with `SNESSetObjective()`.
219
220 This method is the preferred line search for `SNESQN` and `SNESNCG`.
221
222 .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHBISECTION`
223 M*/
SNESLineSearchCreate_CP(SNESLineSearch linesearch)224 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch)
225 {
226 PetscFunctionBegin;
227 linesearch->ops->apply = SNESLineSearchApply_CP;
228 linesearch->ops->destroy = NULL;
229 linesearch->ops->setfromoptions = NULL;
230 linesearch->ops->reset = NULL;
231 linesearch->ops->view = NULL;
232 linesearch->ops->setup = NULL;
233 linesearch->order = SNES_LINESEARCH_ORDER_LINEAR;
234
235 linesearch->max_it = 1;
236 PetscFunctionReturn(PETSC_SUCCESS);
237 }
238