Actual source code: morethuente.c

  1: #include "petscvec.h"
  2: #include "taosolver.h"
  3: #include "private/taolinesearch_impl.h"
  4: #include "morethuente.h"

  6: /*
  7:    This algorithm is taken from More' and Thuente, "Line search algorithms
  8:    with guaranteed sufficient decrease", Argonne National Laboratory, 
  9:    Technical Report MCS-P330-1092.
 10: */

 12: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,
 13:                                  PetscReal *stx, PetscReal *fx, PetscReal *dx,
 14:                                  PetscReal *sty, PetscReal *fy, PetscReal *dy,
 15:                                  PetscReal *stp, PetscReal *fp, PetscReal *dp);

 19: static PetscErrorCode TaoLineSearchDestroy_MT(TaoLineSearch ls)
 20: {
 22:   TAOLINESEARCH_MT_CTX *mt;
 25:   mt = (TAOLINESEARCH_MT_CTX*)(ls->data);
 26:   if (mt->x) {
 27:     PetscObjectDereference((PetscObject)mt->x);  
 28:   }
 29:   if (mt->work) {
 30:     VecDestroy(&mt->work); 
 31:   }
 32:   PetscFree(ls->data);
 33:   ls->data = PETSC_NULL;
 34:   return(0);
 35: }

 39: static PetscErrorCode TaoLineSearchSetFromOptions_MT(TaoLineSearch ls)
 40: {
 42:   
 44:   
 45:   return(0);
 46: }


 51: static PetscErrorCode TaoLineSearchView_MT(TaoLineSearch ls, PetscViewer pv)
 52: {
 54:     PetscBool isascii;
 56:     PetscTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii); 
 57:     if (isascii) {
 58:         PetscViewerASCIIPrintf(pv,"  maxf=%D, ftol=%G, gtol=%G\n",ls->max_funcs, ls->rtol, ls->ftol); 
 59:     } else {
 60:       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for MoreThuente TaoLineSearch",((PetscObject)pv)->type_name);
 61:     }
 62:     return(0);
 63: }

 67: /* @ TaoApply_LineSearch - This routine takes step length of 1.0.

 69:    Input Parameters:
 70: +  tao - TaoSolver context
 71: .  X - current iterate (on output X contains new iterate, X + step*S)
 72: .  f - objective function evaluated at X
 73: .  G - gradient evaluated at X
 74: -  D - search direction


 77:    Info is set to 0.

 79: @ */

 81: static PetscErrorCode TaoLineSearchApply_MT(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s)
 82: {
 84:     TAOLINESEARCH_MT_CTX *mt;
 85:     
 86:     PetscReal    xtrapf = 4.0;
 87:     PetscReal   finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
 88:     PetscReal    dgx, dgy, dg, dg2, fx, fy, stx, sty, dgtest;
 89:     PetscReal ftest1=0.0, ftest2=0.0;
 90:     PetscInt  i, stage1,n1,n2,nn1,nn2;
 91:     PetscReal bstepmin1, bstepmin2, bstepmax;
 92:     PetscBool g_computed=PETSC_FALSE;; /* to prevent extra gradient computation */


101:     /* comm,type,size checks are done in interface TaoLineSearchApply */
102:     mt = (TAOLINESEARCH_MT_CTX*)(ls->data);

104:     ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;

106:     /* Check work vector */
107:     if (!mt->work) {
108:       VecDuplicate(x,&mt->work); 
109:       mt->x = x;
110:       PetscObjectReference((PetscObject)mt->x); 
111:     }
112:     /* If x has changed, then recreate work */
113:     else if (x != mt->x) { 
114:       VecDestroy(&mt->work); 
115:       VecDuplicate(x,&mt->work); 
116:       PetscObjectDereference((PetscObject)mt->x); 
117:       mt->x = x;
118:       PetscObjectReference((PetscObject)mt->x); 
119:     }


122:     if (ls->bounded) {
123:         /* Compute step length needed to make all variables equal a bound */
124:         /* Compute the smallest steplength that will make one nonbinding variable
125:            equal the bound */
126:         VecGetLocalSize(ls->upper,&n1); 
127:         VecGetLocalSize(mt->x, &n2); 
128:         VecGetSize(ls->upper,&nn1); 
129:         VecGetSize(mt->x,&nn2); 
130:         if (n1 != n2 || nn1 != nn2) {
131:             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Variable vector not compatible with bounds vector");
132:         }
133:         VecScale(s,-1.0); 
134:         VecBoundGradientProjection(s,x,ls->lower,ls->upper,s); 
135:         VecScale(s,-1.0); 
136:         VecStepBoundInfo(x,ls->lower,ls->upper,s,&bstepmin1,&bstepmin2,&bstepmax);
137:         
138:         ls->stepmax = PetscMin(bstepmax,1.0e15);
139:     }

141:     VecDot(g,s,&dginit);
142:     
143:     if (PetscIsInfOrNanReal(dginit)) {
144:       PetscInfo1(ls,"Initial Line Search step * g is Inf or Nan (%G)\n",dginit); 
145:       ls->reason=TAOLINESEARCH_FAILED_INFORNAN;
146:       return(0);
147:     }
148:     if (dginit >= 0.0) {
149:       PetscInfo1(ls,"Initial Line Search step * g is not descent direction (%G)\n",dginit); 
150:       ls->reason = TAOLINESEARCH_FAILED_ASCENT;
151:       return(0);
152:     }


155:     /* Initialization */
156:     mt->bracket = 0;
157:     stage1 = 1;
158:     finit = *f;
159:     dgtest = ls->ftol * dginit;
160:     width = ls->stepmax - ls->stepmin;
161:     width1 = width * 2.0;
162:     VecCopy(x,mt->work);
163:     /* Variable dictionary:  
164:      stx, fx, dgx - the step, function, and derivative at the best step
165:      sty, fy, dgy - the step, function, and derivative at the other endpoint 
166:                    of the interval of uncertainty
167:      step, f, dg - the step, function, and derivative at the current step */

169:     stx = 0.0;
170:     fx  = finit;
171:     dgx = dginit;
172:     sty = 0.0;
173:     fy  = finit;
174:     dgy = dginit;
175:  
176:     ls->step=ls->initstep;
177:     for (i=0; i< ls->max_funcs; i++) {
178:     /* Set min and max steps to correspond to the interval of uncertainty */
179:       if (mt->bracket) {
180:         ls->stepmin = PetscMin(stx,sty); 
181:         ls->stepmax = PetscMax(stx,sty); 
182:       } 
183:       else {
184:         ls->stepmin = stx;
185:         ls->stepmax = ls->step + xtrapf * (ls->step - stx);
186:       }

188:       /* Force the step to be within the bounds */
189:       ls->step = PetscMax(ls->step,ls->stepmin);
190:       ls->step = PetscMin(ls->step,ls->stepmax);
191:     
192:       /* If an unusual termination is to occur, then let step be the lowest
193:          point obtained thus far */
194:       if ((stx!=0) && (((mt->bracket) && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) ||
195:         ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) ||
196:                        ((ls->nfeval+ls->nfgeval) >= ls->max_funcs - 1) || (mt->infoc == 0))) {
197:         ls->step = stx;
198:       }

200:       VecCopy(x,mt->work); 
201:       VecAXPY(mt->work,ls->step,s);    /* W = X + step*S */

203:       if (ls->bounded) {
204:           VecMedian(ls->lower, mt->work, ls->upper, mt->work);
205:       }
206:       if (ls->usegts) {
207:         TaoLineSearchComputeObjectiveAndGTS(ls,mt->work,f,&dg);  
208:         g_computed=PETSC_FALSE;
209:       } else {
210:         TaoLineSearchComputeObjectiveAndGradient(ls,mt->work,f,g); 
211:         g_computed=PETSC_TRUE;
212:         if (ls->bounded) {
213:           VecDot(g,x,&dg); 
214:           VecDot(g,mt->work,&dg2); 
215:           dg = (dg2 - dg)/ls->step;
216:         } else {
217:           VecDot(g,s,&dg); 
218:         }
219:       }

221:       if (0 == i) {
222:         ls->f_fullstep=*f;
223:       }



227:       if (PetscIsInfOrNanReal(*f) || PetscIsInfOrNanReal(dg)) {
228:         /* User provided compute function generated Not-a-Number, assume 
229:          domain violation and set function value and directional
230:          derivative to infinity. */
231:         *f = TAO_INFINITY;
232:         dg = TAO_INFINITY;
233:       }

235:       ftest1 = finit + ls->step * dgtest;
236:       if (ls->bounded) {
237:           ftest2 = finit + ls->step * dgtest * ls->ftol;
238:       }
239:       /* Convergence testing */
240:       if (((*f - ftest1 <= 1.0e-10 * PetscAbsReal(finit)) && 
241:            (PetscAbsReal(dg) + ls->gtol*dginit <= 0.0))) {
242:         PetscInfo(ls, "Line search success: Sufficient decrease and directional deriv conditions hold\n"); 
243:         ls->reason = TAOLINESEARCH_SUCCESS;
244:         break;
245:       }

247:       /* Check Armijo if beyond the first breakpoint */
248:       if (ls->bounded && (*f <= ftest2) && (ls->step >= bstepmin2)) {
249:           PetscInfo(ls,"Line search success: Sufficient decrease.\n"); 
250:           break;
251:       }

253:       /* Checks for bad cases */
254:       if (((mt->bracket) && (ls->step <= ls->stepmin||ls->step >= ls->stepmax)) || (!mt->infoc)) {
255:         PetscInfo(ls,"Rounding errors may prevent further progress.  May not be a step satisfying\n"); 
256:         PetscInfo(ls,"sufficient decrease and curvature conditions. Tolerances may be too small.\n"); 
257:         ls->reason = TAOLINESEARCH_HALTED_OTHER;
258:         break;
259:       }
260:       if ((ls->step == ls->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
261:         PetscInfo1(ls,"Step is at the upper bound, stepmax (%G)\n",ls->stepmax); 
262:         ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
263:         break;
264:       }
265:       if ((ls->step == ls->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
266:         PetscInfo1(ls,"Step is at the lower bound, stepmin (%G)\n",ls->stepmin); 
267:         ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
268:         break;
269:       }
270:       if ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol*ls->stepmax)){
271:         PetscInfo1(ls,"Relative width of interval of uncertainty is at most rtol (%G)\n",ls->rtol); 
272:         ls->reason = TAOLINESEARCH_HALTED_RTOL;
273:         break;
274:       }

276:       /* In the first stage, we seek a step for which the modified function
277:          has a nonpositive value and nonnegative derivative */
278:       if ((stage1) && (*f <= ftest1) && (dg >= dginit * PetscMin(ls->ftol, ls->gtol))) {
279:         stage1 = 0;
280:       }

282:       /* A modified function is used to predict the step only if we
283:          have not obtained a step for which the modified function has a 
284:          nonpositive function value and nonnegative derivative, and if a
285:          lower function value has been obtained but the decrease is not
286:          sufficient */

288:       if ((stage1) && (*f <= fx) && (*f > ftest1)) {
289:         fm   = *f - ls->step * dgtest;        /* Define modified function */
290:         fxm  = fx - stx * dgtest;                /* and derivatives */
291:         fym  = fy - sty * dgtest;
292:         dgm  = dg - dgtest;
293:         dgxm = dgx - dgtest;
294:         dgym = dgy - dgtest;

296:         /* if (dgxm * (ls->step - stx) >= 0.0) */
297:         /* Update the interval of uncertainty and compute the new step */
298:         Tao_mcstep(ls,&stx,&fxm,&dgxm,&sty,&fym,&dgym,&ls->step,&fm,&dgm);
299:       
300:         fx  = fxm + stx * dgtest;        /* Reset the function and */
301:         fy  = fym + sty * dgtest;        /* gradient values */
302:         dgx = dgxm + dgtest; 
303:         dgy = dgym + dgtest; 
304:       } 
305:       else {
306:         /* Update the interval of uncertainty and compute the new step */
307:         Tao_mcstep(ls,&stx,&fx,&dgx,&sty,&fy,&dgy,&ls->step,f,&dg);
308:       }
309:     
310:       /* Force a sufficient decrease in the interval of uncertainty */
311:       if (mt->bracket) {
312:         if (PetscAbsReal(sty - stx) >= 0.66 * width1) ls->step = stx + 0.5*(sty - stx);
313:         width1 = width;
314:         width = PetscAbsReal(sty - stx);
315:       }
316:     }
317:     if ((ls->nfeval+ls->nfgeval) > ls->max_funcs) {
318:       PetscInfo2(ls,"Number of line search function evals (%D) > maximum (%D)\n",(ls->nfeval+ls->nfgeval),ls->max_funcs); 
319:       ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
320:     }
321:   
322:     /* Finish computations */
323:     PetscInfo2(ls,"%D function evals in line search, step = %G\n",(ls->nfeval+ls->nfgeval),ls->step); 
324:     
325:     /* Set new solution vector and compute gradient if needed */
326:     VecCopy(mt->work,x); 
327:     if (!g_computed) {
328:       TaoLineSearchComputeGradient(ls,mt->work,g); 
329:     }

