Actual source code: linesearchnleqerr.c
petsc-3.6.1 2015-08-06
1: #include <petsc/private/linesearchimpl.h> /*I "petscsnes.h" I*/
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";
24: static PetscErrorCode SNESLineSearchReset_NLEQERR(SNESLineSearch linesearch)
25: {
26: SNESLineSearch_NLEQERR *nleqerr;
28: nleqerr = (SNESLineSearch_NLEQERR*)linesearch->data;
29: nleqerr->mu_curr = 0.0;
30: nleqerr->norm_delta_x_prev = -1.0;
31: nleqerr->norm_bar_delta_x_prev = -1.0;
32: return(0);
33: }
37: static PetscErrorCode SNESLineSearchApply_NLEQERR(SNESLineSearch linesearch)
38: {
39: PetscBool changed_y,changed_w;
40: PetscErrorCode ierr;
41: Vec X,F,Y,W,G;
42: SNES snes;
43: PetscReal fnorm, xnorm, ynorm, gnorm, wnorm;
44: PetscReal lambda, minlambda, stol;
45: PetscViewer monitor;
46: PetscInt max_its, count, snes_iteration;
47: PetscReal theta, mudash, lambdadash;
48: SNESLineSearch_NLEQERR *nleqerr;
49: KSPConvergedReason kspreason;
52: PetscCitationsRegister(NLEQERR_citation, &NLEQERR_cited);
54: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);
55: SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);
56: SNESLineSearchGetLambda(linesearch, &lambda);
57: SNESLineSearchGetSNES(linesearch, &snes);
58: SNESLineSearchGetMonitor(linesearch, &monitor);
59: SNESLineSearchGetTolerances(linesearch,&minlambda,NULL,NULL,NULL,NULL,&max_its);
60: SNESGetTolerances(snes,NULL,NULL,&stol,NULL,NULL);
61: nleqerr = (SNESLineSearch_NLEQERR*)linesearch->data;
63: /* reset the state of the Lipschitz estimates */
64: SNESGetIterationNumber(snes, &snes_iteration);
65: if (snes_iteration == 0) {
66: SNESLineSearchReset_NLEQERR(linesearch);
67: }
69: /* precheck */
70: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
71: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
73: VecNormBegin(Y, NORM_2, &ynorm);
74: VecNormBegin(X, NORM_2, &xnorm);
75: VecNormEnd(Y, NORM_2, &ynorm);
76: VecNormEnd(X, NORM_2, &xnorm);
78: /* Note: Y is *minus* the Newton step. For whatever reason PETSc doesn't solve with the minus on
79: the RHS. */
81: if (ynorm == 0.0) {
82: if (monitor) {
83: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
84: PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n");
85: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
86: }
87: VecCopy(X,W);
88: VecCopy(F,G);
89: SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);
90: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
91: return(0);
92: }
94: /* At this point, we've solved the Newton system for delta_x, and we assume that
95: its norm is greater than the solution tolerance (otherwise we wouldn't be in
96: here). So let's go ahead and estimate the Lipschitz constant.
98: W contains bar_delta_x_prev at this point. */
100: if (monitor) {
101: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
102: PetscViewerASCIIPrintf(monitor," Line search: norm of Newton step: %14.12e\n", (double) ynorm);
103: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
104: }
106: /* this needs information from a previous iteration, so can't do it on the first one */
107: if (nleqerr->norm_delta_x_prev > 0 && nleqerr->norm_bar_delta_x_prev > 0) {
108: VecWAXPY(G, +1.0, Y, W); /* bar_delta_x - delta_x; +1 because Y is -delta_x */
109: VecNormBegin(G, NORM_2, &gnorm);
110: VecNormEnd(G, NORM_2, &gnorm);
112: nleqerr->mu_curr = nleqerr->lambda_prev * (nleqerr->norm_delta_x_prev * nleqerr->norm_bar_delta_x_prev) / (gnorm * ynorm);
113: lambda = PetscMin(1.0, nleqerr->mu_curr);
115: if (monitor) {
116: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
117: PetscViewerASCIIPrintf(monitor," Line search: Lipschitz estimate: %14.12e; lambda: %14.12e\n", (double) nleqerr->mu_curr, (double) lambda);
118: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
119: }
120: } else {
121: lambda = linesearch->damping;
122: }
124: /* The main while loop of the algorithm.
125: At the end of this while loop, G should have the accepted new X in it. */
127: count = 0;
128: while (PETSC_TRUE) {
129: if (monitor) {
130: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
131: PetscViewerASCIIPrintf(monitor," Line search: entering iteration with lambda: %14.12e\n", lambda);
132: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
133: }
135: /* Check that we haven't performed too many iterations */
136: count += 1;
137: if (count >= max_its) {
138: if (monitor) {
139: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
140: PetscViewerASCIIPrintf(monitor," Line search: maximum iterations reached\n");
141: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
142: }
143: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
144: return(0);
145: }
147: /* Now comes the Regularity Test. */
148: if (lambda <= minlambda) {
149: /* This isn't what is suggested by Deuflhard, but it works better in my experience */
150: if (monitor) {
151: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
152: PetscViewerASCIIPrintf(monitor," Line search: lambda has reached lambdamin, taking full Newton step\n");
153: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
154: }
155: lambda = 1.0;
156: VecWAXPY(G, -lambda, Y, X);
158: /* and clean up the state for next time */
159: SNESLineSearchReset_NLEQERR(linesearch);
160: break;
161: }
163: /* Compute new trial iterate */
164: VecWAXPY(W, -lambda, Y, X);
165: 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: KSPSolve(snes->ksp, G, W);
169: KSPGetConvergedReason(snes->ksp, &kspreason);
170: if (kspreason < 0) {
171: PetscInfo(snes,"Solution for \\bar{delta x}^{k+1} failed.");
172: }
174: /* W now contains -bar_delta_x_curr. */
176: VecNorm(W, NORM_2, &wnorm);
177: if (monitor) {
178: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
179: PetscViewerASCIIPrintf(monitor," Line search: norm of simplified Newton update: %14.12e\n", (double) wnorm);
180: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
181: }
183: /* compute the monitoring quantities theta and mudash. */
185: theta = wnorm / ynorm;
187: VecWAXPY(G, -(1.0 - lambda), Y, W);
188: VecNorm(G, NORM_2, &gnorm);
190: mudash = (0.5 * ynorm * lambda * lambda) / gnorm;
192: /* Check for termination of the linesearch */
193: if (theta >= 1.0) {
194: /* need to go around again with smaller lambda */
195: if (monitor) {
196: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
197: PetscViewerASCIIPrintf(monitor," Line search: monotonicity check failed, ratio: %14.12e\n", (double) theta);
198: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
199: }
200: lambda = PetscMin(mudash, 0.5 * lambda);
201: lambda = PetscMax(lambda, minlambda);
202: /* continue through the loop, i.e. go back to regularity test */
203: } else {
204: /* linesearch terminated */
205: lambdadash = PetscMin(1.0, mudash);
207: if (lambdadash == 1.0 && lambda == 1.0 && wnorm <= stol) {
208: /* store the updated state, X - Y - W, in G:
209: I need to keep W for the next linesearch */
210: VecCopy(X, G);
211: VecAXPY(G, -1.0, Y);
212: VecAXPY(G, -1.0, W);
213: break;
214: }
216: /* Deuflhard suggests to add the following:
217: else if (lambdadash >= 4.0 * lambda) {
218: lambda = lambdadash;
219: }
220: to continue through the loop, i.e. go back to regularity test.
221: I deliberately exclude this, as I have practical experience of this
222: getting stuck in infinite loops (on e.g. an Allen--Cahn problem). */
224: else {
225: /* accept iterate without adding on, i.e. don't use bar_delta_x;
226: again, I need to keep W for the next linesearch */
227: VecWAXPY(G, -lambda, Y, X);
228: break;
229: }
230: }
231: }
233: if (linesearch->ops->viproject) {
234: (*linesearch->ops->viproject)(snes, G);
235: }
237: /* W currently contains -bar_delta_u. Scale it so that it contains bar_delta_u. */
238: VecScale(W, -1.0);
240: /* postcheck */
241: SNESLineSearchPostCheck(linesearch,X,Y,G,&changed_y,&changed_w);
242: if (changed_y || changed_w) {
243: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_USER);
244: PetscInfo(snes,"Changing the search direction here doesn't make sense.\n");
245: return(0);
246: }
248: /* copy the solution and information from this iteration over */
249: nleqerr->norm_delta_x_prev = ynorm;
250: nleqerr->norm_bar_delta_x_prev = wnorm;
251: nleqerr->lambda_prev = lambda;
253: VecCopy(G, X);
254: SNESComputeFunction(snes, X, F);
255: VecNorm(X, NORM_2, &xnorm);
256: VecNorm(F, NORM_2, &fnorm);
257: SNESLineSearchSetLambda(linesearch, lambda);
258: SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm);
259: return(0);
260: }
264: PetscErrorCode SNESLineSearchView_NLEQERR(SNESLineSearch linesearch, PetscViewer viewer)
265: {
266: PetscErrorCode ierr;
267: PetscBool iascii;
268: SNESLineSearch_NLEQERR *nleqerr;
271: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
272: nleqerr = (SNESLineSearch_NLEQERR*)linesearch->data;
273: if (iascii) {
274: PetscViewerASCIIPrintf(viewer, " NLEQ-ERR affine-covariant linesearch");
275: PetscViewerASCIIPrintf(viewer, " current local Lipschitz estimate omega=%e\n", (double)nleqerr->mu_curr);
276: }
277: return(0);
278: }
282: static PetscErrorCode SNESLineSearchDestroy_NLEQERR(SNESLineSearch linesearch)
283: {
287: PetscFree(linesearch->data);
288: return(0);
289: }
293: /*MC
294: SNESLINESEARCHNLEQERR - Error-oriented affine-covariant globalised Newton algorithm of Deuflhard (2011).
296: This linesearch is intended for Newton-type methods which are affine covariant. Affine covariance
297: means that Newton's method will give the same iterations for F(x) = 0 and AF(x) = 0 for a nonsingular
298: matrix A. This is a fundamental property; the philosophy of this linesearch is that globalisations
299: of Newton's method should carefully preserve it.
301: For a discussion of the theory behind this algorithm, see
303: @book{deuflhard2011,
304: title={Newton Methods for Nonlinear Problems},
305: author={Deuflhard, P.},
306: volume={35},
307: year={2011},
308: publisher={Springer-Verlag},
309: address={Berlin, Heidelberg}
310: }
312: Pseudocode is given on page 148.
314: Options Database Keys:
315: + -snes_linesearch_damping<1.0> - initial step length
316: - -snes_linesearch_minlambda<1e-12> - minimum step length allowed
318: Contributed by Patrick Farrell <patrick.farrell@maths.ox.ac.uk>
320: Level: advanced
322: .keywords: SNES, SNESLineSearch, damping
324: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
325: M*/
326: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NLEQERR(SNESLineSearch linesearch)
327: {
329: SNESLineSearch_NLEQERR *nleqerr;
330: PetscErrorCode ierr;
333: linesearch->ops->apply = SNESLineSearchApply_NLEQERR;
334: linesearch->ops->destroy = SNESLineSearchDestroy_NLEQERR;
335: linesearch->ops->setfromoptions = NULL;
336: linesearch->ops->reset = SNESLineSearchReset_NLEQERR;
337: linesearch->ops->view = SNESLineSearchView_NLEQERR;
338: linesearch->ops->setup = NULL;
340: PetscNewLog(linesearch,&nleqerr);
342: linesearch->data = (void*)nleqerr;
343: linesearch->max_its = 40;
344: SNESLineSearchReset_NLEQERR(linesearch);
345: return(0);
346: }