Actual source code: morethuente.c

  1: /*$Id$*/

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

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

 10:      subroutine mcstep

 12:      the purpose of mcstep is to compute a safeguarded step for
 13:      a linesearch and to update an interval of uncertainty for
 14:      a minimizer of the function.

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

 23:      the subroutine statement is

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

 28:      where

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

 36:        sty, fy, and dy are variables which specify the step,
 37:          the function, and the derivative at the other endpoint of
 38:          the interval of uncertainty. On output these parameters are
 39:          updated appropriately.

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

 46:        bracket is a logical variable which specifies if a minimizer
 47:          has been bracketed.  If the minimizer has not been bracketed
 48:          then on input bracket must be set false.  If the minimizer
 49:          is bracketed then on output bracket is set true.

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

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

 59:      subprograms called

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

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

 66: */

 70: static int TaoStep_LineSearch(TAO_SOLVER tao,
 71:                               double *stx, double *fx, double *dx,
 72:                               double *sty, double *fy, double *dy,
 73:                               double *stp, double *fp, double *dp)
 74: {
 75:   TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
 76:   double gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
 77:   TaoInt bound;

 79:   TaoFunctionBegin;

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

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

 90:   if (*fp > *fx) {
 91:     // Case 1: a higher function value.
 92:     // The minimum is bracketed. If the cubic step is closer
 93:     // to stx than the quadratic step, the cubic step is taken,
 94:     // else the average of the cubic and quadratic steps is taken.

 96:     neP->infoc = 1;
 97:     bound = 1;
 98:     theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
 99:     s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dx));
100:     s = TaoMax(s,TaoAbsDouble(*dp));
101:     gamma1 = s*sqrt(pow(theta/s,2.0) - (*dx/s)*(*dp/s));
102:     if (*stp < *stx) gamma1 = -gamma1;
103:     /* Can p be 0?  Check */
104:     p = (gamma1 - *dx) + theta;
105:     q = ((gamma1 - *dx) + gamma1) + *dp;
106:     r = p/q;
107:     stpc = *stx + r*(*stp - *stx);
108:     stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);

110:     if (TaoAbsDouble(stpc-*stx) < TaoAbsDouble(stpq-*stx)) {
111:       stpf = stpc;
112:     } 
113:     else {
114:       stpf = stpc + 0.5*(stpq - stpc);
115:     }
116:     neP->bracket = 1;
117:   }
118:   else if (sgnd < 0.0) {
119:     // Case 2: A lower function value and derivatives of
120:     // opposite sign. The minimum is bracketed. If the cubic
121:     // step is closer to stx than the quadratic (secant) step,
122:     // the cubic step is taken, else the quadratic step is taken.

124:     neP->infoc = 2;
125:     bound = 0;
126:     theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
127:     s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dx));
128:     s = TaoMax(s,TaoAbsDouble(*dp));
129:     gamma1 = s*sqrt(pow(theta/s,2.0) - (*dx/s)*(*dp/s));
130:     if (*stp > *stx) gamma1 = -gamma1;
131:     p = (gamma1 - *dp) + theta;
132:     q = ((gamma1 - *dp) + gamma1) + *dx;
133:     r = p/q;
134:     stpc = *stp + r*(*stx - *stp);
135:     stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);

137:     if (TaoAbsDouble(stpc-*stp) > TaoAbsDouble(stpq-*stp)) {
138:       stpf = stpc;
139:     }
140:     else {
141:       stpf = stpq;
142:     }
143:     neP->bracket = 1;
144:   }
145:   else if (TaoAbsDouble(*dp) < TaoAbsDouble(*dx)) {
146:     // Case 3: A lower function value, derivatives of the
147:     // same sign, and the magnitude of the derivative decreases.
148:     // The cubic step is only used if the cubic tends to infinity
149:     // in the direction of the step or if the minimum of the cubic
150:     // is beyond stp. Otherwise the cubic step is defined to be
151:     // either stepmin or stepmax. The quadratic (secant) step is also
152:     // computed and if the minimum is bracketed then the the step
153:     // closest to stx is taken, else the step farthest away is taken.

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

161:     // The case gamma1 = 0 only arises if the cubic does not tend
162:     // to infinity in the direction of the step.
163:     gamma1 = s*sqrt(TaoMax(0.0,pow(theta/s,2.0) - (*dx/s)*(*dp/s)));
164:     if (*stp > *stx) gamma1 = -gamma1;
165:     p = (gamma1 - *dp) + theta;
166:     q = (gamma1 + (*dx - *dp)) + gamma1;
167:     r = p/q;
168:     if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r*(*stx - *stp);
169:     else if (*stp > *stx)        stpc = neP->stepmax;
170:     else                         stpc = neP->stepmin;
171:     stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);