331:     return(0);
332: }

337: PetscErrorCode TaoLineSearchCreate_MT(TaoLineSearch ls)
338: {
340:     TAOLINESEARCH_MT_CTX *ctx;
343:     PetscNewLog(ls,TAOLINESEARCH_MT_CTX,&ctx); 
344:     ctx->bracket=0;
345:     ctx->infoc=1;
346:     ls->data = (void*)ctx;
347:     ls->initstep = 1.0;
348:     ls->ops->setup=0; 
349:     ls->ops->apply=TaoLineSearchApply_MT;
350:     ls->ops->view =TaoLineSearchView_MT;
351:     ls->ops->destroy=TaoLineSearchDestroy_MT;
352:     ls->ops->setfromoptions=TaoLineSearchSetFromOptions_MT;

354:     return(0);
355: }



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

364:      subroutine mcstep

366:      the purpose of mcstep is to compute a safeguarded step for
367:      a linesearch and to update an interval of uncertainty for
368:      a minimizer of the function.

370:      the parameter stx contains the step with the least function
371:      value. the parameter stp contains the current step. it is
372:      assumed that the derivative at stx is negative in the
373:      direction of the step. if bracket is set true then a
374:      minimizer has been bracketed in an interval of uncertainty
375:      with endpoints stx and sty.

