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