173:     if (neP->bracket) {
174:       if (TaoAbsDouble(*stp-stpc) < TaoAbsDouble(*stp-stpq)) {
175:         stpf = stpc;
176:       } 
177:       else {
178:         stpf = stpq;
179:       }
180:     }
181:     else {
182:       if (TaoAbsDouble(*stp-stpc) > TaoAbsDouble(*stp-stpq)) {
183:         stpf = stpc;
184:       }
185:       else {
186:         stpf = stpq;
187:       }
188:     }
189:   }
190:   else {
191:     // Case 4: A lower function value, derivatives of the
192:     // same sign, and the magnitude of the derivative does
193:     // not decrease. If the minimum is not bracketed, the step
194:     // is either stpmin or stpmax, else the cubic step is taken.

196:     neP->infoc = 4;
197:     bound = 0;
198:     if (neP->bracket) {
199:       theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
200:       s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dy));
201:       s = TaoMax(s,TaoAbsDouble(*dp));
202:       gamma1 = s*sqrt(pow(theta/s,2.0) - (*dy/s)*(*dp/s));
203:       if (*stp > *sty) gamma1 = -gamma1;
204:       p = (gamma1 - *dp) + theta;
205:       q = ((gamma1 - *dp) + gamma1) + *dy;
206:       r = p/q;
207:       stpc = *stp + r*(*sty - *stp);
208:       stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);

210:       stpf = stpc;
211:     } 
212:     else if (*stp > *stx) {
213:       stpf = neP->stepmax;
214:     } 
215:     else {
216:       stpf = neP->stepmin;
217:     }
218:   }
219:   
220:   // Update the interval of uncertainty.  This update does not
221:   // depend on the new step or the case analysis above.

223:   if (*fp > *fx) {
224:     *sty = *stp;
225:     *fy = *fp;
226:     *dy = *dp;
227:   } 
228:   else {
229:     if (sgnd < 0.0) {
230:       *sty = *stx;
231:       *fy = *fx;
232:       *dy = *dx;
233:     }
234:     *stx = *stp;
235:     *fx = *fp;
236:     *dx = *dp;
237:   }
238:   
239:   // Compute the new step and safeguard it.
240:   stpf = TaoMin(neP->stepmax,stpf);
241:   stpf = TaoMax(neP->stepmin,stpf);
242:   *stp = stpf;
243:   if (neP->bracket && bound) {
244:     if (*sty > *stx) {
245:       *stp = TaoMin(*stx+0.66*(*sty-*stx),*stp);
246:     }
247:     else {
248:       *stp = TaoMax(*stx+0.66*(*sty-*stx),*stp);
249:     }
250:   }
251:   TaoFunctionReturn(0);
252: }

