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: }