377:      the subroutine statement is

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

382:      where

384:        stx, fx, and dx are variables which specify the step,
385:          the function, and the derivative at the best step obtained
386:          so far. The derivative must be negative in the direction
387:          of the step, that is, dx and stp-stx must have opposite
388:          signs. On output these parameters are updated appropriately.

390:        sty, fy, and dy are variables which specify the step,
391:          the function, and the derivative at the other endpoint of
392:          the interval of uncertainty. On output these parameters are
393:          updated appropriately.

395:        stp, fp, and dp are variables which specify the step,
396:          the function, and the derivative at the current step.
397:          If bracket is set true then on input stp must be
398:          between stx and sty. On output stp is set to the new step.

400:        bracket is a logical variable which specifies if a minimizer
401:          has been bracketed.  If the minimizer has not been bracketed
402:          then on input bracket must be set false.  If the minimizer
403:          is bracketed then on output bracket is set true.

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

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

413:      subprograms called

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

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

420: */

424: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,
425:                               PetscReal *stx, PetscReal *fx, PetscReal *dx,
426:                               PetscReal *sty, PetscReal *fy, PetscReal *dy,
427:                               PetscReal *stp, PetscReal *fp, PetscReal *dp)
428: {
429:   TAOLINESEARCH_MT_CTX *mtP = (TAOLINESEARCH_MT_CTX *) ls->data;
430:   PetscReal gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
431:   PetscInt bound;


435:   /* Check the input parameters for errors */
436:   mtP->infoc = 0;
437:   if (mtP->bracket && (*stp <= PetscMin(*stx,*sty) || (*stp >= PetscMax(*stx,*sty)))) SETERRQ(PETSC_COMM_SELF,1,"bad stp in bracket");
438:   if (*dx * (*stp-*stx) >= 0.0) SETERRQ(PETSC_COMM_SELF,1,"dx * (stp-stx) >= 0.0");
439:   if (ls->stepmax < ls->stepmin) SETERRQ(PETSC_COMM_SELF,1,"stepmax > stepmin");

441:   /* Determine if the derivatives have opposite sign */
442:   sgnd = *dp * (*dx / PetscAbsReal(*dx));

444:   if (*fp > *fx) {
445:     /* Case 1: a higher function value.
446:      The minimum is bracketed. If the cubic step is closer
447:      to stx than the quadratic step, the cubic step is taken,
448:      else the average of the cubic and quadratic steps is taken. */

450:     mtP->infoc = 1;
451:     bound = 1;
452:     theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
453:     s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
454:     s = PetscMax(s,PetscAbsReal(*dp));
455:     gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
456:     if (*stp < *stx) gamma1 = -gamma1;
457:     /* Can p be 0?  Check */
458:     p = (gamma1 - *dx) + theta;
459:     q = ((gamma1 - *dx) + gamma1) + *dp;
460:     r = p/q;
461:     stpc = *stx + r*(*stp - *stx);
462:     stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);

464:     if (PetscAbsReal(stpc-*stx) < PetscAbsReal(stpq-*stx)) {
465:       stpf = stpc;
466:     } 
467:     else {
468:       stpf = stpc + 0.5*(stpq - stpc);
469:     }
470:     mtP->bracket = 1;
471:   }
472:   else if (sgnd < 0.0) {
473:     /* Case 2: A lower function value and derivatives of
474:      opposite sign. The minimum is bracketed. If the cubic
475:      step is closer to stx than the quadratic (secant) step,
476:      the cubic step is taken, else the quadratic step is taken. */

478:     mtP->infoc = 2;
479:     bound = 0;
480:     theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
481:     s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
482:     s = PetscMax(s,PetscAbsReal(*dp));
483:     gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
484:     if (*stp > *stx) gamma1 = -gamma1;
485:     p = (gamma1 - *dp) + theta;
486:     q = ((gamma1 - *dp) + gamma1) + *dx;
487:     r = p/q;
488:     stpc = *stp + r*(*stx - *stp);
489:     stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);