256: /* @ 
257:    TaoApply_LineSearch - This routine performs a line search algorithm.

259:    Input Parameters:
260: +  tao - TAO_SOLVER context
261: .  X - current iterate (on output X contains new iterate, X + step*S)
262: .  S - search direction
263: .  f - objective function evaluated at X
264: .  G - gradient evaluated at X
265: .  W - work vector
266: -  step - initial estimate of step length

268:    Output parameters:
269: +  f - objective function evaluated at new iterate, X + step*S
270: .  G - gradient evaluated at new iterate, X + step*S
271: .  X - new iterate
272: .  gnorm - 2-norm of G
273: -  step - final step length


276:    Info is set to one of:
277: +   0 - the line search succeeds; the sufficient decrease
278:    condition and the directional derivative condition hold

280:    negative number if an input parameter is invalid
281: .   -1 -  step < 0 
282: .   -2 -  ftol < 0 
283: .   -3 -  rtol < 0 
284: .   -4 -  gtol < 0 
285: .   -5 -  stepmin < 0 
286: .   -6 -  stepmax < stepmin
287: -   -7 -  maxfev < 0

289:    positive number > 1 if the line search otherwise terminates
290: +    2 -  Relative width of the interval of uncertainty is 
291:          at most rtol.
292: .    3 -  Maximum number of function evaluations (maxfev) has 
293:          been reached.
294: .    4 -  Step is at the lower bound, stepmin.
295: .    5 -  Step is at the upper bound, stepmax.
296: .    6 -  Rounding errors may prevent further progress. 
297:          There may not be a step that satisfies the 
298:          sufficient decrease and curvature conditions.  
299:          Tolerances may be too small.
300: -    7 -  Search direction is not a descent direction.

302:    Notes:
303:    This algorithm is taken from More' and Thuente, "Line search algorithms
304:    with guaranteed sufficient decrease", Argonne National Laboratory, 
305:    Technical Report MCS-P330-1092.

307:    Notes:
308:    This routine is used within the following TAO unconstrained minimization 
309:    solvers: Newton linesearch (tao_nls), limited memory variable metric
310:    (tao_lmvm), and conjugate gradient (tao_cg).

312:    Level: advanced

314: .keywords: TAO_SOLVER, linesearch
315: @ */
316: static int TaoApply_LineSearch(TAO_SOLVER tao,TaoVec* X,TaoVec* G,TaoVec* S,TaoVec* W,double *f, double *f_full,
317:                         double *step,TaoInt *info2,void*ctx)
318: {
319:   TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
320:   double    xtrapf = 4.0;
321:   double    finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
322:   double    dgx, dgy, dg, fx, fy, stx, sty, dgtest, ftest1=0.0;
323:   int       info;
324:   TaoInt  i, stage1;

326: #if defined(PETSC_USE_COMPLEX)
327:   PetscScalar    cdginit, cdg, cstep = 0.0;
328: #endif

330:   TaoFunctionBegin;
331:   /* neP->stepmin - lower bound for step */
332:   /* neP->stepmax - upper bound for step */
333:   /* neP->rtol           - relative tolerance for an acceptable step */
334:   /* neP->ftol           - tolerance for sufficient decrease condition */
335:   /* neP->gtol           - tolerance for curvature condition */
336:   /* neP->nfev           - number of function evaluations */
337:   /* neP->maxfev  - maximum number of function evaluations */

339:   *info2 = 0;

341:   /* Check input parameters for errors */
342:   if (*step < 0.0) {
343:     info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: step (%g) < 0\n",*step); CHKERRQ(info);
344:     *info2 = -1;
345:   } 
346:   
347:   if (neP->ftol < 0.0) {
348:     info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: ftol (%g) < 0\n",neP->ftol); CHKERRQ(info);
349:     *info2 = -2;
350:   } 
351:   
352:   if (neP->rtol < 0.0) {
353:     info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: rtol (%g) < 0\n",neP->rtol); CHKERRQ(info);
354:     *info2 = -3;
355:   } 

357:   if (neP->gtol < 0.0) {
358:     info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: gtol (%g) < 0\n",neP->gtol); CHKERRQ(info);
359:     *info2 = -4;
360:   } 

362:   if (neP->stepmin < 0.0) {
363:     info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: stepmin (%g) < 0\n",neP->stepmin); CHKERRQ(info);
364:     *info2 = -5;
365:   } 

367:   if (neP->stepmax < neP->stepmin) {
368:     info = PetscInfo2(tao, "TaoApply_LineSearch: Line search error: stepmax (%g) < stepmin (%g)\n", neP->stepmax,neP->stepmin); CHKERRQ(info);
369:     *info2 = -6;
370:   }
371:   
372:   if (neP->maxfev < 0) {
373:     info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: maxfev (%d) < 0\n",neP->maxfev); CHKERRQ(info);
374:     *info2 = -7;
375:   }

377:   if (TaoInfOrNaN(*f)) {
378:     info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: function (%g) inf or nan\n", *f); CHKERRQ(info);
379:     *info2 = -8;
380:   }

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

384: #if defined(PETSC_USE_COMPLEX)
385:   info = G->Dot(S,&cdginit);CHKERRQ(info); dginit = TaoReal(cdginit);
386: #else
387:   info = G->Dot(S,&dginit);CHKERRQ(info);
388: #endif

390:   if (TaoInfOrNaN(dginit)) {
391:     info = PetscInfo1(tao,"TaoApply_LineSearch: Line search error: dginit (%g) inf or nan\n", dginit); CHKERRQ(info);
392:     *info2 = -9;
393:   }

395:   if (dginit >= 0.0) {
396:     info = PetscInfo(tao,"TaoApply_LineSearch:Search direction not a descent direction\n"); CHKERRQ(info);
397:     *info2 = -10;
398:   }

400:   if (*info2) {
401:     TaoFunctionReturn(0);
402:   }

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

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

439:     /* Force the step to be within the bounds */
440:     *step = TaoMax(*step,neP->stepmin);
441:     *step = TaoMin(*step,neP->stepmax);
442:     
443:     /* If an unusual termination is to occur, then let step be the lowest
444:        point obtained thus far */
445:     if (((neP->bracket) && (*step <= neP->stepmin || *step >= neP->stepmax)) ||
446:         ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol * neP->stepmax)) ||
447:         (neP->nfev >= neP->maxfev - 1) || (neP->infoc == 0)) {
448:       *step = stx;
449:     }

451: #if defined(PETSC_USE_COMPLEX)
452:     cstep = *step;
453:     info = X->Waxpby(cstep,S,1.0,W);CHKERRQ(info);
454: #else
455:     info = X->Waxpby(*step,S,1.0,W);CHKERRQ(info);         /* X = W + step*S */
456: #endif
457:     tao->current_step=*step;
458:     info = TaoComputeMeritFunctionGradient(tao,X,f,G);CHKERRQ(info);
459:     tao->new_search=TAO_FALSE;
460:     if (0 == i) {
461:       *f_full = *f;
462:     }

464:     neP->nfev++;
465: #if defined(PETSC_USE_COMPLEX)
466:     info = G->Dot(S,&cdg);CHKERRQ(info); dg = TaoReal(cdg);
467: #else
468:     info = G->Dot(S,&dg);CHKERRQ(info);                /* dg = G^T S */
469: #endif

471:     if (TaoInfOrNaN(*f) || TaoInfOrNaN(dg)) {
472:       // User provided compute function generated Not-a-Number, assume 
473:       // domain violation and set function value and directional
474:       // derivative to infinity.
475:       *f = TAO_INFINITY;
476:       dg = TAO_INFINITY;
477:     }

479:     ftest1 = finit + *step * dgtest;

481:     /* Convergence testing */
482:     if (((*f - ftest1 <= 1.0e-10 * fabs(finit)) && 
483:          (TaoAbsDouble(dg) + neP->gtol*dginit <= 0.0))) {
484:       info = PetscInfo(tao, "TaoApply_LineSearch:Line search success: Sufficient decrease and directional deriv conditions hold\n"); CHKERRQ(info);
485:       *info2 = 0; 
486:       break;
487:     }

