Actual source code: linesearchcp.c

      1: #include <petsc/private/linesearchimpl.h>
      2: #include <petscsnes.h>

      4: static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch)
      5: {
      6:   PetscBool   changed_y, changed_w;
      7:   Vec         X, Y, F, W;
      8:   SNES        snes;
      9:   PetscReal   xnorm, ynorm, gnorm, minlambda, maxlambda, rtol, atol, ltol;
     10:   PetscReal   lambda, lambda_old, lambda_update, delLambda;
     11:   PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s;
     12:   PetscInt    i, max_it;
     13:   PetscViewer monitor;

     15:   PetscFunctionBegin;
     16:   PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL));
     17:   PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
     18:   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
     19:   PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
     20:   PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxlambda, &rtol, &atol, &ltol, &max_it));
     21:   PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
     22:   PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));

     24:   /* precheck */
     25:   PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
     26:   lambda_old = 0.0;

     28:   /* evaluate initial point */
     29:   if (linesearch->ops->vidirderiv) {
     30:     PetscCall((*linesearch->ops->vidirderiv)(snes, F, X, Y, &fty_old));
     31:   } else {
     32:     PetscCall(VecDot(F, Y, &fty_old));
     33:   }
     34:   if (PetscAbsScalar(fty_old) < atol * ynorm) {
     35:     if (monitor) {
     36:       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     37:       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search terminated at initial point because dot(F,Y) = %g < atol*||y|| = %g\n", (double)PetscAbsScalar(fty_old), (double)(atol * ynorm)));
     38:       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     39:     }
     40:     PetscCall(SNESSetConvergedReason(linesearch->snes, SNES_CONVERGED_FNORM_ABS));
     41:     PetscFunctionReturn(PETSC_SUCCESS);
     42:   }
     43:   fty_init = fty_old;

     45:   for (i = 0; i < max_it; i++) {
     46:     /* compute the norm at lambda */
     47:     PetscCall(VecWAXPY(W, -lambda, Y, X));
     48:     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
     49:     PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
     50:     if (linesearch->ops->vidirderiv) {
     51:       PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty));
     52:     } else {
     53:       PetscCall(VecDot(F, Y, &fty));
     54:     }

     56:     /* compute change of lambda */
     57:     delLambda = lambda - lambda_old;

     59:     /* check change of lambda tolerance */
     60:     if (PetscAbsReal(delLambda) < ltol) {
     61:       if (monitor) {
     62:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     63:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: abs(dlambda) = %g < ltol = %g\n", (double)PetscAbsReal(delLambda), (double)ltol));
     64:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     65:       }
     66:       break;
     67:     }

     69:     /* check relative tolerance */
     70:     if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) {
     71:       if (monitor) {
     72:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     73:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: abs(fty/fty_init) = %g <= rtol  = %g\n", (double)(PetscAbsScalar(fty) / PetscAbsScalar(fty_init)), (double)rtol));
     74:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     75:       }
     76:       break;
     77:     }

     79:     /* check absolute tolerance */
     80:     if (PetscAbsScalar(fty) < atol * ynorm && i > 0) {
     81:       if (monitor) {
     82:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     83:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: abs(fty)/||y|| = %g <= atol = %g\n", (double)(PetscAbsScalar(fty) / ynorm), (double)atol));
     84:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     85:       }
     86:       break;
     87:     }

     89:     /* print iteration information */
     90:     if (monitor) {
     91:       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     92:       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: lambdas = [%g, %g], ftys = [%g, %g]\n", (double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old)));
     93:       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     94:     }

     96:     /* approximate the second derivative */
     97:     if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) {
     98:       /* first order finite difference approximation */
     99:       s = (fty - fty_old) / delLambda;
    100:     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
    101:       /* evaluate function at midpoint 0.5 * (lambda + lambda_old) */
    102:       PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X));
    103:       if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
    104:       PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
    105:       if (linesearch->ops->vidirderiv) {
    106:         PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1));
    107:       } else {
    108:         PetscCall(VecDot(F, Y, &fty_mid1));
    109:       }

    111:       /* second order finite difference approximation */
    112:       s = (3. * fty - 4. * fty_mid1 + fty_old) / delLambda;

    114:     } else {
    115:       /* evaluate function at midpoint 0.5 * (lambda + lambda_old) */
    116:       PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X));
    117:       if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
    118:       PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
    119:       if (linesearch->ops->vidirderiv) {
    120:         PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1));
    121:       } else {
    122:         PetscCall(VecDot(F, Y, &fty_mid1));
    123:       }

    125:       /* evaluate function at lambda + 0.5 * (lambda - lambda_old) */
    126:       PetscCall(VecWAXPY(W, -(lambda + 0.5 * (lambda - lambda_old)), Y, X));
    127:       if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
    128:       PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
    129:       if (linesearch->ops->vidirderiv) {
    130:         PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid2));
    131:       } else {
    132:         PetscCall(VecDot(F, Y, &fty_mid2));
    133:       }

    135:       /* third order finite difference approximation */
    136:       s = (2. * fty_mid2 + 3. * fty - 6. * fty_mid1 + fty_old) / (3. * delLambda);
    137:     }

    139:     /* compute secant update (modifying the search direction if necessary) */
    140:     if (PetscRealPart(s) > 0.) s = -s;
    141:     if (s == 0.0) {
    142:       if (monitor) {
    143:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    144:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: del2Fnrm = 0\n"));
    145:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    146:       }
    147:       break;
    148:     }
    149:     lambda_update = lambda - PetscRealPart(fty / s);

    151:     /* if step is too small, go the opposite direction */
    152:     if (lambda_update < minlambda) lambda_update = lambda + PetscRealPart(fty / s);
    153:     /* if secant method would step out of bounds, exit with the respective bound */
    154:     if (lambda_update > maxlambda) {
    155:       lambda = maxlambda;
    156:       break;
    157:     }

    159:     /* if lambda is NaN or Inf, do not accept update but exit with previous lambda */
    160:     if (PetscIsInfOrNanReal(lambda_update)) {
    161:       if (monitor) {
    162:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    163:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: lambda_update is NaN or Inf\n"));
    164:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    165:       }
    166:       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
    167:       break;
    168:     }

    170:     /* compute the new state of the line search */
    171:     lambda_old = lambda;
    172:     lambda     = lambda_update;
    173:     fty_old    = fty;
    174:   }

    176:   /* construct the solution */
    177:   PetscCall(VecWAXPY(W, -lambda, Y, X));
    178:   if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
    179:   /* postcheck */
    180:   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
    181:   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
    182:   if (changed_y) {
    183:     if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
    184:     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
    185:   }
    186:   PetscCall(VecCopy(W, X));
    187:   PetscCall((*linesearch->ops->snesfunc)(snes, X, F));

    189:   PetscCall(SNESLineSearchComputeNorms(linesearch));
    190:   PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));

    192:   if (monitor) {
    193:     PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    194:     PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm));
    195:     PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    196:   }
    197:   PetscFunctionReturn(PETSC_SUCCESS);
    198: }

    200: /*MC
    201:    SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some
    202:    artificial $G(x)$ for which the `SNESFunctionFn` $F(x) = grad G(x)$. Therefore, this line search seeks
    203:    to find roots of the directional derivative via a secant method, that is $F(x_k - \lambda Y_k) \cdot Y_k / ||Y_k|| = 0$.

    205:    Options Database Keys:
    206: +  -snes_linesearch_minlambda <1e\-12> - the minimum acceptable `lambda` (scaling of solution update)
    207: .  -snes_linesearch_maxlambda <1.0>    - the algorithm ensures that `lambda` is never larger than this value
    208: .  -snes_linesearch_damping <1.0>      - initial `lambda` on entry to the line search
    209: .  -snes_linesearch_order <1>          - order of the approximation in the secant method, must be 1, 2, or 3
    210: .  -snes_linesearch_max_it <1>         - the maximum number of secant iterations performed
    211: .  -snes_linesearch_rtol <1e\-8>       - relative tolerance for the directional derivative
    212: .  -snes_linesearch_atol <1e\-15>      - absolute tolerance for the directional derivative
    213: -  -snes_linesearch_ltol <1e\-8>       - minimum absolute change in `lambda` allowed

    215:    Level: advanced

    217:    Notes:
    218:    This method does NOT use the objective function if it is provided with `SNESSetObjective()`.

    220:    This method is the preferred line search for `SNESQN` and `SNESNCG`.

    222: .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHBISECTION`
    223: M*/
    224: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch)
    225: {
    226:   PetscFunctionBegin;
    227:   linesearch->ops->apply          = SNESLineSearchApply_CP;
    228:   linesearch->ops->destroy        = NULL;
    229:   linesearch->ops->setfromoptions = NULL;
    230:   linesearch->ops->reset          = NULL;
    231:   linesearch->ops->view           = NULL;
    232:   linesearch->ops->setup          = NULL;
    233:   linesearch->order               = SNES_LINESEARCH_ORDER_LINEAR;

    235:   linesearch->max_it = 1;
    236:   PetscFunctionReturn(PETSC_SUCCESS);
    237: }