Actual source code: linesearchnleqerr.c
1: #include <petsc/private/linesearchimpl.h>
2: #include <petsc/private/snesimpl.h>
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;
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";
22: static PetscErrorCode SNESLineSearchReset_NLEQERR(SNESLineSearch linesearch)
23: {
24: SNESLineSearch_NLEQERR *nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data;
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: }
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_its, count, snes_iteration;
42: PetscReal theta, mudash, lambdadash;
43: SNESLineSearch_NLEQERR *nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data;
44: KSPConvergedReason kspreason;
46: PetscFunctionBegin;
47: PetscCall(PetscCitationsRegister(NLEQERR_citation, &NLEQERR_cited));
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_its));
55: PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL));
57: /* reset the state of the Lipschitz estimates */
58: PetscCall(SNESGetIterationNumber(snes, &snes_iteration));
59: if (!snes_iteration) PetscCall(SNESLineSearchReset_NLEQERR(linesearch));
61: /* precheck */
62: PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
63: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
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));
70: /* Note: Y is *minus* the Newton step. For whatever reason PETSc doesn't solve with the minus on the RHS. */
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: }
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.
89: W contains bar_delta_x_prev at this point. */
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: }
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));
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);
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: }
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. */
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: }
126: /* Check that we haven't performed too many iterations */
127: count += 1;
128: if (count >= max_its) {
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: }
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 lambdamin, taking full Newton step\n"));
144: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
145: }
146: lambda = 1.0;
147: PetscCall(VecWAXPY(G, -lambda, Y, X));
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: }
163: /* Compute new trial iterate */
164: PetscCall(VecWAXPY(W, -lambda, Y, X));
165: PetscCall(SNESComputeFunction(snes, W, G));
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"));
172: /* W now contains -bar_delta_x_curr. */
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: }
181: /* compute the monitoring quantities theta and mudash. */
183: theta = wnorm / ynorm;
185: PetscCall(VecWAXPY(G, -(1.0 - lambda), Y, W));
186: PetscCall(VecNorm(G, NORM_2, &gnorm));
188: mudash = (0.5 * ynorm * lambda * lambda) / gnorm;
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);
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: }
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). */
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: }
230: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, G));
232: /* W currently contains -bar_delta_u. Scale it so that it contains bar_delta_u. */
233: PetscCall(VecScale(W, -1.0));
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: }
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;
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: PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
253: PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm < 0 ? PETSC_INFINITY : ynorm));
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: static PetscErrorCode SNESLineSearchView_NLEQERR(SNESLineSearch linesearch, PetscViewer viewer)
258: {
259: PetscBool iascii;
260: SNESLineSearch_NLEQERR *nleqerr;
262: PetscFunctionBegin;
263: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
264: nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data;
265: if (iascii) {
266: PetscCall(PetscViewerASCIIPrintf(viewer, " NLEQ-ERR affine-covariant linesearch"));
267: PetscCall(PetscViewerASCIIPrintf(viewer, " current local Lipschitz estimate omega=%e\n", (double)nleqerr->mu_curr));
268: }
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode SNESLineSearchDestroy_NLEQERR(SNESLineSearch linesearch)
273: {
274: PetscFunctionBegin;
275: PetscCall(PetscFree(linesearch->data));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*MC
280: SNESLINESEARCHNLEQERR - Error-oriented affine-covariant globalised Newton algorithm of Deuflhard {cite}`deuflhard2011`
282: This linesearch is intended for Newton-type methods which are affine covariant. Affine covariance
283: means that Newton's method will give the same iterations for F(x) = 0 and AF(x) = 0 for any nonsingular
284: matrix A. This is a fundamental property; the philosophy of this linesearch is that globalisations
285: of Newton's method should carefully preserve it.
287: Options Database Keys:
288: + -snes_linesearch_damping<1.0> - initial step length
289: - -snes_linesearch_minlambda<1e-12> - minimum step length allowed
291: Level: advanced
293: Note:
294: Contributed by Patrick Farrell <patrick.farrell@maths.ox.ac.uk>
296: .seealso: [](ch_snes), `SNESLineSearch`, `SNES`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
297: M*/
298: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NLEQERR(SNESLineSearch linesearch)
299: {
300: SNESLineSearch_NLEQERR *nleqerr;
302: PetscFunctionBegin;
303: linesearch->ops->apply = SNESLineSearchApply_NLEQERR;
304: linesearch->ops->destroy = SNESLineSearchDestroy_NLEQERR;
305: linesearch->ops->setfromoptions = NULL;
306: linesearch->ops->reset = SNESLineSearchReset_NLEQERR;
307: linesearch->ops->view = SNESLineSearchView_NLEQERR;
308: linesearch->ops->setup = NULL;
310: PetscCall(PetscNew(&nleqerr));
312: linesearch->data = (void *)nleqerr;
313: linesearch->max_its = 40;
314: PetscCall(SNESLineSearchReset_NLEQERR(linesearch));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }