Actual source code: morethuente.c

      1: /*$Id: morethuente.c 1.35 05/05/10 15:21:38-05:00 sarich@zorak.(none) $*/

      3: #include "morethuente.h"  /*I "tao_solver.h
      4: " I*/

      6: static int TaoApply_LineSearch(TAO_SOLVER,TaoVec*,TaoVec*,TaoVec*,TaoVec*,
      7:                                double*,double*,int*,void*);

      9: /*
     10:      The subroutine mcstep is taken from the work of Jorge Nocedal.
     11:      this is a variant of More' and Thuente's routine.

     13:      subroutine mcstep

     15:      the purpose of mcstep is to compute a safeguarded step for
     16:      a linesearch and to update an interval of uncertainty for
     17:      a minimizer of the function.

     19:      the parameter stx contains the step with the least function
     20:      value. the parameter stp contains the current step. it is
     21:      assumed that the derivative at stx is negative in the
     22:      direction of the step. if bracket is set true then a
     23:      minimizer has been bracketed in an interval of uncertainty
     24:      with endpoints stx and sty.

     26:      the subroutine statement is

     28:      subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
     29:                        stpmin,stpmax,info)

     31:      where

     33:        stx, fx, and dx are variables which specify the step,
     34:          the function, and the derivative at the best step obtained
     35:          so far. the derivative must be negative in the direction
     36:          of the step, that is, dx and stp-stx must have opposite
     37:          signs. on output these parameters are updated appropriately.

     39:        sty, fy, and dy are variables which specify the step,
     40:          the function, and the derivative at the other endpoint of
     41:          the interval of uncertainty. on output these parameters are
     42:          updated appropriately.

     44:        stp, fp, and dp are variables which specify the step,
     45:          the function, and the derivative at the current step.
     46:          if bracket is set true then on input stp must be
     47:          between stx and sty. on output stp is set to the new step.

     49:        bracket is a logical variable which specifies if a minimizer
     50:          has been bracketed. if the minimizer has not been bracketed
     51:          then on input bracket must be set false. if the minimizer
     52:          is bracketed then on output bracket is set true.

     54:        stpmin and stpmax are input variables which specify lower
     55:          and upper bounds for the step.

     57:        info is an integer output variable set as follows:
     58:          if info = 1,2,3,4,5, then the step has been computed
     59:          according to one of the five cases below. otherwise
     60:          info = 0, and this indicates improper input parameters.

     62:      subprograms called

     64:        fortran-supplied ... abs,max,min,sqrt

     66:      argonne national laboratory. minpack project. june 1983
     67:      jorge j. more', david j. thuente

     69:  */

     73: int TaoStep_LineSearch(TAO_SOLVER tao,double *stx,double *fx,double *dx,
     74:                        double *sty,double *fy,double *dy,double *stp,
     75:                        double *fp,double *dp)
     76: {
     77:   TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
     78:   double            gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
     79:   double            two = 2.0, zero = 0.0;
     80:   int               bound;

     82:   TaoFunctionBegin;
     83:   /* Check the input parameters for errors */
     84:   neP->infoc = 0;
     85:   if (neP->bracket && (*stp <= TaoMin(*stx,*sty) || (*stp >= TaoMax(*stx,*sty))))
     86:     SETERRQ(1,"bad stp in bracket");
     87:   if (*dx * (*stp-*stx) >= zero) SETERRQ(1,"dx * (stp-stx) >= 0");
     88:   if (neP->stepmax < neP->stepmin) SETERRQ(1,"stepmax > stepmin");

     90:   /* Determine if the derivatives have opposite sign */
     91:   sgnd = *dp * (*dx/TaoAbsDouble(*dx));

     93:   /* Case 1: a higher function value.
     94:      the minimum is bracketed. if the cubic step is closer
     95:      to stx than the quadratic step, the cubic step is taken,
     96:      else the average of the cubic and quadratic steps is taken.
     97:    */
     98:   if (*fp > *fx) {
     99:     neP->infoc = 1;
    100:     bound = 1;
    101:     theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
    102:     s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dx));
    103:     s = TaoMax(s,TaoAbsDouble(*dp));
    104:     gamma1 = s*sqrt(pow(theta/s,two) - (*dx/s)*(*dp/s));
    105:     if (*stp < *stx) gamma1 = -gamma1;
    106:     /* Can p be 0?  Check */
    107:     p = (gamma1 - *dx) + theta;
    108:     q = ((gamma1 - *dx) + gamma1) + *dp;
    109:     r = p/q;
    110:     stpc = *stx + r*(*stp - *stx);
    111:     stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);
    112:     if (TaoAbsDouble(stpc-*stx) < TaoAbsDouble(stpq-*stx)) {
    113:       stpf = stpc;
    114:     } else {
    115:       stpf = stpc + 0.5*(stpq - stpc);
    116:     }
    117:     neP->bracket = 1;
    118:   }
    119:   /*
    120:      Case 2: A lower function value and derivatives of
    121:      opposite sign. the minimum is bracketed. if the cubic
    122:      step is closer to stx than the quadratic (secant) step,
    123:      the cubic step is taken, else the quadratic step is taken.
    124:   */
    125:   else if (sgnd < zero) {
    126:     neP->infoc = 2;
    127:     bound = 0;
    128:     theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
    129:     s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dx));
    130:     s = TaoMax(s,TaoAbsDouble(*dp));
    131:     gamma1 = s*sqrt(pow(theta/s,two) - (*dx/s)*(*dp/s));
    132:     if (*stp > *stx) gamma1 = -gamma1;
    133:     p = (gamma1 - *dp) + theta;
    134:     q = ((gamma1 - *dp) + gamma1) + *dx;
    135:     r = p/q;
    136:     stpc = *stp + r*(*stx - *stp);
    137:     stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);
    138:     if (TaoAbsDouble(stpc-*stp) > TaoAbsDouble(stpq-*stp)) stpf = stpc;
    139:     else                                                       stpf = stpq;
    140:     neP->bracket = 1;
    141:   }

    143: /*   Case 3: A lower function value, derivatives of the
    144:      same sign, and the magnitude of the derivative decreases.
    145:      the cubic step is only used if the cubic tends to infinity
    146:      in the direction of the step or if the minimum of the cubic
    147:      is beyond stp. otherwise the cubic step is defined to be
    148:      either stepmin or stepmax. the quadratic (secant) step is also
    149:      computed and if the minimum is bracketed then the the step
    150:      closest to stx is taken, else the step farthest away is taken.
    151:  */

    153:   else if (TaoAbsDouble(*dp) < TaoAbsDouble(*dx)) {
    154:     neP->infoc = 3;
    155:     bound = 1;
    156:     theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
    157:     s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dx));
    158:     s = TaoMax(s,TaoAbsDouble(*dp));

    160:     /* The case gamma1 = 0 only arises if the cubic does not tend
    161:        to infinity in the direction of the step. */
    162:     gamma1 = s*sqrt(TaoMax(zero,pow(theta/s,two) - (*dx/s)*(*dp/s)));
    163:     if (*stp > *stx) gamma1 = -gamma1;
    164:     p = (gamma1 - *dp) + theta;
    165:     q = (gamma1 + (*dx - *dp)) + gamma1;
    166:     r = p/q;
    167:     if (r < zero && gamma1 != zero) stpc = *stp + r*(*stx - *stp);
    168:     else if (*stp > *stx)        stpc = neP->stepmax;
    169:     else                         stpc = neP->stepmin;
    170:     stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
    171:     if (neP->bracket) {
    172:       if (TaoAbsDouble(*stp-stpc) < TaoAbsDouble(*stp-stpq)) stpf = stpc;
    173:       else                                                       stpf = stpq;
    174:     }
    175:     else {
    176:       if (TaoAbsDouble(*stp-stpc) > TaoAbsDouble(*stp-stpq)) stpf = stpc;
    177:       else                                                       stpf = stpq;
    178:     }
    179:   }

    181: /*   Case 4: A lower function value, derivatives of the
    182:      same sign, and the magnitude of the derivative does
    183:      not decrease. if the minimum is not bracketed, the step
    184:      is either stpmin or stpmax, else the cubic step is taken.
    185:  */
    186:   else {
    187:     neP->infoc = 4;
    188:     bound = 0;
    189:     if (neP->bracket) {
    190:       theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
    191:       s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dy));
    192:       s = TaoMax(s,TaoAbsDouble(*dp));
    193:       gamma1 = s*sqrt(pow(theta/s,two) - (*dy/s)*(*dp/s));
    194:       if (*stp > *sty) gamma1 = -gamma1;
    195:       p = (gamma1 - *dp) + theta;
    196:       q = ((gamma1 - *dp) + gamma1) + *dy;
    197:       r = p/q;
    198:       stpc = *stp + r*(*sty - *stp);
    199:       stpf = stpc;
    200:     } else if (*stp > *stx) {
    201:       stpf = neP->stepmax;
    202:     } else {
    203:       stpf = neP->stepmin;
    204:     }
    205:   }

    207:   /* Update the interval of uncertainty.  This update does not
    208:      depend on the new step or the case analysis above. */

    210:   if (*fp > *fx) {
    211:     *sty = *stp;
    212:     *fy = *fp;
    213:     *dy = *dp;
    214:   } else {
    215:     if (sgnd < zero) {
    216:       *sty = *stx;
    217:       *fy = *fx;
    218:       *dy = *dx;
    219:     }
    220:     *stx = *stp;
    221:     *fx = *fp;
    222:     *dx = *dp;
    223:   }

    225:   /* Compute the new step and safeguard it */
    226:   stpf = TaoMin(neP->stepmax,stpf);
    227:   stpf = TaoMax(neP->stepmin,stpf);
    228:   *stp = stpf;
    229:   if (neP->bracket && bound) {
    230:     if (*sty > *stx) *stp = TaoMin(*stx+0.66*(*sty-*stx),*stp);
    231:     else             *stp = TaoMax(*stx+0.66*(*sty-*stx),*stp);
    232:   }
    233:   TaoFunctionReturn(0);
    234: }

    236: /* ---------------------------------------------------------- */
    239: int TaoDestroy_LineSearch(TAO_SOLVER tao, void *ctx)
    240: {
    241:   int  info;

    243:   TaoFunctionBegin;
    244:   info = TaoFree(tao->linectx);CHKERRQ(info);
    245:   TaoFunctionReturn(0);
    246: }
    247: /*------------------------------------------------------------*/
    250: int TaoSetOptions_LineSearch(TAO_SOLVER tao, void*linectx)
    251: {
    252:   TAO_LINESEARCH *ctx = (TAO_LINESEARCH *)linectx;
    253:   int            info;

    255:   TaoFunctionBegin;
    256:   info = TaoOptionsHead("More-Thuente line line search options for unconstrained minimization");CHKERRQ(info);
    257:   info = TaoOptionInt("-tao_ls_maxfev","max function evals in line search","",ctx->maxfev,&ctx->maxfev,0);CHKERRQ(info);
    258:   info = TaoOptionDouble("-tao_ls_ftol","tol for sufficient decrease","",ctx->ftol,&ctx->ftol,0);CHKERRQ(info);
    259:   info = TaoOptionDouble("-tao_ls_gtol","tol for curvature condition","",ctx->gtol,&ctx->gtol,0);CHKERRQ(info);
    260:   info = TaoOptionDouble("-tao_ls_rtol","relative tol for acceptable step","",ctx->rtol,&ctx->rtol,0);CHKERRQ(info);
    261:   info = TaoOptionDouble("-tao_ls_stepmin","lower bound for step","",ctx->stepmin,&ctx->stepmin,0);CHKERRQ(info);
    262:   info = TaoOptionDouble("-tao_ls_stepmax","upper bound for step","",ctx->stepmax,&ctx->stepmax,0);CHKERRQ(info);
    263:   info = TaoOptionsTail();CHKERRQ(info);
    264:   TaoFunctionReturn(0);
    265: }
    266: /*------------------------------------------------------------*/

    270: int TaoView_LineSearch(TAO_SOLVER tao, void*ctx)
    271: {
    272:   TAO_LINESEARCH *ls = (TAO_LINESEARCH *)ctx;
    273:   int            info;

    275:   TaoFunctionBegin;

    277:   info=TaoPrintInt(tao,"    More'-Thuente Line search: maxf=%d,",ls->maxfev);CHKERRQ(info);
    278:   info=TaoPrintDouble(tao," ftol=%g,",ls->ftol);CHKERRQ(info);
    279:   info=TaoPrintDouble(tao," rtol=%g,",ls->rtol);CHKERRQ(info);
    280:   info=TaoPrintDouble(tao," gtol=%g\n",ls->gtol);CHKERRQ(info);

    282:   TaoFunctionReturn(0);
    283: }

    286: /* ---------------------------------------------------------- */
    289: /* @
    290:    TaoApply_LineSearch - This routine performs a line search algorithm.

    292:    Input Parameters:
    293: +  tao - TAO_SOLVER context
    294: .  X - current iterate (on output X contains new iterate, X + step*S)
    295: .  S - search direction
    296: .  f - objective function evaluated at X
    297: .  G - gradient evaluated at X
    298: .  W - work vector
    299: -  step - initial estimate of step length

    301:    Output parameters:
    302: +  f - objective function evaluated at new iterate, X + step*S
    303: .  G - gradient evaluated at new iterate, X + step*S
    304: .  X - new iterate
    305: .  gnorm - 2-norm of G
    306: -  step - final step length

    309:    Info is set to one of:
    310: +   0 - the line search succeeds; the sufficient decrease
    311:    condition and the directional derivative condition hold

    313:    negative number if an input parameter is invalid
    314: .   -1 -  step < 0
    315: .   -2 -  ftol < 0
    316: .   -3 -  rtol < 0
    317: .   -4 -  gtol < 0
    318: .   -5 -  stepmin < 0
    319: .   -6 -  stepmax < stepmin
    320: -   -7 -  maxfev < 0

    322:    positive number > 1 if the line search otherwise terminates
    323: +    2 -  Relative width of the interval of uncertainty is
    324:          at most rtol.
    325: .    3 -  Maximum number of function evaluations (maxfev) has
    326:          been reached.
    327: .    4 -  Step is at the lower bound, stepmin.
    328: .    5 -  Step is at the upper bound, stepmax.
    329: .    6 -  Rounding errors may prevent further progress.
    330:          There may not be a step that satisfies the
    331:          sufficient decrease and curvature conditions.
    332:          Tolerances may be too small.
    333: -    7 -  Search direction is not a descent direction.

    335:    Notes:
    336:    This algorithm is taken from More' and Thuente, "Line search algorithms
    337:    with guaranteed sufficient decrease", Argonne National Laboratory,
    338:    Technical Report MCS-P330-1092.

    340:    Notes:
    341:    This routine is used within the following TAO unconstrained minimization
    342:    solvers: Newton linesearch (tao_nls), limited memory variable metric
    343:    (tao_lmvm), and conjugate gradient (tao_cg).

    345:    Level: advanced

    347: .keywords: TAO_SOLVER, linesearch
    348: @ */
    349: static int TaoApply_LineSearch(TAO_SOLVER tao,TaoVec* X,TaoVec* G,TaoVec* S,TaoVec* W,double *f,
    350:                         double *step,int *info2,void*ctx)
    351: {
    352:   TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
    353:   double    zero = 0.0, two = 2.0, p5 = 0.5, p66 = 0.66, xtrapf = 4.0;
    354:   double    finit, width, width1, dginit,fm, fxm, fym, dgm, dgxm, dgym;
    355:   double    dgx, dgy, dg, fx, fy, stx, sty, dgtest, ftest1=0.0;
    356:   int       info, i, stage1;

    358: #if defined(PETSC_USE_COMPLEX)
    359:   PetscScalar    cdginit, cdg, cstep = 0.0;
    360: #endif

    362:   TaoFunctionBegin;
    363:   /* neP->stepmin - lower bound for step */
    364:   /* neP->stepmax - upper bound for step */
    365:   /* neP->rtol           - relative tolerance for an acceptable step */
    366:   /* neP->ftol           - tolerance for sufficient decrease condition */
    367:   /* neP->gtol           - tolerance for curvature condition */
    368:   /* neP->nfev           - number of function evaluations */
    369:   /* neP->maxfev  - maximum number of function evaluations */

    371:   /* Check input parameters for errors */
    372:   if (*step < zero) {
    373:     info = PetscLogInfo((tao,"TaoApply_LineSearch:Line search error: step (%g) < 0\n",*step)); CHKERRQ(info);
    374:     *info2 = -1; TaoFunctionReturn(0);
    375:   } else if (neP->ftol < zero) {
    376:     info = PetscLogInfo((tao,"TaoApply_LineSearch:Line search error: ftol (%g) < 0\n",neP->ftol)); CHKERRQ(info);
    377:     *info2 = -2; TaoFunctionReturn(0);
    378:   } else if (neP->rtol < zero) {
    379:     info = PetscLogInfo((tao,"TaoApply_LineSearch:Line search error: rtol (%g) < 0\n",neP->rtol)); CHKERRQ(info);
    380:     *info2 = -3; TaoFunctionReturn(0);
    381:   } else if (neP->gtol < zero) {
    382:     info = PetscLogInfo((tao,"TaoApply_LineSearch:Line search error: gtol (%g) < 0\n",neP->gtol)); CHKERRQ(info);
    383:     *info2 = -4; TaoFunctionReturn(0);
    384:   } else if (neP->stepmin < zero) {
    385:     info = PetscLogInfo((tao,"TaoApply_LineSearch:Line search error: stepmin (%g) < 0\n",neP->stepmin)); CHKERRQ(info);
    386:     *info2 = -5; TaoFunctionReturn(0);
    387:   } else if (neP->stepmax < neP->stepmin) {
    388:     info = PetscLogInfo((tao,"TaoApply_LineSearch:Line search error: stepmax (%g) < stepmin (%g)\n", neP->stepmax,neP->stepmin)); CHKERRQ(info);
    389:     *info2 = -6; TaoFunctionReturn(0);
    390:   } else if (neP->maxfev < zero) {
    391:     info = PetscLogInfo((tao,"TaoApply_LineSearch:Line search error: maxfev (%d) < 0\n",neP->maxfev)); CHKERRQ(info);
    392:     *info2 = -7; TaoFunctionReturn(0);
    393:   }

    395:   /* Check that search direction is a descent direction */

    397: #if defined(PETSC_USE_COMPLEX)
    398:   info = G->Dot(S,&cdginit);CHKERRQ(info); dginit = TaoReal(cdginit);
    399: #else
    400:   info = G->Dot(S,&dginit);CHKERRQ(info);
    401: #endif

    403:   if (dginit >= zero) {
    404:     info = PetscLogInfo((tao,"TaoApply_LineSearch:Search direction not a descent direction\n")); CHKERRQ(info);
    405:     *info2 = 7; TaoFunctionReturn(0);
    406:   }

    408:   /* Initialization */
    409:   neP->bracket = 0;
    410:   *info2          = 0;
    411:   stage1  = 1;
    412:   finit   = *f;
    413:   dgtest  = neP->ftol * dginit;
    414:   width   = neP->stepmax - neP->stepmin;
    415:   width1  = width * two;
    416:   info = W->CopyFrom(X);CHKERRQ(info);
    417:   /* Variable dictionary:
    418:      stx, fx, dgx - the step, function, and derivative at the best step
    419:      sty, fy, dgy - the step, function, and derivative at the other endpoint
    420:                    of the interval of uncertainty
    421:      step, f, dg - the step, function, and derivative at the current step */

    423:   stx = zero;
    424:   fx  = finit;
    425:   dgx = dginit;
    426:   sty = zero;
    427:   fy  = finit;
    428:   dgy = dginit;
    429:
    430:   neP->nfev = 0;
    431:   for (i=0; i< neP->maxfev; i++) {
    432:     /* Set min and max steps to correspond to the interval of uncertainty */
    433:     if (neP->bracket) {
    434:       neP->stepmin = TaoMin(stx,sty);
    435:       neP->stepmax = TaoMax(stx,sty);
    436:     } else {
    437:       neP->stepmin = stx;
    438:       neP->stepmax = *step + xtrapf * (*step - stx);
    439:     }

    441:     /* Force the step to be within the bounds */
    442:     *step = TaoMax(*step,neP->stepmin);
    443:     *step = TaoMin(*step,neP->stepmax);

    445:     /* If an unusual termination is to occur, then let step be the lowest
    446:        point obtained thus far */
    447:     if (((neP->bracket) && (*step <= neP->stepmin || *step >= neP->stepmax)) ||
    448:         ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol * neP->stepmax)) ||
    449:         (neP->nfev >= neP->maxfev - 1) || (neP->infoc == 0))
    450:       *step = stx;
    451:
    452: #if defined(PETSC_USE_COMPLEX)
    453:     cstep = *step;
    454:     info = X->Waxpby(cstep,S,1.0,W);CHKERRQ(info);
    455: #else
    456:     info = X->Waxpby(*step,S,1.0,W);CHKERRQ(info);         /* X = W + step*S */
    457: #endif
    458:     info = TaoComputeMeritFunctionGradient(tao,X,f,G);CHKERRQ(info);
    459:     neP->nfev++;
    460: #if defined(PETSC_USE_COMPLEX)
    461:     info = G->Dot(S,&cdg);CHKERRQ(info); dg = TaoReal(cdg);
    462: #else
    463:     info = G->Dot(S,&dg);CHKERRQ(info);                /* dg = G^T S */
    464: #endif
    465:     ftest1 = finit + *step * dgtest;

    467:     /* Convergence testing */
    468:     if (((neP->bracket) && (*step <= neP->stepmin||*step >= neP->stepmax)) || (!neP->infoc)) {
    469:       info = PetscLogInfo((tao,"TaoApply_LineSearch:Rounding errors may prevent further progress.  May not be a step satisfying\n")); CHKERRQ(info);
    470:       info = PetscLogInfo((tao,"TaoApply_LineSearch:sufficient decrease and curvature conditions. Tolerances may be too small.\n")); CHKERRQ(info);
    471:       *info2 = 6; break;
    472:     }
    473:     if ((*step == neP->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
    474:       info = PetscLogInfo((tao,"TaoApply_LineSearch:Step is at the upper bound, stepmax (%g)\n",neP->stepmax)); CHKERRQ(info);
    475:       *info2 = 5; break;
    476:     }
    477:     if ((*step == neP->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
    478:       info = PetscLogInfo((tao,"TaoApply_LineSearch:Step is at the lower bound, stepmin (%g)\n",neP->stepmin)); CHKERRQ(info);
    479:       *info2 = 4; break;
    480:     }
    481:     if (neP->nfev >= neP->maxfev) {
    482:       info = PetscLogInfo((tao,"TaoApply_LineSearch:Number of line search function evals (%d) > maximum (%d)\n",neP->nfev,neP->maxfev)); CHKERRQ(info);
    483:       *info2 = 3; break;
    484:     }
    485:     if ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol*neP->stepmax)){
    486:       info = PetscLogInfo((tao,"TaoApply_LineSearch:Relative width of interval of uncertainty is at most rtol (%g)\n",neP->rtol)); CHKERRQ(info);
    487:       *info2 = 2; break;
    488:     }
    489:     if ((*f <= ftest1) && (TaoAbsDouble(dg) <= neP->gtol*(-dginit))) {
    490:       info = PetscLogInfo((tao,"TaoApply_LineSearch:Line search success: Sufficient decrease and directional deriv conditions hold\n")); CHKERRQ(info);
    491:       *info2 = 0; break;
    492:     } else {
    493:       if (*f > ftest1){
    494:         info = PetscLogInfo((tao,"TaoApply_LineSearch: Sufficient decrease only if %g < %g\n",*f,ftest1)); CHKERRQ(info);
    495:       }
    496:       if (TaoAbsDouble(dg) > neP->gtol*(-dginit)){
    497:         info = PetscLogInfo((tao,"TaoApply_LineSearch: Sufficient directional deriv only if %g < %g * %g\n",TaoAbsDouble(dg),neP->gtol,(-dginit))); CHKERRQ(info);
    498:       }
    499:     }

    501:     /* In the first stage, we seek a step for which the modified function
    502:         has a nonpositive value and nonnegative derivative */
    503:     if ((stage1) && (*f <= ftest1) && (dg >= dginit * TaoMin(neP->ftol, neP->gtol)))
    504:       stage1 = 0;

    506:     /* A modified function is used to predict the step only if we
    507:        have not obtained a step for which the modified function has a
    508:        nonpositive function value and nonnegative derivative, and if a
    509:        lower function value has been obtained but the decrease is not
    510:        sufficient */

    512:     if ((stage1) && (*f <= fx) && (*f > ftest1)) {
    513:       fm   = *f - *step * dgtest;        /* Define modified function */
    514:       fxm  = fx - stx * dgtest;                /* and derivatives */
    515:       fym  = fy - sty * dgtest;
    516:       dgm  = dg - dgtest;
    517:       dgxm = dgx - dgtest;
    518:       dgym = dgy - dgtest;

    520:       /* Update the interval of uncertainty and compute the new step */
    521:       info = TaoStep_LineSearch(tao,&stx,&fxm,&dgxm,&sty,&fym,&dgym,step,&fm,&dgm);CHKERRQ(info);

    523:       fx  = fxm + stx * dgtest;        /* Reset the function and */
    524:       fy  = fym + sty * dgtest;        /* gradient values */
    525:       dgx = dgxm + dgtest;
    526:       dgy = dgym + dgtest;
    527:     } else {
    528:       /* Update the interval of uncertainty and compute the new step */
    529:       info = TaoStep_LineSearch(tao,&stx,&fx,&dgx,&sty,&fy,&dgy,step,f,&dg);CHKERRQ(info);
    530:     }

    532:    /* Force a sufficient decrease in the interval of uncertainty */
    533:    if (neP->bracket) {
    534:      if (TaoAbsDouble(sty - stx) >= p66 * width1) *step = stx + p5*(sty - stx);
    535:        width1 = width;
    536:        width = TaoAbsDouble(sty - stx);
    537:      }
    538:    }

    540:   /* Finish computations */
    541:   info = PetscLogInfo((tao,"TaoApply_LineSearch:%d function evals in line search, step = %10.4f\n",neP->nfev,*step)); CHKERRQ(info);
    542:   /*
    543:   info = G->Norm2(dginitt);CHKERRQ(info);
    544:   */
    545:   TaoFunctionReturn(0);
    546: }

    548: /* ---------------------------------------------------------- */

    552: /*@C
    553:    TaoCreateMoreThuenteLineSearch - Create a line search

    555:    Input Parameters:
    556: +  tao - TAO_SOLVER context
    557: .  fftol - the sufficient descent parameter , greater than 0.
    558: -  ggtol - the curvature tolerance, greater than 0, less than 1.

    561:    Note:
    562:    If either fftol or ggtol is 0, default parameters will be used.

    564:    Note:
    565:    This algorithm is taken from More' and Thuente, "Line search algorithms
    566:    with guaranteed sufficient decrease", Argonne National Laboratory,
    567:    Technical Report MCS-P330-1092.

    569:    Note:
    570:    This line search enforces the strong Wolfe conditions for unconstrained
    571:    optimization.  This routine is used within the following TAO unconstrained
    572:    minimization solvers: Newton linesearch (tao_nls), limited memory variable
    573:    metric (tao_lmvm), and nonlinear conjugate gradient methods.

    575:    Level: developer

    577: .keywords: TAO_SOLVER, linesearch
    578: @*/
    579: int TaoCreateMoreThuenteLineSearch(TAO_SOLVER tao, double fftol, double ggtol)
    580: {
    581:   int info;
    582:   TAO_LINESEARCH *neP;

    584:   TaoFunctionBegin;

    586:   info = TaoNew(TAO_LINESEARCH,&neP); CHKERRQ(info);
    587:   info = PetscLogObjectMemory(tao,sizeof(TAO_LINESEARCH)); CHKERRQ(info);

    589:   if (fftol<=0) fftol=0.001;
    590:   if (ggtol<=0) ggtol=0.9;

    592:   neP->ftol                  = fftol;
    593:   neP->rtol                  = 1.0e-10;
    594:   neP->gtol                  = ggtol;
    595:   neP->stepmin                  = 1.0e-20;
    596:   neP->stepmax                  = 1.0e+20;
    597:   neP->nfev                  = 0;
    598:   neP->bracket                  = 0;
    599:   neP->infoc              = 1;
    600:   neP->maxfev                  = 30;

    602:   info = TaoSetLineSearch(tao,0,
    603:                           TaoSetOptions_LineSearch,
    604:                           TaoApply_LineSearch,
    605:                           TaoView_LineSearch,
    606:                           TaoDestroy_LineSearch,
    607:                           (void *) neP);CHKERRQ(info);

    609:   TaoFunctionReturn(0);
    610: }

    615: /*

    618: */