1 #include <petsc/private/linesearchimpl.h> /*I "petscsnes.h" I*/
2 #include <petsc/private/snesimpl.h>
3
4 typedef struct {
5 PetscReal norm_delta_x_prev; /* norm of previous update */
6 PetscReal norm_bar_delta_x_prev; /* norm of previous bar update */
7 PetscReal mu_curr; /* current local Lipschitz estimate */
8 PetscReal lambda_prev; /* previous step length: for some reason SNESLineSearchGetLambda returns 1 instead of the previous step length */
9 } SNESLineSearch_NLEQERR;
10
11 static PetscBool NLEQERR_cited = PETSC_FALSE;
12 static const char NLEQERR_citation[] = "@book{deuflhard2011,\n"
13 " title = {Newton Methods for Nonlinear Problems},\n"
14 " author = {Peter Deuflhard},\n"
15 " volume = 35,\n"
16 " year = 2011,\n"
17 " isbn = {978-3-642-23898-7},\n"
18 " doi = {10.1007/978-3-642-23899-4},\n"
19 " publisher = {Springer-Verlag},\n"
20 " address = {Berlin, Heidelberg}\n}\n";
21
SNESLineSearchReset_NLEQERR(SNESLineSearch linesearch)22 static PetscErrorCode SNESLineSearchReset_NLEQERR(SNESLineSearch linesearch)
23 {
24 SNESLineSearch_NLEQERR *nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data;
25
26 PetscFunctionBegin;
27 nleqerr->mu_curr = 0.0;
28 nleqerr->norm_delta_x_prev = -1.0;
29 nleqerr->norm_bar_delta_x_prev = -1.0;
30 PetscFunctionReturn(PETSC_SUCCESS);
31 }
32
SNESLineSearchApply_NLEQERR(SNESLineSearch linesearch)33 static PetscErrorCode SNESLineSearchApply_NLEQERR(SNESLineSearch linesearch)
34 {
35 PetscBool changed_y, changed_w;
36 Vec X, F, Y, W, G;
37 SNES snes;
38 PetscReal fnorm, xnorm, ynorm, gnorm, wnorm;
39 PetscReal lambda, minlambda, stol;
40 PetscViewer monitor;
41 PetscInt max_it, count, snes_iteration;
42 PetscReal theta, mudash, lambdadash;
43 SNESLineSearch_NLEQERR *nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data;
44 KSPConvergedReason kspreason;
45
46 PetscFunctionBegin;
47 PetscCall(PetscCitationsRegister(NLEQERR_citation, &NLEQERR_cited));
48
49 PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G));
50 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm));
51 PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
52 PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
53 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
54 PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, NULL, NULL, NULL, NULL, &max_it));
55 PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL));
56
57 /* reset the state of the Lipschitz estimates */
58 PetscCall(SNESGetIterationNumber(snes, &snes_iteration));
59 if (!snes_iteration) PetscCall(SNESLineSearchReset_NLEQERR(linesearch));
60
61 /* precheck */
62 PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
63 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
64
65 PetscCall(VecNormBegin(Y, NORM_2, &ynorm));
66 PetscCall(VecNormBegin(X, NORM_2, &xnorm));
67 PetscCall(VecNormEnd(Y, NORM_2, &ynorm));
68 PetscCall(VecNormEnd(X, NORM_2, &xnorm));
69
70 /* Note: Y is *minus* the Newton step. For whatever reason PETSc doesn't solve with the minus on the RHS. */
71
72 if (ynorm == 0.0) {
73 if (monitor) {
74 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
75 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Initial direction and size is 0\n"));
76 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
77 }
78 PetscCall(VecCopy(X, W));
79 PetscCall(VecCopy(F, G));
80 PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
81 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
82 PetscFunctionReturn(PETSC_SUCCESS);
83 }
84
85 /* At this point, we've solved the Newton system for delta_x, and we assume that
86 its norm is greater than the solution tolerance (otherwise we wouldn't be in
87 here). So let's go ahead and estimate the Lipschitz constant.
88
89 W contains bar_delta_x_prev at this point. */
90
91 if (monitor) {
92 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
93 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: norm of Newton step: %14.12e\n", (double)ynorm));
94 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
95 }
96
97 /* this needs information from a previous iteration, so can't do it on the first one */
98 if (nleqerr->norm_delta_x_prev > 0 && nleqerr->norm_bar_delta_x_prev > 0) {
99 PetscCall(VecWAXPY(G, +1.0, Y, W)); /* bar_delta_x - delta_x; +1 because Y is -delta_x */
100 PetscCall(VecNormBegin(G, NORM_2, &gnorm));
101 PetscCall(VecNormEnd(G, NORM_2, &gnorm));
102
103 nleqerr->mu_curr = nleqerr->lambda_prev * (nleqerr->norm_delta_x_prev * nleqerr->norm_bar_delta_x_prev) / (gnorm * ynorm);
104 lambda = PetscMin(1.0, nleqerr->mu_curr);
105
106 if (monitor) {
107 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
108 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Lipschitz estimate: %14.12e; lambda: %14.12e\n", (double)nleqerr->mu_curr, (double)lambda));
109 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
110 }
111 } else {
112 lambda = linesearch->damping;
113 }
114
115 /* The main while loop of the algorithm.
116 At the end of this while loop, G should have the accepted new X in it. */
117
118 count = 0;
119 while (PETSC_TRUE) {
120 if (monitor) {
121 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
122 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: entering iteration with lambda: %14.12e\n", (double)lambda));
123 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
124 }
125
126 /* Check that we haven't performed too many iterations */
127 count += 1;
128 if (count >= max_it) {
129 if (monitor) {
130 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
131 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: maximum iterations reached\n"));
132 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
133 }
134 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
135 PetscFunctionReturn(PETSC_SUCCESS);
136 }
137
138 /* Now comes the Regularity Test. */
139 if (lambda <= minlambda) {
140 /* This isn't what is suggested by Deuflhard, but it works better in my experience */
141 if (monitor) {
142 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
143 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambda has reached minlambda, taking full Newton step\n"));
144 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
145 }
146 lambda = 1.0;
147 PetscCall(VecWAXPY(G, -lambda, Y, X));
148
149 /* and clean up the state for next time */
150 PetscCall(SNESLineSearchReset_NLEQERR(linesearch));
151 /*
152 The clang static analyzer detected a problem here; once the loop is broken the values
153 nleqerr->norm_delta_x_prev = ynorm;
154 nleqerr->norm_bar_delta_x_prev = wnorm;
155 are set, but wnorm has not even been computed.
156 I don't know if this is the correct fix but by setting ynorm and wnorm to -1.0 at
157 least the linesearch object is kept in the state set by the SNESLineSearchReset_NLEQERR() call above
158 */
159 ynorm = wnorm = -1.0;
160 break;
161 }
162
163 /* Compute new trial iterate */
164 PetscCall(VecWAXPY(W, -lambda, Y, X));
165 PetscCall(SNESComputeFunction(snes, W, G));
166
167 /* Solve linear system for bar_delta_x_curr: old Jacobian, new RHS. Note absence of minus sign, compared to Deuflhard, in keeping with PETSc convention */
168 PetscCall(KSPSolve(snes->ksp, G, W));
169 PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason));
170 if (kspreason < 0) PetscCall(PetscInfo(snes, "Solution for \\bar{delta x}^{k+1} failed.\n"));
171
172 /* W now contains -bar_delta_x_curr. */
173
174 PetscCall(VecNorm(W, NORM_2, &wnorm));
175 if (monitor) {
176 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
177 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: norm of simplified Newton update: %14.12e\n", (double)wnorm));
178 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
179 }
180
181 /* compute the monitoring quantities theta and mudash. */
182
183 theta = wnorm / ynorm;
184
185 PetscCall(VecWAXPY(G, -(1.0 - lambda), Y, W));
186 PetscCall(VecNorm(G, NORM_2, &gnorm));
187
188 mudash = (0.5 * ynorm * lambda * lambda) / gnorm;
189
190 /* Check for termination of the linesearch */
191 if (theta >= 1.0) {
192 /* need to go around again with smaller lambda */
193 if (monitor) {
194 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
195 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: monotonicity check failed, ratio: %14.12e\n", (double)theta));
196 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
197 }
198 lambda = PetscMin(mudash, 0.5 * lambda);
199 lambda = PetscMax(lambda, minlambda);
200 /* continue through the loop, i.e. go back to regularity test */
201 } else {
202 /* linesearch terminated */
203 lambdadash = PetscMin(1.0, mudash);
204
205 if (lambdadash == 1.0 && lambda == 1.0 && wnorm <= stol) {
206 /* store the updated state, X - Y - W, in G:
207 I need to keep W for the next linesearch */
208 PetscCall(VecWAXPY(G, -1.0, Y, X));
209 PetscCall(VecAXPY(G, -1.0, W));
210 break;
211 }
212
213 /* Deuflhard suggests to add the following:
214 else if (lambdadash >= 4.0 * lambda) {
215 lambda = lambdadash;
216 }
217 to continue through the loop, i.e. go back to regularity test.
218 I deliberately exclude this, as I have practical experience of this
219 getting stuck in infinite loops (on e.g. an Allen--Cahn problem). */
220
221 else {
222 /* accept iterate without adding on, i.e. don't use bar_delta_x;
223 again, I need to keep W for the next linesearch */
224 PetscCall(VecWAXPY(G, -lambda, Y, X));
225 break;
226 }
227 }
228 }
229
230 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, G));
231
232 /* W currently contains -bar_delta_u. Scale it so that it contains bar_delta_u. */
233 PetscCall(VecScale(W, -1.0));
234
235 /* postcheck */
236 PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, G, &changed_y, &changed_w));
237 if (changed_y || changed_w) {
238 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_USER));
239 PetscCall(PetscInfo(snes, "Changing the search direction here doesn't make sense.\n"));
240 PetscFunctionReturn(PETSC_SUCCESS);
241 }
242
243 /* copy the solution and information from this iteration over */
244 nleqerr->norm_delta_x_prev = ynorm;
245 nleqerr->norm_bar_delta_x_prev = wnorm;
246 nleqerr->lambda_prev = lambda;
247
248 PetscCall(VecCopy(G, X));
249 PetscCall(SNESComputeFunction(snes, X, F));
250 PetscCall(VecNorm(X, NORM_2, &xnorm));
251 PetscCall(VecNorm(F, NORM_2, &fnorm));
252 SNESLineSearchCheckFunctionDomainError(snes, linesearch, fnorm);
253 PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
254 PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm < 0 ? PETSC_INFINITY : ynorm));
255 PetscFunctionReturn(PETSC_SUCCESS);
256 }
257
SNESLineSearchView_NLEQERR(SNESLineSearch linesearch,PetscViewer viewer)258 static PetscErrorCode SNESLineSearchView_NLEQERR(SNESLineSearch linesearch, PetscViewer viewer)
259 {
260 PetscBool isascii;
261 SNESLineSearch_NLEQERR *nleqerr;
262
263 PetscFunctionBegin;
264 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
265 nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data;
266 if (isascii) {
267 PetscCall(PetscViewerASCIIPrintf(viewer, " NLEQ-ERR affine-covariant linesearch"));
268 PetscCall(PetscViewerASCIIPrintf(viewer, " current local Lipschitz estimate omega=%e\n", (double)nleqerr->mu_curr));
269 }
270 PetscFunctionReturn(PETSC_SUCCESS);
271 }
272
SNESLineSearchDestroy_NLEQERR(SNESLineSearch linesearch)273 static PetscErrorCode SNESLineSearchDestroy_NLEQERR(SNESLineSearch linesearch)
274 {
275 PetscFunctionBegin;
276 PetscCall(PetscFree(linesearch->data));
277 PetscFunctionReturn(PETSC_SUCCESS);
278 }
279
280 /*MC
281 SNESLINESEARCHNLEQERR - Error-oriented affine-covariant globalised Newton algorithm of Deuflhard {cite}`deuflhard2011`
282
283 This linesearch is intended for Newton-type methods which are affine covariant. Affine covariance
284 means that Newton's method will give the same iterations for F(x) = 0 and AF(x) = 0 for any nonsingular
285 matrix A. This is a fundamental property; the philosophy of this linesearch is that globalisations
286 of Newton's method should carefully preserve it.
287
288 Options Database Keys:
289 + -snes_linesearch_damping <1.0> - initial `lambda`
290 . -snes_linesearch_max_it <40> - maximum number of iterations for the line search
291 - -snes_linesearch_minlambda <1e\-12> - minimum `lambda` allowed
292
293 Level: advanced
294
295 Note:
296 Contributed by Patrick Farrell <patrick.farrell@maths.ox.ac.uk>
297
298 .seealso: [](ch_snes), `SNESLineSearch`, `SNES`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
299 M*/
SNESLineSearchCreate_NLEQERR(SNESLineSearch linesearch)300 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NLEQERR(SNESLineSearch linesearch)
301 {
302 SNESLineSearch_NLEQERR *nleqerr;
303
304 PetscFunctionBegin;
305 linesearch->ops->apply = SNESLineSearchApply_NLEQERR;
306 linesearch->ops->destroy = SNESLineSearchDestroy_NLEQERR;
307 linesearch->ops->setfromoptions = NULL;
308 linesearch->ops->reset = SNESLineSearchReset_NLEQERR;
309 linesearch->ops->view = SNESLineSearchView_NLEQERR;
310 linesearch->ops->setup = NULL;
311
312 PetscCall(PetscNew(&nleqerr));
313
314 linesearch->data = (void *)nleqerr;
315 linesearch->max_it = 40;
316 PetscCall(SNESLineSearchReset_NLEQERR(linesearch));
317 PetscFunctionReturn(PETSC_SUCCESS);
318 }
319