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;
      8:   Vec            X;
      9:   Vec            F;
     10:   Vec            Y;
     11:   Vec            W;
     12:   SNES           snes;
     13:   PetscReal      gnorm;
     14:   PetscReal      ynorm;
     15:   PetscReal      xnorm;
     16:   PetscReal      steptol, maxstep, rtol, atol, ltol;
     17:   PetscViewer    monitor;
     18:   PetscReal      lambda, lambda_old, lambda_mid, lambda_update, delLambda;
     19:   PetscReal      fnrm, fnrm_old, fnrm_mid;
     20:   PetscReal      delFnrm, delFnrm_old, del2Fnrm;
     21:   PetscInt       i, max_its;
     22:   PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);

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

     33:   SNESGetObjective(snes,&objective,NULL);

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

     45:   for (i = 0; i < max_its; i++) {

     47:     while (PETSC_TRUE) {
     48:       VecCopy(X, W);
     49:       VecAXPY(W, -lambda_mid, Y);
     50:       if (linesearch->ops->viproject) {
     51:         (*linesearch->ops->viproject)(snes, W);
     52:       }
     53:       if (!objective) {
     54:         /* compute the norm at the midpoint */
     55:         (*linesearch->ops->snesfunc)(snes, W, F);
     56:         if (linesearch->ops->vinorm) {
     57:           fnrm_mid = gnorm;
     58:           (*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid);
     59:         } else {
     60:           VecNorm(F,NORM_2,&fnrm_mid);
     61:         }

     63:         /* compute the norm at the new endpoit */
     64:         VecCopy(X, W);
     65:         VecAXPY(W, -lambda, Y);
     66:         if (linesearch->ops->viproject) {
     67:           (*linesearch->ops->viproject)(snes, W);
     68:         }
     69:         (*linesearch->ops->snesfunc)(snes, W, F);
     70:         if (linesearch->ops->vinorm) {
     71:           fnrm = gnorm;
     72:           (*linesearch->ops->vinorm)(snes, F, W, &fnrm);
     73:         } else {
     74:           VecNorm(F,NORM_2,&fnrm);
     75:         }
     76:         fnrm_mid = fnrm_mid*fnrm_mid;
     77:         fnrm = fnrm*fnrm;
     78:       } else {
     79:         /* compute the objective at the midpoint */
     80:         VecCopy(X, W);
     81:         VecAXPY(W, -lambda_mid, Y);
     82:         SNESComputeObjective(snes,W,&fnrm_mid);

     84:         /* compute the objective at the new endpoint */
     85:         VecCopy(X, W);
     86:         VecAXPY(W, -lambda, Y);
     87:         SNESComputeObjective(snes,W,&fnrm);
     88:       }
     89:       if (!PetscIsInfOrNanReal(fnrm)) break;
     90:       if (monitor) {
     91:         PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
     92:         PetscViewerASCIIPrintf(monitor,"    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda);
     93:         PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
     94:       }
     95:       if (lambda <= steptol) {
     96:         SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
     97:         return(0);
     98:       }
     99:       maxstep = .95*lambda; /* forbid the search from ever going back to the "failed" length that generates Nan or Inf */
    100:       lambda  = .5*(lambda + lambda_old);
    101:       lambda_mid = .5*(lambda + lambda_old);
    102:     }

    104:     delLambda   = lambda - lambda_old;
    105:     /* compute f'() at the end points using second order one sided differencing */
    106:     delFnrm     = (3.*fnrm - 4.*fnrm_mid + 1.*fnrm_old) / delLambda;
    107:     delFnrm_old = (-3.*fnrm_old + 4.*fnrm_mid -1.*fnrm) / delLambda;
    108:     /* compute f''() at the midpoint using centered differencing */
    109:     del2Fnrm    = (delFnrm - delFnrm_old) / delLambda;

    111:     if (monitor) {
    112:       PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
    113:       if (!objective) {
    114:         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));
    115:       } else {
    116:         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);
    117:       }
    118:       PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
    119:     }

    121:     /* compute the secant (Newton) update -- always go downhill */
    122:     if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
    123:     else if (del2Fnrm < 0.) lambda_update = lambda + delFnrm / del2Fnrm;
    124:     else break;

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

    128:     if (PetscIsInfOrNanReal(lambda_update)) break;

    130:     if (lambda_update > maxstep) break;

    132:     /* update the endpoints and the midpoint of the bracketed secant region */
    133:     lambda_old = lambda;
    134:     lambda     = lambda_update;
    135:     fnrm_old   = fnrm;
    136:     lambda_mid = 0.5*(lambda + lambda_old);
    137:   }
    138:   /* construct the solution */
    139:   VecCopy(X, W);
    140:   VecAXPY(W, -lambda, Y);
    141:   if (linesearch->ops->viproject) {
    142:     (*linesearch->ops->viproject)(snes, W);
    143:   }

    145:   /* postcheck */
    146:   SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
    147:   if (changed_y) {
    148:     VecAXPY(X, -lambda, Y);
    149:     if (linesearch->ops->viproject) {
    150:       (*linesearch->ops->viproject)(snes, X);
    151:     }
    152:   } else {
    153:     VecCopy(W, X);
    154:   }
    155:   (*linesearch->ops->snesfunc)(snes,X,F);

    157:   SNESLineSearchSetLambda(linesearch, lambda);
    158:   SNESLineSearchComputeNorms(linesearch);
    159:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);

    161:   if (monitor) {
    162:     PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
    163:     PetscViewerASCIIPrintf(monitor,"    Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
    164:     PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
    165:   }
    166:   if (lambda <= steptol) {
    167:     SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
    168:   }
    169:   return(0);
    170: }

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

    175:    Attempts to solve min_lambda f(x + lambda y) using the secant method with the initial bracketing of lambda between [0,damping]. Differences of f()
    176:    are used to approximate the first and second derivative of f() with respect to lambda, f'() and f''(). The secant method is run for maxit iterations.

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

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

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

    188:    Level: advanced

    190:    Developer Notes:
    191:     A better name for this method might be SNESLINESEARCHSECANT, L2 is not descriptive

    193: .seealso: SNESLINESEARCHBT, SNESLINESEARCHCP, SNESLineSearch, SNESLineSearchCreate(), SNESLineSearchSetType()
    194: M*/
    195: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
    196: {
    198:   linesearch->ops->apply          = SNESLineSearchApply_L2;
    199:   linesearch->ops->destroy        = NULL;
    200:   linesearch->ops->setfromoptions = NULL;
    201:   linesearch->ops->reset          = NULL;
    202:   linesearch->ops->view           = NULL;
    203:   linesearch->ops->setup          = NULL;

    205:   linesearch->max_its = 1;
    206:   return(0);
    207: }