Actual source code: linesearchl2.c

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

      4: static PetscErrorCode SNESLineSearchApply_L2(SNESLineSearch linesearch)
      5: {
      6:   PetscBool        changed_y, changed_w;
      7:   Vec              X;
      8:   Vec              F;
      9:   Vec              Y;
     10:   Vec              W;
     11:   SNES             snes;
     12:   PetscReal        gnorm;
     13:   PetscReal        ynorm;
     14:   PetscReal        xnorm;
     15:   PetscReal        steptol, maxstep, rtol, atol, ltol;
     16:   PetscViewer      monitor;
     17:   PetscReal        lambda, lambda_old, lambda_mid, lambda_update, delLambda;
     18:   PetscReal        fnrm, fnrm_old, fnrm_mid;
     19:   PetscReal        delFnrm, delFnrm_old, del2Fnrm;
     20:   PetscInt         i, max_its;
     21:   SNESObjectiveFn *objective;

     23:   PetscFunctionBegin;
     24:   PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL));
     25:   PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
     26:   PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
     27:   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
     28:   PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
     29:   PetscCall(SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, &ltol, &max_its));
     30:   PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));

     32:   PetscCall(SNESGetObjective(snes, &objective, NULL));

     34:   /* precheck */
     35:   PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
     36:   lambda_old = 0.0;
     37:   if (!objective) {
     38:     fnrm_old = gnorm * gnorm;
     39:   } else {
     40:     PetscCall(SNESComputeObjective(snes, X, &fnrm_old));
     41:   }
     42:   lambda_mid = 0.5 * (lambda + lambda_old);

     44:   for (i = 0; i < max_its; i++) {
     45:     while (PETSC_TRUE) {
     46:       PetscCall(VecWAXPY(W, -lambda_mid, Y, X));
     47:       if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
     48:       if (!objective) {
     49:         /* compute the norm at the midpoint */
     50:         PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
     51:         if (linesearch->ops->vinorm) {
     52:           fnrm_mid = gnorm;
     53:           PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid));
     54:         } else {
     55:           PetscCall(VecNorm(F, NORM_2, &fnrm_mid));
     56:         }

     58:         /* compute the norm at the new endpoint */
     59:         PetscCall(VecWAXPY(W, -lambda, Y, X));
     60:         if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
     61:         PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
     62:         if (linesearch->ops->vinorm) {
     63:           fnrm = gnorm;
     64:           PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm));
     65:         } else {
     66:           PetscCall(VecNorm(F, NORM_2, &fnrm));
     67:         }
     68:         fnrm_mid = fnrm_mid * fnrm_mid;
     69:         fnrm     = fnrm * fnrm;
     70:       } else {
     71:         /* compute the objective at the midpoint */
     72:         PetscCall(SNESComputeObjective(snes, W, &fnrm_mid));

     74:         /* compute the objective at the new endpoint */
     75:         PetscCall(VecWAXPY(W, -lambda, Y, X));
     76:         PetscCall(SNESComputeObjective(snes, W, &fnrm));
     77:       }
     78:       if (!PetscIsInfOrNanReal(fnrm)) break;
     79:       if (monitor) {
     80:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     81:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda));
     82:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     83:       }
     84:       if (lambda <= steptol) {
     85:         PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
     86:         PetscFunctionReturn(PETSC_SUCCESS);
     87:       }
     88:       maxstep    = .95 * lambda; /* forbid the search from ever going back to the "failed" length that generates Nan or Inf */
     89:       lambda     = .5 * (lambda + lambda_old);
     90:       lambda_mid = .5 * (lambda + lambda_old);
     91:     }

     93:     delLambda = lambda - lambda_old;
     94:     /* compute f'() at the end points using second order one sided differencing */
     95:     delFnrm     = (3. * fnrm - 4. * fnrm_mid + 1. * fnrm_old) / delLambda;
     96:     delFnrm_old = (-3. * fnrm_old + 4. * fnrm_mid - 1. * fnrm) / delLambda;
     97:     /* compute f''() at the midpoint using centered differencing */
     98:     del2Fnrm = (delFnrm - delFnrm_old) / delLambda;

    100:     if (monitor) {
    101:       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    102:       if (!objective) {
    103:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: lambdas = [%g, %g, %g], fnorms = [%g, %g, %g]\n", (double)lambda, (double)lambda_mid, (double)lambda_old, (double)PetscSqrtReal(fnrm), (double)PetscSqrtReal(fnrm_mid), (double)PetscSqrtReal(fnrm_old)));
    104:       } else {
    105:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: lambdas = [%g, %g, %g], obj = [%g, %g, %g]\n", (double)lambda, (double)lambda_mid, (double)lambda_old, (double)fnrm, (double)fnrm_mid, (double)fnrm_old));
    106:       }
    107:       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    108:     }

    110:     /* compute the secant (Newton) update -- always go downhill */
    111:     if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
    112:     else if (del2Fnrm < 0.) lambda_update = lambda + delFnrm / del2Fnrm;
    113:     else break;

    115:     if (lambda_update < steptol) lambda_update = 0.5 * (lambda + lambda_old);

    117:     if (PetscIsInfOrNanReal(lambda_update)) break;

    119:     if (lambda_update > maxstep) break;

    121:     /* update the endpoints and the midpoint of the bracketed secant region */
    122:     lambda_old = lambda;
    123:     lambda     = lambda_update;
    124:     fnrm_old   = fnrm;
    125:     lambda_mid = 0.5 * (lambda + lambda_old);
    126:   }
    127:   /* construct the solution */
    128:   PetscCall(VecWAXPY(W, -lambda, Y, X));
    129:   if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));

    131:   /* postcheck */
    132:   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
    133:   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
    134:   if (changed_y) {
    135:     if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
    136:     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
    137:   }
    138:   PetscCall(VecCopy(W, X));
    139:   PetscCall((*linesearch->ops->snesfunc)(snes, X, F));

    141:   PetscCall(SNESLineSearchComputeNorms(linesearch));

    143:   if (monitor) {
    144:     PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &gnorm, NULL));
    145:     PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    146:     PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search terminated: lambda = %g, fnorm = %g\n", (double)lambda, (double)gnorm));
    147:     PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    148:   }
    149:   if (lambda <= steptol) PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
    150:   PetscFunctionReturn(PETSC_SUCCESS);
    151: }

    153: /*MC
    154:    SNESLINESEARCHL2 - Secant search in the L2 norm of the function or the objective function, if it is provided with `SNESSetObjective()`.

    156:    Attempts to solve $ \min_{\lambda} f(x + \lambda y) $ using the secant method with the initial bracketing of $ \lambda $ between [0,damping].
    157:    Differences of $f()$ are used to approximate the first and second derivative of $f()$ with respect to
    158:    $\lambda$, $f'()$ and $f''()$. The secant method is run for `maxit` iterations.

    160:    When an objective function is provided $f(w)$ is the objective function otherwise $f(w) = ||F(w)||^2$.
    161:    $x$ is the current step and $y$ is the search direction.

    163:    This has no checks on whether the secant method is actually converging.

    165:    Options Database Keys:
    166: +  -snes_linesearch_max_it <maxit>        - maximum number of iterations, default is 1
    167: .  -snes_linesearch_maxstep <length>      - the algorithm insures that a step length is never longer than this value
    168: .  -snes_linesearch_damping <damping>     - initial step is scaled back by this factor, default is 1.0
    169: -  -snes_linesearch_minlambda <minlambda> - minimum allowable lambda

    171:    Level: advanced

    173:    Developer Note:
    174:    A better name for this method might be `SNESLINESEARCHSECANT`, L2 is not descriptive

    176: .seealso: [](ch_snes), `SNESLINESEARCHBT`, `SNESLINESEARCHCP`, `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
    177: M*/
    178: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
    179: {
    180:   PetscFunctionBegin;
    181:   linesearch->ops->apply          = SNESLineSearchApply_L2;
    182:   linesearch->ops->destroy        = NULL;
    183:   linesearch->ops->setfromoptions = NULL;
    184:   linesearch->ops->reset          = NULL;
    185:   linesearch->ops->view           = NULL;
    186:   linesearch->ops->setup          = NULL;

    188:   linesearch->max_its = 1;
    189:   PetscFunctionReturn(PETSC_SUCCESS);
    190: }