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:   PetscFunctionBegin;
     27:   nleqerr->mu_curr               = 0.0;
     28:   nleqerr->norm_delta_x_prev     = -1.0;
     29:   nleqerr->norm_bar_delta_x_prev = -1.0;
     30:   PetscFunctionReturn(PETSC_SUCCESS);
     31: }

     33: static PetscErrorCode SNESLineSearchApply_NLEQERR(SNESLineSearch linesearch)
     34: {
     35:   PetscBool               changed_y, changed_w;
     36:   Vec                     X, F, Y, W, G;
     37:   SNES                    snes;
     38:   PetscReal               fnorm, xnorm, ynorm, gnorm, wnorm;
     39:   PetscReal               lambda, minlambda, stol;
     40:   PetscViewer             monitor;
     41:   PetscInt                max_it, count, snes_iteration;
     42:   PetscReal               theta, mudash, lambdadash;
     43:   SNESLineSearch_NLEQERR *nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data;
     44:   KSPConvergedReason      kspreason;

     46:   PetscFunctionBegin;
     47:   PetscCall(PetscCitationsRegister(NLEQERR_citation, &NLEQERR_cited));

     49:   PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G));
     50:   PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm));
     51:   PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
     52:   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
     53:   PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
     54:   PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, NULL, NULL, NULL, NULL, &max_it));
     55:   PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL));

     57:   /* reset the state of the Lipschitz estimates */
     58:   PetscCall(SNESGetIterationNumber(snes, &snes_iteration));
     59:   if (!snes_iteration) PetscCall(SNESLineSearchReset_NLEQERR(linesearch));

     61:   /* precheck */
     62:   PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
     63:   PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));

     65:   PetscCall(VecNormBegin(Y, NORM_2, &ynorm));
     66:   PetscCall(VecNormBegin(X, NORM_2, &xnorm));
     67:   PetscCall(VecNormEnd(Y, NORM_2, &ynorm));
     68:   PetscCall(VecNormEnd(X, NORM_2, &xnorm));

     70:   /* Note: Y is *minus* the Newton step. For whatever reason PETSc doesn't solve with the minus on  the RHS. */

     72:   if (ynorm == 0.0) {
     73:     if (monitor) {
     74:       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     75:       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Initial direction and size is 0\n"));
     76:       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     77:     }
     78:     PetscCall(VecCopy(X, W));
     79:     PetscCall(VecCopy(F, G));
     80:     PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
     81:     PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
     82:     PetscFunctionReturn(PETSC_SUCCESS);
     83:   }

     85:   /* At this point, we've solved the Newton system for delta_x, and we assume that
     86:      its norm is greater than the solution tolerance (otherwise we wouldn't be in
     87:      here). So let's go ahead and estimate the Lipschitz constant.

     89:      W contains bar_delta_x_prev at this point. */

     91:   if (monitor) {
     92:     PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     93:     PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: norm of Newton step: %14.12e\n", (double)ynorm));
     94:     PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     95:   }

     97:   /* this needs information from a previous iteration, so can't do it on the first one */
     98:   if (nleqerr->norm_delta_x_prev > 0 && nleqerr->norm_bar_delta_x_prev > 0) {
     99:     PetscCall(VecWAXPY(G, +1.0, Y, W)); /* bar_delta_x - delta_x; +1 because Y is -delta_x */
    100:     PetscCall(VecNormBegin(G, NORM_2, &gnorm));
    101:     PetscCall(VecNormEnd(G, NORM_2, &gnorm));

    103:     nleqerr->mu_curr = nleqerr->lambda_prev * (nleqerr->norm_delta_x_prev * nleqerr->norm_bar_delta_x_prev) / (gnorm * ynorm);
    104:     lambda           = PetscMin(1.0, nleqerr->mu_curr);

    106:     if (monitor) {
    107:       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    108:       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Lipschitz estimate: %14.12e; lambda: %14.12e\n", (double)nleqerr->mu_curr, (double)lambda));
    109:       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    110:     }
    111:   } else {
    112:     lambda = linesearch->damping;
    113:   }

    115:   /* The main while loop of the algorithm.
    116:      At the end of this while loop, G should have the accepted new X in it. */

    118:   count = 0;
    119:   while (PETSC_TRUE) {
    120:     if (monitor) {
    121:       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    122:       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: entering iteration with lambda: %14.12e\n", (double)lambda));
    123:       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    124:     }

    126:     /* Check that we haven't performed too many iterations */
    127:     count += 1;
    128:     if (count >= max_it) {
    129:       if (monitor) {
    130:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    131:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: maximum iterations reached\n"));
    132:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    133:       }
    134:       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
    135:       PetscFunctionReturn(PETSC_SUCCESS);
    136:     }

    138:     /* Now comes the Regularity Test. */
    139:     if (lambda <= minlambda) {
    140:       /* This isn't what is suggested by Deuflhard, but it works better in my experience */
    141:       if (monitor) {
    142:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    143:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: lambda has reached minlambda, taking full Newton step\n"));
    144:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    145:       }
    146:       lambda = 1.0;
    147:       PetscCall(VecWAXPY(G, -lambda, Y, X));

    149:       /* and clean up the state for next time */
    150:       PetscCall(SNESLineSearchReset_NLEQERR(linesearch));
    151:       /*
    152:          The clang static analyzer detected a problem here; once the loop is broken the values
    153:          nleqerr->norm_delta_x_prev     = ynorm;
    154:          nleqerr->norm_bar_delta_x_prev = wnorm;
    155:          are set, but wnorm has not even been computed.
    156:          I don't know if this is the correct fix but by setting ynorm and wnorm to -1.0 at
    157:          least the linesearch object is kept in the state set by the SNESLineSearchReset_NLEQERR() call above
    158:       */
    159:       ynorm = wnorm = -1.0;
    160:       break;
    161:     }

    163:     /* Compute new trial iterate */
    164:     PetscCall(VecWAXPY(W, -lambda, Y, X));
    165:     PetscCall(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:     PetscCall(KSPSolve(snes->ksp, G, W));
    169:     PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason));
    170:     if (kspreason < 0) PetscCall(PetscInfo(snes, "Solution for \\bar{delta x}^{k+1} failed.\n"));

    172:     /* W now contains -bar_delta_x_curr. */

    174:     PetscCall(VecNorm(W, NORM_2, &wnorm));
    175:     if (monitor) {
    176:       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    177:       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: norm of simplified Newton update: %14.12e\n", (double)wnorm));
    178:       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    179:     }

    181:     /* compute the monitoring quantities theta and mudash. */

    183:     theta = wnorm / ynorm;

    185:     PetscCall(VecWAXPY(G, -(1.0 - lambda), Y, W));
    186:     PetscCall(VecNorm(G, NORM_2, &gnorm));

    188:     mudash = (0.5 * ynorm * lambda * lambda) / gnorm;

    190:     /* Check for termination of the linesearch */
    191:     if (theta >= 1.0) {
    192:       /* need to go around again with smaller lambda */
    193:       if (monitor) {
    194:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    195:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: monotonicity check failed, ratio: %14.12e\n", (double)theta));
    196:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    197:       }
    198:       lambda = PetscMin(mudash, 0.5 * lambda);
    199:       lambda = PetscMax(lambda, minlambda);
    200:       /* continue through the loop, i.e. go back to regularity test */
    201:     } else {
    202:       /* linesearch terminated */
    203:       lambdadash = PetscMin(1.0, mudash);

    205:       if (lambdadash == 1.0 && lambda == 1.0 && wnorm <= stol) {
    206:         /* store the updated state, X - Y - W, in G:
    207:            I need to keep W for the next linesearch */
    208:         PetscCall(VecWAXPY(G, -1.0, Y, X));
    209:         PetscCall(VecAXPY(G, -1.0, W));
    210:         break;
    211:       }

    213:       /* Deuflhard suggests to add the following:
    214:       else if (lambdadash >= 4.0 * lambda) {
    215:         lambda = lambdadash;
    216:       }
    217:       to continue through the loop, i.e. go back to regularity test.
    218:       I deliberately exclude this, as I have practical experience of this
    219:       getting stuck in infinite loops (on e.g. an Allen--Cahn problem). */

    221:       else {
    222:         /* accept iterate without adding on, i.e. don't use bar_delta_x;
    223:            again, I need to keep W for the next linesearch */
    224:         PetscCall(VecWAXPY(G, -lambda, Y, X));
    225:         break;
    226:       }
    227:     }
    228:   }

    230:   if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, G));

    232:   /* W currently contains -bar_delta_u. Scale it so that it contains bar_delta_u. */
    233:   PetscCall(VecScale(W, -1.0));

    235:   /* postcheck */
    236:   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, G, &changed_y, &changed_w));
    237:   if (changed_y || changed_w) {
    238:     PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_USER));
    239:     PetscCall(PetscInfo(snes, "Changing the search direction here doesn't make sense.\n"));
    240:     PetscFunctionReturn(PETSC_SUCCESS);
    241:   }

    243:   /* copy the solution and information from this iteration over */
    244:   nleqerr->norm_delta_x_prev     = ynorm;
    245:   nleqerr->norm_bar_delta_x_prev = wnorm;
    246:   nleqerr->lambda_prev           = lambda;

    248:   PetscCall(VecCopy(G, X));
    249:   PetscCall(SNESComputeFunction(snes, X, F));
    250:   PetscCall(VecNorm(X, NORM_2, &xnorm));
    251:   PetscCall(VecNorm(F, NORM_2, &fnorm));
    252:   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
    253:   PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm < 0 ? PETSC_INFINITY : ynorm));
    254:   PetscFunctionReturn(PETSC_SUCCESS);
    255: }

    257: static PetscErrorCode SNESLineSearchView_NLEQERR(SNESLineSearch linesearch, PetscViewer viewer)
    258: {
    259:   PetscBool               isascii;
    260:   SNESLineSearch_NLEQERR *nleqerr;

    262:   PetscFunctionBegin;
    263:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    264:   nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data;
    265:   if (isascii) {
    266:     PetscCall(PetscViewerASCIIPrintf(viewer, "  NLEQ-ERR affine-covariant linesearch"));
    267:     PetscCall(PetscViewerASCIIPrintf(viewer, "  current local Lipschitz estimate omega=%e\n", (double)nleqerr->mu_curr));
    268:   }
    269:   PetscFunctionReturn(PETSC_SUCCESS);
    270: }

    272: static PetscErrorCode SNESLineSearchDestroy_NLEQERR(SNESLineSearch linesearch)
    273: {
    274:   PetscFunctionBegin;
    275:   PetscCall(PetscFree(linesearch->data));
    276:   PetscFunctionReturn(PETSC_SUCCESS);
    277: }

    279: /*MC
    280:    SNESLINESEARCHNLEQERR - Error-oriented affine-covariant globalised Newton algorithm of Deuflhard {cite}`deuflhard2011`

    282:    This linesearch is intended for Newton-type methods which are affine covariant. Affine covariance
    283:    means that Newton's method will give the same iterations for F(x) = 0 and AF(x) = 0 for any nonsingular
    284:    matrix A. This is a fundamental property; the philosophy of this linesearch is that globalisations
    285:    of Newton's method should carefully preserve it.

    287:    Options Database Keys:
    288: +  -snes_linesearch_damping <1.0>      - initial `lambda`
    289: .  -snes_linesearch_max_it <40>        - maximum number of iterations for the line search
    290: -  -snes_linesearch_minlambda <1e\-12> - minimum `lambda` allowed

    292:    Level: advanced

    294:    Note:
    295:    Contributed by Patrick Farrell <patrick.farrell@maths.ox.ac.uk>

    297: .seealso: [](ch_snes), `SNESLineSearch`, `SNES`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
    298: M*/
    299: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NLEQERR(SNESLineSearch linesearch)
    300: {
    301:   SNESLineSearch_NLEQERR *nleqerr;

    303:   PetscFunctionBegin;
    304:   linesearch->ops->apply          = SNESLineSearchApply_NLEQERR;
    305:   linesearch->ops->destroy        = SNESLineSearchDestroy_NLEQERR;
    306:   linesearch->ops->setfromoptions = NULL;
    307:   linesearch->ops->reset          = SNESLineSearchReset_NLEQERR;
    308:   linesearch->ops->view           = SNESLineSearchView_NLEQERR;
    309:   linesearch->ops->setup          = NULL;

    311:   PetscCall(PetscNew(&nleqerr));

    313:   linesearch->data   = (void *)nleqerr;
    314:   linesearch->max_it = 40;
    315:   PetscCall(SNESLineSearchReset_NLEQERR(linesearch));
    316:   PetscFunctionReturn(PETSC_SUCCESS);
    317: }