Actual source code: linesearchbisection.c

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

      4: static PetscErrorCode SNESLineSearchApply_Bisection(SNESLineSearch linesearch)
      5: {
      6:   PetscBool   changed_y, changed_w;
      7:   Vec         X, F, Y, W, G;
      8:   SNES        snes;
      9:   PetscReal   ynorm;
     10:   PetscReal   lambda_left, lambda, lambda_right, lambda_old, fnorm;
     11:   PetscScalar fty_left, fty, fty_initial;
     12:   PetscViewer monitor;
     13:   PetscReal   rtol, atol, ltol;
     14:   PetscInt    it, max_it;

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

     24:   /* pre-check */
     25:   PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));

     27:   /* compute ynorm to normalize search direction */
     28:   PetscCall(VecNorm(Y, NORM_2, &ynorm));

     30:   /* initialize interval for bisection */
     31:   lambda_left  = 0.0;
     32:   lambda_right = lambda;

     34:   /* compute fty at left end of interval */
     35:   if (linesearch->ops->vidirderiv) {
     36:     PetscCall((*linesearch->ops->vidirderiv)(snes, F, X, Y, &fty_left));
     37:   } else {
     38:     PetscCall(VecDot(F, Y, &fty_left));
     39:   }
     40:   fty_initial = fty_left;

     42:   /* compute fty at right end of interval (initial lambda) */
     43:   PetscCall(VecWAXPY(W, -lambda, Y, X));
     44:   if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
     45:   PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
     46:   if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
     47:     PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations during line search!\n"));
     48:     snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
     49:     PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
     50:     PetscFunctionReturn(PETSC_SUCCESS);
     51:   }
     52:   if (linesearch->ops->vidirderiv) {
     53:     PetscCall((*linesearch->ops->vidirderiv)(snes, G, W, Y, &fty));
     54:   } else {
     55:     PetscCall(VecDot(G, Y, &fty));
     56:   }
     57:   /* check whether sign changes in interval */
     58:   if (!PetscIsInfOrNanScalar(fty) && (PetscRealPart(fty_left * fty) > 0.0)) {
     59:     /* no change of sign: accept full step */
     60:     if (monitor) {
     61:       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     62:       PetscCall(PetscViewerASCIIPrintf(monitor, "      Line search: sign of fty does not change in step interval, accepting full step\n"));
     63:       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     64:     }
     65:   } else {
     66:     /* change of sign: iteratively bisect interval */
     67:     lambda_old = 0.0;
     68:     it         = 0;

     70:     while (PETSC_TRUE) {
     71:       /* check for infinity or NaN */
     72:       if (PetscIsInfOrNanScalar(fty)) {
     73:         if (monitor) {
     74:           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     75:           PetscCall(PetscViewerASCIIPrintf(monitor, "      Line search fty is infinity or NaN!\n"));
     76:           PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     77:         }
     78:         PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_CONV_FAILED, "infinity or NaN in function evaluation");
     79:         if (snes->functiondomainerror) {
     80:           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
     81:           snes->functiondomainerror = PETSC_FALSE;
     82:         } else PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
     83:         PetscFunctionReturn(PETSC_SUCCESS);
     84:         break;
     85:       }

     87:       /* check absolute tolerance */
     88:       if (PetscAbsScalar(fty) <= atol * ynorm) {
     89:         if (monitor) {
     90:           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     91:           PetscCall(PetscViewerASCIIPrintf(monitor, "      Line search: abs(fty)/||y|| = %g <= atol = %g\n", (double)(PetscAbsScalar(fty) / ynorm), (double)atol));
     92:           PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     93:         }
     94:         break;
     95:       }

     97:       /* check relative tolerance */
     98:       if (PetscAbsScalar(fty) / PetscAbsScalar(fty_initial) <= rtol) {
     99:         if (monitor) {
    100:           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    101:           PetscCall(PetscViewerASCIIPrintf(monitor, "      Line search: abs(fty/fty_initial) = %g <= rtol  = %g\n", (double)(PetscAbsScalar(fty) / PetscAbsScalar(fty_initial)), (double)rtol));
    102:           PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    103:         }
    104:         break;
    105:       }

    107:       /* check maximum number of iterations */
    108:       if (it > max_it) {
    109:         if (monitor) {
    110:           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    111:           PetscCall(PetscViewerASCIIPrintf(monitor, "      Line search: maximum iterations reached\n"));
    112:           PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    113:         }
    114:         PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
    115:         break;
    116:       }

    118:       /* check change of lambda tolerance */
    119:       if (PetscAbsReal(lambda - lambda_old) < ltol) {
    120:         if (monitor) {
    121:           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    122:           PetscCall(PetscViewerASCIIPrintf(monitor, "      Line search: abs(dlambda) = %g < ltol = %g\n", (double)PetscAbsReal(lambda - lambda_old), (double)ltol));
    123:           PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    124:         }
    125:         break;
    126:       }

    128:       /* determine direction of bisection (not necessary for 0th iteration) */
    129:       if (it > 0) {
    130:         if (PetscRealPart(fty * fty_left) <= 0.0) {
    131:           lambda_right = lambda;
    132:         } else {
    133:           lambda_left = lambda;
    134:           /* also update fty_left for direction check in next iteration */
    135:           fty_left = fty;
    136:         }
    137:       }

    139:       /* bisect interval */
    140:       lambda_old = lambda;
    141:       lambda     = 0.5 * (lambda_left + lambda_right);

    143:       /* compute fty at new lambda */
    144:       PetscCall(VecWAXPY(W, -lambda, Y, X));
    145:       if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
    146:       PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
    147:       if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
    148:         PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations during line search!\n"));
    149:         snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
    150:         PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
    151:         PetscFunctionReturn(PETSC_SUCCESS);
    152:       }
    153:       if (linesearch->ops->vidirderiv) {
    154:         PetscCall((*linesearch->ops->vidirderiv)(snes, G, W, Y, &fty));
    155:       } else {
    156:         PetscCall(VecDot(G, Y, &fty));
    157:       }

    159:       /* print iteration information */
    160:       if (monitor) {
    161:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    162:         PetscCall(PetscViewerASCIIPrintf(monitor, "      %3" PetscInt_FMT " Line search: fty/||y|| = %g, lambda = %g\n", it, (double)(PetscRealPart(fty) / ynorm), (double)lambda));
    163:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    164:       }

    166:       /* count up iteration */
    167:       it++;
    168:     }
    169:   }

    171:   /* post-check */
    172:   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
    173:   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
    174:   if (changed_y) {
    175:     if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
    176:     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
    177:   }

    179:   /* update solution*/
    180:   PetscCall(VecCopy(W, X));
    181:   PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
    182:   PetscCall(SNESLineSearchComputeNorms(linesearch));
    183:   PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &fnorm, NULL));
    184:   SNESLineSearchCheckFunctionDomainError(snes, linesearch, fnorm);
    185:   PetscFunctionReturn(PETSC_SUCCESS);
    186: }

    188: /*MC
    189:    SNESLINESEARCHBISECTION - Bisection line search.
    190:    Similar to the critical point line search, `SNESLINESEARCHCP`, the bisection line search assumes that there exists some $G(x)$ for which the `SNESFunctionFn` $F(x) = grad G(x)$.
    191:    This line search seeks to find the root of the directional derivative, that is $F(x_k - \lambda Y_k) \cdot Y_k / ||Y_k|| = 0$, along the search direction $Y_k$ through bisection.

    193:    Options Database Keys:
    194: +  -snes_linesearch_max_it 50   - maximum number of bisection iterations for the line search
    195: .  -snes_linesearch_damping 1.0 - initial `lambda` on entry to the line search
    196: .  -snes_linesearch_rtol 1e\-8  - relative tolerance for the directional derivative
    197: .  -snes_linesearch_atol 1e\-6  - absolute tolerance for the directional derivative
    198: -  -snes_linesearch_ltol 1e\-6  - minimum absolute change in `lambda` allowed (this is an alternative to setting a maximum number of iterations)

    200:    Level: intermediate

    202:    Notes:
    203:    `lambda` is the scaling of the search direction (vector) that is computed by this algorithm.
    204:    If there is no change of sign in the directional derivative from $\lambda=0$ to the initial `lambda` (the damping), then the initial `lambda` will be used.
    205:    Hence, this line search will always give a `lambda` in the interval $[0, damping]$.
    206:    This method does NOT use the objective function if it is provided with `SNESSetObjective()`.

    208: .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHCP`
    209: M*/
    210: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Bisection(SNESLineSearch linesearch)
    211: {
    212:   PetscFunctionBegin;
    213:   linesearch->ops->apply          = SNESLineSearchApply_Bisection;
    214:   linesearch->ops->destroy        = NULL;
    215:   linesearch->ops->setfromoptions = NULL;
    216:   linesearch->ops->reset          = NULL;
    217:   linesearch->ops->view           = NULL;
    218:   linesearch->ops->setup          = NULL;

    220:   /* set default option values */
    221:   linesearch->max_it = 50;
    222:   linesearch->rtol   = 1e-8;
    223:   linesearch->atol   = 1e-6;
    224:   linesearch->ltol   = 1e-6;
    225:   PetscFunctionReturn(PETSC_SUCCESS);
    226: }