491:     if (PetscAbsReal(stpc-*stp) > PetscAbsReal(stpq-*stp)) {
492:       stpf = stpc;
493:     }
494:     else {
495:       stpf = stpq;
496:     }
497:     mtP->bracket = 1;
498:   }
499:   else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
500:     /* Case 3: A lower function value, derivatives of the
501:      same sign, and the magnitude of the derivative decreases.
502:      The cubic step is only used if the cubic tends to infinity
503:      in the direction of the step or if the minimum of the cubic
504:      is beyond stp. Otherwise the cubic step is defined to be
505:      either stepmin or stepmax. The quadratic (secant) step is also
506:      computed and if the minimum is bracketed then the the step
507:      closest to stx is taken, else the step farthest away is taken. */

509:     mtP->infoc = 3;
510:     bound = 1;
511:     theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
512:     s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
513:     s = PetscMax(s,PetscAbsReal(*dp));

515:     /* The case gamma1 = 0 only arises if the cubic does not tend
516:        to infinity in the direction of the step. */
517:     gamma1 = s*PetscSqrtScalar(PetscMax(0.0,PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s)));
518:     if (*stp > *stx) gamma1 = -gamma1;
519:     p = (gamma1 - *dp) + theta;
520:     q = (gamma1 + (*dx - *dp)) + gamma1;
521:     r = p/q;
522:     if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r*(*stx - *stp);
523:     else if (*stp > *stx)        stpc = ls->stepmax;
524:     else                         stpc = ls->stepmin;
525:     stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);