489:     /* Checks for bad cases */
490:     if (((neP->bracket) && (*step <= neP->stepmin||*step >= neP->stepmax)) || (!neP->infoc)) {
491:       info = PetscInfo(tao,"TaoApply_LineSearch:Rounding errors may prevent further progress.  May not be a step satisfying\n"); CHKERRQ(info);
492:       info = PetscInfo(tao,"TaoApply_LineSearch:sufficient decrease and curvature conditions. Tolerances may be too small.\n"); CHKERRQ(info);
493:       *info2 = 6; break;
494:     }
495:     if ((*step == neP->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
496:       info = PetscInfo1(tao,"TaoApply_LineSearch:Step is at the upper bound, stepmax (%g)\n",neP->stepmax); CHKERRQ(info);
497:       *info2 = 5; break;
498:     }
499:     if ((*step == neP->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
500:       info = PetscInfo1(tao,"TaoApply_LineSearch:Step is at the lower bound, stepmin (%g)\n",neP->stepmin); CHKERRQ(info);
501:       *info2 = 4; break;
502:     }
503:     if (neP->nfev >= neP->maxfev) {
504:       info = PetscInfo2(tao,"TaoApply_LineSearch:Number of line search function evals (%d) > maximum (%d)\n",neP->nfev,neP->maxfev); CHKERRQ(info);
505:       *info2 = 3; break;
506:     }
507:     if ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol*neP->stepmax)){
508:       info = PetscInfo1(tao,"TaoApply_LineSearch:Relative width of interval of uncertainty is at most rtol (%g)\n",neP->rtol); CHKERRQ(info);
509:       *info2 = 2; break;
510:     }

512:     /* In the first stage, we seek a step for which the modified function
513:        has a nonpositive value and nonnegative derivative */
514:     if ((stage1) && (*f <= ftest1) && (dg >= dginit * TaoMin(neP->ftol, neP->gtol))) {
515:       stage1 = 0;
516:     }

518:     /* A modified function is used to predict the step only if we
519:        have not obtained a step for which the modified function has a 
520:        nonpositive function value and nonnegative derivative, and if a
521:        lower function value has been obtained but the decrease is not
522:        sufficient */

524:     if ((stage1) && (*f <= fx) && (*f > ftest1)) {
525:       fm   = *f - *step * dgtest;        /* Define modified function */
526:       fxm  = fx - stx * dgtest;                /* and derivatives */
527:       fym  = fy - sty * dgtest;
528:       dgm  = dg - dgtest;
529:       dgxm = dgx - dgtest;
530:       dgym = dgy - dgtest;
531:       
532:       /* Update the interval of uncertainty and compute the new step */
533:       info = TaoStep_LineSearch(tao,&stx,&fxm,&dgxm,&sty,&fym,&dgym,step,&fm,&dgm);CHKERRQ(info);
534:       
535:       fx  = fxm + stx * dgtest;        /* Reset the function and */
536:       fy  = fym + sty * dgtest;        /* gradient values */
537:       dgx = dgxm + dgtest; 
538:       dgy = dgym + dgtest; 
539:     } 
540:     else {
541:       /* Update the interval of uncertainty and compute the new step */
542:       info = TaoStep_LineSearch(tao,&stx,&fx,&dgx,&sty,&fy,&dgy,step,f,&dg);CHKERRQ(info);
543:     }
544:     
545:     /* Force a sufficient decrease in the interval of uncertainty */
546:     if (neP->bracket) {
547:       if (TaoAbsDouble(sty - stx) >= 0.66 * width1) *step = stx + 0.5*(sty - stx);
548:       width1 = width;
549:       width = TaoAbsDouble(sty - stx);
550:     }
551:   }
552:   
553:   /* Finish computations */
554:   info = PetscInfo2(tao,"TaoApply_LineSearch:%d function evals in line search, step = %10.4f\n",neP->nfev,*step); CHKERRQ(info);
555:   TaoFunctionReturn(0);
556: }

