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