527:     if (mtP->bracket) {
528:       if (PetscAbsReal(*stp-stpc) < PetscAbsReal(*stp-stpq)) {
529:         stpf = stpc;
530:       } 
531:       else {
532:         stpf = stpq;
533:       }
534:     }
535:     else {
536:       if (PetscAbsReal(*stp-stpc) > PetscAbsReal(*stp-stpq)) {
537:         stpf = stpc;
538:       }
539:       else {
540:         stpf = stpq;
541:       }
542:     }
543:   }
544:   else {
545:     /* Case 4: A lower function value, derivatives of the
546:        same sign, and the magnitude of the derivative does
547:        not decrease. If the minimum is not bracketed, the step
548:        is either stpmin or stpmax, else the cubic step is taken. */

550:     mtP->infoc = 4;
551:     bound = 0;
552:     if (mtP->bracket) {
553:       theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
554:       s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dy));
555:       s = PetscMax(s,PetscAbsReal(*dp));
556:       gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dy/s)*(*dp/s));
557:       if (*stp > *sty) gamma1 = -gamma1;
558:       p = (gamma1 - *dp) + theta;
559:       q = ((gamma1 - *dp) + gamma1) + *dy;
560:       r = p/q;
561:       stpc = *stp + r*(*sty - *stp);
562:       stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);

564:       stpf = stpc;
565:     } 
566:     else if (*stp > *stx) {
567:       stpf = ls->stepmax;
568:     } 
569:     else {
570:       stpf = ls->stepmin;
571:     }
572:   }
573:   
574:   /* Update the interval of uncertainty.  This update does not
575:      depend on the new step or the case analysis above. */

577:   if (*fp > *fx) {
578:     *sty = *stp;
579:     *fy = *fp;
580:     *dy = *dp;
581:   } 
582:   else {
583:     if (sgnd < 0.0) {
584:       *sty = *stx;
585:       *fy = *fx;
586:       *dy = *dx;
587:     }
588:     *stx = *stp;
589:     *fx = *fp;
590:     *dx = *dp;
591:   }
592:   
593:   /* Compute the new step and safeguard it. */
594:   stpf = PetscMin(ls->stepmax,stpf);
595:   stpf = PetscMax(ls->stepmin,stpf);
596:   *stp = stpf;
597:   if (mtP->bracket && bound) {
598:     if (*sty > *stx) {
599:       *stp = PetscMin(*stx+0.66*(*sty-*stx),*stp);
600:     }
601:     else {
602:       *stp = PetscMax(*stx+0.66*(*sty-*stx),*stp);
603:     }
604:   }
605:   return(0);
606: }