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