560: /* @ 
561:    TaoApply_BoundLineSearch - This routine performs a line search algorithm.

563:    Input Parameters:
564: +  tao - TAO_SOLVER context
565: .  X - current iterate (on output X contains new iterate, X + step*S)
566: .  S - search direction
567: .  f - objective function evaluated at X
568: .  G - gradient evaluated at X
569: .  W - work vector
570: -  step - initial estimate of step length

572:    Output parameters:
573: +  f - objective function evaluated at new iterate, X + step*S
574: .  G - gradient evaluated at new iterate, X + step*S
575: .  X - new iterate
576: .  gnorm - 2-norm of G
577: -  step - final step length


580:    Info is set to one of:
581: +   0 - the line search succeeds; the sufficient decrease
582:    condition and the directional derivative condition hold

584:    negative number if an input parameter is invalid
585: .   -1 -  step < 0 
586: .   -2 -  ftol < 0 
587: .   -3 -  rtol < 0 
588: .   -4 -  gtol < 0 
589: .   -5 -  stepmin < 0 
590: .   -6 -  stepmax < stepmin
591: -   -7 -  maxfev < 0

593:    positive number > 1 if the line search otherwise terminates
594: +    2 -  Relative width of the interval of uncertainty is 
595:          at most rtol.
596: .    3 -  Maximum number of function evaluations (maxfev) has 
597:          been reached.
598: .    4 -  Step is at the lower bound, stepmin.
599: .    5 -  Step is at the upper bound, stepmax.
600: .    6 -  Rounding errors may prevent further progress. 
601:          There may not be a step that satisfies the 
602:          sufficient decrease and curvature conditions.  
603:          Tolerances may be too small.
604: -    7 -  Search direction is not a descent direction.

606:    Notes:
607:    This algorithm is is a modification of the algorithm by More' and
608:    Thuente.  The modifications concern bounds.  This algorithm steps
609:    in the direction passed into this routine.  This point get
610:    projected back into the feasible set.  In the context of bound
611:    constrained optimization, there may not be a point in the piecewise
612:    linear path that satisfies the Wolfe conditions.  When the active
613:    set is changing, decrease in the objective function may be
614:    sufficient to terminate this line search.

616:    Note: Much of this coded is identical to the More' Thuente line search.
617:    Variations to the code are commented.

619:    Notes:
620:    This routine is used within the following TAO bound constrained minimization
621:    solvers: Newton linesearch (tao_tron) and limited memory variable metric
622:    (tao_blmvm).

624:    Level: advanced

626: .keywords: TAO_SOLVER, linesearch
627: @ */
628: static int TaoApply_BoundLineSearch(TAO_SOLVER tao,TaoVec* X,TaoVec* G,TaoVec* S,TaoVec* W,double *f, double *f_full,
629:                         double *step,TaoInt *info2,void*ctx)
630: {
631:   TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
632:   TaoVec *XL,*XU;
633:   double    xtrapf = 4.0;
634:   double    finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
635:   double    dgx, dgy, dg, fx, fy, stx, sty, dgtest, ftest1=0.0, ftest2=0.0;
636:   double    bstepmin1, bstepmin2, bstepmax;
637:   double    dg1, dg2;
638:   int       info;
639:   TaoInt    i, stage1;

641: #if defined(PETSC_USE_COMPLEX)
642:   PetscScalar    cdginit, cdg, cstep = 0.0;
643: #endif

645:   TaoFunctionBegin;
646:   /* neP->stepmin - lower bound for step */
647:   /* neP->stepmax - upper bound for step */
648:   /* neP->rtol           - relative tolerance for an acceptable step */
649:   /* neP->ftol           - tolerance for sufficient decrease condition */
650:   /* neP->gtol           - tolerance for curvature condition */
651:   /* neP->nfev           - number of function evaluations */
652:   /* neP->maxfev  - maximum number of function evaluations */

654:   /* Check input parameters for errors */
655:   if (*step < 0.0) {
656:     info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: step (%g) < 0\n",*step); CHKERRQ(info);
657:     *info2 = -1; TaoFunctionReturn(0);
658:   } 
659:   else if (neP->ftol < 0.0) {
660:     info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: ftol (%g) < 0\n",neP->ftol); CHKERRQ(info);
661:     *info2 = -2; TaoFunctionReturn(0);
662:   } 
663:   else if (neP->rtol < 0.0) {
664:     info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: rtol (%g) < 0\n",neP->rtol); CHKERRQ(info);
665:     *info2 = -3; TaoFunctionReturn(0);
666:   } 
667:   else if (neP->gtol < 0.0) {
668:     info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: gtol (%g) < 0\n",neP->gtol); CHKERRQ(info);
669:     *info2 = -4; TaoFunctionReturn(0);
670:   } 
671:   else if (neP->stepmin < 0.0) {
672:     info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: stepmin (%g) < 0\n",neP->stepmin); CHKERRQ(info);
673:     *info2 = -5; TaoFunctionReturn(0);
674:   }
675:   else if (neP->stepmax < neP->stepmin) {
676:     info = PetscInfo2(tao,"TaoApply_BoundLineSearch:Line search error: stepmax (%g) < stepmin (%g)\n", neP->stepmax,neP->stepmin); CHKERRQ(info);
677:     *info2 = -6; TaoFunctionReturn(0);
678:   } 
679:   else if (neP->maxfev < 0.0) {
680:     info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: maxfev (%d) < 0\n",neP->maxfev); CHKERRQ(info);
681:     *info2 = -7; TaoFunctionReturn(0);
682:   }

684:   /* Compute step length needed to make all variables equal a bound */ 
685:   /* Compute the smallest steplength that will make one nonbinding  */
686:   /* variable equal the bound */ 
687:   info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
688:   info = S->Negate(); CHKERRQ(info);
689:   info = S->BoundGradientProjection(S,XL,X,XU); CHKERRQ(info);
690:   info = S->Negate(); CHKERRQ(info);
691:   info = X->StepBoundInfo(XL,XU,S,&bstepmin1,&bstepmin2,&bstepmax); CHKERRQ(info);
692:   neP->stepmax=TaoMin(bstepmax,1.0e+15);

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

696: #if defined(PETSC_USE_COMPLEX)
697:   info = G->Dot(S,&cdginit);CHKERRQ(info); dginit = TaoReal(cdginit);
698: #else
699:   info = G->Dot(S,&dginit);CHKERRQ(info);
700: #endif

702:   if (dginit >= 0.0) {
703:     info = PetscInfo(tao,"TaoApply_BoundLineSearch:Search direction not a descent direction\n"); CHKERRQ(info);
704:     *info2 = 7; TaoFunctionReturn(0);
705:   }

707:   /* Initialization */
708:   neP->bracket = 0;
709:   *info2  = 0;
710:   stage1  = 1;
711:   finit   = *f;
712:   dgtest  = neP->ftol * dginit;
713:   width   = neP->stepmax - neP->stepmin;
714:   width1  = width * 2.0;
715:   info = W->CopyFrom(X);CHKERRQ(info);
716:   /* Variable dictionary:  
717:      stx, fx, dgx - the step, function, and derivative at the best step
718:      sty, fy, dgy - the step, function, and derivative at the other endpoint 
719:                    of the interval of uncertainty
720:      step, f, dg - the step, function, and derivative at the current step */

722:   stx = 0.0;
723:   fx  = finit;
724:   dgx = dginit;
725:   sty = 0.0;
726:   fy  = finit;
727:   dgy = dginit;
728:  
729:   neP->nfev = 0;
730:   tao->new_search=TAO_TRUE;
731:   for (i=0; i< neP->maxfev; i++) {
732:     /* Set min and max steps to correspond to the interval of uncertainty */
733:     if (neP->bracket) {
734:       neP->stepmin = TaoMin(stx,sty); 
735:       neP->stepmax = TaoMax(stx,sty); 
736:     } else {
737:       neP->stepmin = stx;
738:       neP->stepmax = *step + xtrapf * (*step - stx);
739:     }

741:     /* Force the step to be within the bounds */
742:     *step = TaoMax(*step,neP->stepmin);
743:     *step = TaoMin(*step,neP->stepmax);

745:     /* If an unusual termination is to occur, then let step be the lowest
746:        point obtained thus far */
747:     if (((neP->bracket) && (*step <= neP->stepmin || *step >= neP->stepmax)) ||
748:         ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol * neP->stepmax)) ||
749:         (neP->nfev >= neP->maxfev - 1) || (neP->infoc == 0)) {
750:       *step = stx;
751:     }
752:  
753: #if defined(PETSC_USE_COMPLEX)
754:     cstep = *step;
755:     info = X->Waxpby(cstep,S,1.0,W);CHKERRQ(info);
756: #else
757:     info = X->Waxpby(*step,S,1.0,W);CHKERRQ(info);         /* X = W + step*S */
758: #endif

