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: */