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: }