760:     info=X->Median(XL,X,XU);CHKERRQ(info);
761:     tao->current_step=*step;
762:     info = TaoComputeMeritFunctionGradient(tao,X,f,G);CHKERRQ(info);
763:     tao->new_search=TAO_FALSE;
764:     if (0 == i) {
765:       *f_full = *f;
766:     }

768:     neP->nfev++;
769: #if defined(PETSC_USE_COMPLEX)
770:     info = G->Dot(S,&cdg);CHKERRQ(info); dg = TaoReal(cdg);
771: #else
772:     info = G->Dot(W,&dg1);CHKERRQ(info);                /* dg = G^T S */
773:     info = G->Dot(X,&dg2);CHKERRQ(info);                /* dg = G^T S */
774:     dg = (dg2-dg1) / (*step);
775: #endif

777:     if ((*f != *f) || (dg1 != dg1) || (dg2 != dg2) || (dg != dg)) {
778:       // User provided compute function generated Not-a-Number, assume
779:       // domain violation and set function value and directional
780:       // derivative to infinity.
781:       *f = TAO_INFINITY;
782:       dg = TAO_INFINITY;
783:       dg1 = TAO_INFINITY;
784:       dg2 = TAO_INFINITY;
785:     }

787:     ftest1 = finit + (*step) * dgtest;
788:     ftest2 = finit + (*step) * dgtest * neP->ftol;        // Armijo

790:     /* Convergence testing */
791:     if ((*f <= ftest1) && (TaoAbsDouble(dg) <= neP->gtol*(-dginit))) {
792:       info = PetscInfo(tao, "TaoApply_BoundLineSearch:Line search success: Sufficient decrease and directional deriv conditions hold\n"); CHKERRQ(info);
793:       *info2 = 0; 
794:       break;
795:     }

797:     /* Check Armijo if beyond the first breakpoint */
798:     if ((*f <= ftest2) && (*step >= bstepmin2)) {
799:       info = PetscInfo(tao,"TaoApply_BoundLineSearch:Line search success: Sufficient decrease\n"); CHKERRQ(info);
800:       *info2 = 0; 
801:       break;
802:     }

804:     /* Checks for bad cases */
805:     if (((neP->bracket) && (*step <= neP->stepmin||*step >= neP->stepmax)) || (!neP->infoc)) {
806:       info = PetscInfo(tao,"TaoApply_LineSearch:Rounding errors may prevent further progress.  May not be a step satisfying\n"); CHKERRQ(info);
807:       info = PetscInfo(tao,"TaoApply_BoundLineSearch:sufficient decrease and curvature conditions. Tolerances may be too small.\n"); CHKERRQ(info);
808:       *info2 = 6; break;
809:     }
810:     if ((*step == neP->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
811:       info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Step is at the upper bound, stepmax (%g)\n",neP->stepmax); CHKERRQ(info);
812:       *info2 = 5; break;
813:     }
814:     if ((*step == neP->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
815:       info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Step is at the lower bound, stepmin (%g)\n",neP->stepmin); CHKERRQ(info);
816:       *info2 = 4; break;
817:     }
818:     if (neP->nfev >= neP->maxfev) {
819:       info = PetscInfo2(tao,"TaoApply_BoundLineSearch:Number of line search function evals (%d) > maximum (%d)\n",neP->nfev,neP->maxfev); CHKERRQ(info);
820:       *info2 = 3; break;
821:     }
822:     if ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol*neP->stepmax)){
823:       info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Relative width of interval of uncertainty is at most rtol (%g)\n",neP->rtol); CHKERRQ(info);
824:       *info2 = 2; break;
825:     }

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

832:     /* A modified function is used to predict the step only if we
833:        have not obtained a step for which the modified function has a 
834:        nonpositive function value and nonnegative derivative, and if a
835:        lower function value has been obtained but the decrease is not
836:        sufficient */

838:     if ((stage1) && (*f <= fx) && (*f > ftest1)) {
839:       fm   = *f - *step * dgtest;        /* Define modified function */
840:       fxm  = fx - stx * dgtest;                /* and derivatives */
841:       fym  = fy - sty * dgtest;
842:       dgm  = dg - dgtest;
843:       dgxm = dgx - dgtest;
844:       dgym = dgy - dgtest;

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

849:       fx  = fxm + stx * dgtest;        /* Reset the function and */
850:       fy  = fym + sty * dgtest;        /* gradient values */
851:       dgx = dgxm + dgtest; 
852:       dgy = dgym + dgtest; 
853:     } else {
854:       /* Update the interval of uncertainty and compute the new step */
855:       info = TaoStep_LineSearch(tao,&stx,&fx,&dgx,&sty,&fy,&dgy,step,f,&dg);CHKERRQ(info);
856:     }

858:     /* Force a sufficient decrease in the interval of uncertainty */
859:     if (neP->bracket) {
860:       if (TaoAbsDouble(sty - stx) >= 0.66 * width1) *step = stx + 0.5*(sty - stx);
861:       width1 = width;
862:       width = TaoAbsDouble(sty - stx);
863:     }
864:   }

866:   /* Finish computations */
867:   info = PetscInfo2(tao,"TaoApply_BoundLineSearch:%d function evals in line search, step = %10.4e\n",neP->nfev,*step); CHKERRQ(info);
868:   TaoFunctionReturn(0);
869: }

873: static int TaoDestroy_LineSearch(TAO_SOLVER tao, void *ctx)
874: {
875:   int  info;

877:   TaoFunctionBegin;
878:   info = TaoFree(tao->linectx);CHKERRQ(info);
879:   TaoFunctionReturn(0);
880: }

884: static int TaoSetOptions_LineSearch(TAO_SOLVER tao, void*linectx)
885: {
886:   TAO_LINESEARCH *ctx = (TAO_LINESEARCH *)linectx;
887:   int            info;

889:   TaoFunctionBegin;
890:   info = TaoOptionsHead("More-Thuente line line search options for minimization");CHKERRQ(info);
891:   info = TaoOptionInt("-tao_ls_maxfev","max function evals in line search","",ctx->maxfev,&ctx->maxfev,0);CHKERRQ(info);
892:   info = TaoOptionDouble("-tao_ls_ftol","tol for sufficient decrease","",ctx->ftol,&ctx->ftol,0);CHKERRQ(info);
893:   info = TaoOptionDouble("-tao_ls_gtol","tol for curvature condition","",ctx->gtol,&ctx->gtol,0);CHKERRQ(info);
894:   info = TaoOptionDouble("-tao_ls_rtol","relative tol for acceptable step","",ctx->rtol,&ctx->rtol,0);CHKERRQ(info);
895:   info = TaoOptionDouble("-tao_ls_stepmin","lower bound for step","",ctx->stepmin,&ctx->stepmin,0);CHKERRQ(info);
896:   info = TaoOptionDouble("-tao_ls_stepmax","upper bound for step","",ctx->stepmax,&ctx->stepmax,0);CHKERRQ(info);
897:   info = TaoOptionsTail();CHKERRQ(info);
898:   TaoFunctionReturn(0);
899: }

903: static int TaoView_LineSearch(TAO_SOLVER tao, void*ctx)
904: {
905:   TAO_LINESEARCH *ls = (TAO_LINESEARCH *)ctx;
906:   int            info;

908:   TaoFunctionBegin;
909:   info=TaoPrintInt(tao,"    More'-Thuente Line search: maxf=%d,",ls->maxfev);CHKERRQ(info);
910:   info=TaoPrintDouble(tao," ftol=%g,",ls->ftol);CHKERRQ(info);
911:   info=TaoPrintDouble(tao," rtol=%g,",ls->rtol);CHKERRQ(info);
912:   info=TaoPrintDouble(tao," gtol=%g\n",ls->gtol);CHKERRQ(info);
913:   TaoFunctionReturn(0);
914: }

918: /*@C
919:    TaoCreateMoreThuenteLineSearch - Create a line search

921:    Input Parameters:
922: +  tao - TAO_SOLVER context
923: .  fftol - the sufficient descent parameter , greater than 0.
924: -  ggtol - the curvature tolerance, greater than 0, less than 1.


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

930:    Note:
931:    This algorithm is taken from More' and Thuente, "Line search algorithms
932:    with guaranteed sufficient decrease", Argonne National Laboratory, 
933:    Technical Report MCS-P330-1092.

935:    Note:
936:    This line search enforces the strong Wolfe conditions for unconstrained
937:    optimization.  This routine is used within the following TAO unconstrained
938:    minimization solvers: Newton linesearch (tao_nls), limited memory variable 
939:    metric (tao_lmvm), and nonlinear conjugate gradient methods.

941:    Level: developer

943: .keywords: TAO_SOLVER, linesearch
944: @*/
945: int TaoCreateMoreThuenteLineSearch(TAO_SOLVER tao, double fftol, double ggtol)
946: {
947:   int info;
948:   TAO_LINESEARCH *neP;

950:   TaoFunctionBegin;

952:   info = TaoNew(TAO_LINESEARCH,&neP); CHKERRQ(info);
953:   info = PetscLogObjectMemory(tao,sizeof(TAO_LINESEARCH)); CHKERRQ(info);

955:   if (fftol<=0) fftol=0.0001;
956:   if (ggtol<=0) ggtol=0.9;

958:   neP->ftol                  = fftol;
959:   neP->rtol                  = 1.0e-10;
960:   neP->gtol                  = ggtol;
961:   neP->stepmin                  = 1.0e-20;
962:   neP->stepmax                  = 1.0e+20;
963:   neP->nfev                  = 0; 
964:   neP->bracket                  = 0; 
965:   neP->infoc              = 1;
966:   neP->maxfev                  = 30;

968:   info = TaoSetLineSearch(tao,0,
969:                           TaoSetOptions_LineSearch,
970:                           TaoApply_LineSearch,
971:                           TaoView_LineSearch,
972:                           TaoDestroy_LineSearch,
973:                           (void *) neP);CHKERRQ(info);

975:   TaoFunctionReturn(0);
976: }

980: /*@C
981:    TaoCreateMoreThuenteBoundLineSearch - Create a line search

983:    Input Parameters:
984: +  tao - TAO_SOLVER context
985: .  fftol - the sufficient descent parameter , greater than 0.
986: -  ggtol - the curvature tolerance, greater than 0, less than 1.


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

992:    Note:
993:    This algorithm is is a modification of the algorithm by More' and Thuente.
994:    The modifications concern bounds.  This algorithm steps in the direction
995:    passed into this routine.  This point get projected back into the feasible set.
996:    It tries to satisfy the Wolfe conditions, but in the context of bound constrained
997:    optimization, there may not be a point in the piecewise linear 
998:    path that satisfies the Wolfe conditions.  When the active set
999:    is changing, decrease in the objective function may be sufficient
1000:    to terminate this line search.

1002:    Note:
1003:    This routine is used within the following TAO bound constrained
1004:    minimization solvers: Newton trust region (tao_tron) and limited memory variable 
1005:    metric (tao_blmvm).

1007:    Level: developer

1009: .keywords: TAO_SOLVER, linesearch
1010: @*/
1011: int TaoCreateMoreThuenteBoundLineSearch(TAO_SOLVER tao, double fftol, double ggtol)
1012: {
1013:   int info;
1014:   TAO_LINESEARCH *neP;

1016:   TaoFunctionBegin;

1018:   info = TaoNew(TAO_LINESEARCH,&neP); CHKERRQ(info);
1019:   info = PetscLogObjectMemory(tao,sizeof(TAO_LINESEARCH)); CHKERRQ(info);

1021:   if (fftol<=0) fftol=0.0001;
1022:   if (ggtol<=0) ggtol=0.9;

1024:   neP->ftol                  = fftol;
1025:   neP->rtol                  = 1.0e-10;
1026:   neP->gtol                  = ggtol;
1027:   neP->stepmin                  = 1.0e-20;
1028:   neP->stepmax                  = 1.0e+20;
1029:   neP->nfev                  = 0; 
1030:   neP->bracket                  = 0; 
1031:   neP->infoc              = 1;
1032:   neP->maxfev                  = 30;

1034:   info = TaoSetLineSearch(tao,0,
1035:                           TaoSetOptions_LineSearch,
1036:                           TaoApply_BoundLineSearch,
1037:                           TaoView_LineSearch,
1038:                           TaoDestroy_LineSearch,
1039:                           (void *) neP);CHKERRQ(info);

1041:   TaoFunctionReturn(0);
1042: }