1: #include <petsc/private/taolinesearchimpl.h>
2: #include <../src/tao/linesearch/impls/morethuente/morethuente.h>
4: /*
5: This algorithm is taken from More' and Thuente, "Line search algorithms
6: with guaranteed sufficient decrease", Argonne National Laboratory,
7: Technical Report MCS-P330-1092.
8: */
10: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,PetscReal *stx,PetscReal *fx,PetscReal *dx,PetscReal *sty,PetscReal *fy,PetscReal *dy,PetscReal *stp,PetscReal *fp,PetscReal *dp);
14: static PetscErrorCode TaoLineSearchDestroy_MT(TaoLineSearch ls) 15: {
16: PetscErrorCode ierr;
17: TaoLineSearch_MT *mt;
21: mt = (TaoLineSearch_MT*)(ls->data);
22: if (mt->x) {
23: PetscObjectDereference((PetscObject)mt->x);
24: }
25: VecDestroy(&mt->work);
26: PetscFree(ls->data);
27: return(0);
28: }
32: static PetscErrorCode TaoLineSearchSetFromOptions_MT(PetscOptions *PetscOptionsObject,TaoLineSearch ls) 33: {
36: return(0);
37: }
41: static PetscErrorCode TaoLineSearchView_MT(TaoLineSearch ls, PetscViewer pv) 42: {
44: PetscBool isascii;
47: PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii);
48: if (isascii) {
49: PetscViewerASCIIPrintf(pv," maxf=%D, ftol=%g, gtol=%g\n",ls->max_funcs,(double)ls->rtol,(double)ls->ftol);
50: }
51: return(0);
52: }
56: /* @ TaoApply_LineSearch - This routine takes step length of 1.0.
58: Input Parameters:
59: + tao - Tao context
60: . X - current iterate (on output X contains new iterate, X + step*S)
61: . f - objective function evaluated at X
62: . G - gradient evaluated at X
63: - D - search direction
66: Info is set to 0.
68: @ */
70: static PetscErrorCode TaoLineSearchApply_MT(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s) 71: {
72: PetscErrorCode ierr;
73: TaoLineSearch_MT *mt;
75: PetscReal xtrapf = 4.0;
76: PetscReal finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
77: PetscReal dgx, dgy, dg, dg2, fx, fy, stx, sty, dgtest;
78: PetscReal ftest1=0.0, ftest2=0.0;
79: PetscInt i, stage1,n1,n2,nn1,nn2;
80: PetscReal bstepmin1, bstepmin2, bstepmax;
81: PetscBool g_computed=PETSC_FALSE; /* to prevent extra gradient computation */
90: /* comm,type,size checks are done in interface TaoLineSearchApply */
91: mt = (TaoLineSearch_MT*)(ls->data);
92: ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
94: /* Check work vector */
95: if (!mt->work) {
96: VecDuplicate(x,&mt->work);
97: mt->x = x;
98: PetscObjectReference((PetscObject)mt->x);
99: } else if (x != mt->x) {
100: VecDestroy(&mt->work);
101: VecDuplicate(x,&mt->work);
102: PetscObjectDereference((PetscObject)mt->x);
103: mt->x = x;
104: PetscObjectReference((PetscObject)mt->x);
105: }
107: if (ls->bounded) {
108: /* Compute step length needed to make all variables equal a bound */
109: /* Compute the smallest steplength that will make one nonbinding variable
110: equal the bound */
111: VecGetLocalSize(ls->upper,&n1);
112: VecGetLocalSize(mt->x, &n2);
113: VecGetSize(ls->upper,&nn1);
114: VecGetSize(mt->x,&nn2);
115: if (n1 != n2 || nn1 != nn2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Variable vector not compatible with bounds vector");
116: VecScale(s,-1.0);
117: VecBoundGradientProjection(s,x,ls->lower,ls->upper,s);
118: VecScale(s,-1.0);
119: VecStepBoundInfo(x,s,ls->lower,ls->upper,&bstepmin1,&bstepmin2,&bstepmax);
120: ls->stepmax = PetscMin(bstepmax,1.0e15);
121: }
123: VecDot(g,s,&dginit);
124: if (PetscIsInfOrNanReal(dginit)) {
125: PetscInfo1(ls,"Initial Line Search step * g is Inf or Nan (%g)\n",(double)dginit);
126: ls->reason=TAOLINESEARCH_FAILED_INFORNAN;
127: return(0);
128: }
129: if (dginit >= 0.0) {
130: PetscInfo1(ls,"Initial Line Search step * g is not descent direction (%g)\n",(double)dginit);
131: ls->reason = TAOLINESEARCH_FAILED_ASCENT;
132: return(0);
133: }
136: /* Initialization */
137: mt->bracket = 0;
138: stage1 = 1;
139: finit = *f;
140: dgtest = ls->ftol * dginit;
141: width = ls->stepmax - ls->stepmin;
142: width1 = width * 2.0;
143: VecCopy(x,mt->work);
144: /* Variable dictionary:
145: stx, fx, dgx - the step, function, and derivative at the best step
146: sty, fy, dgy - the step, function, and derivative at the other endpoint
147: of the interval of uncertainty
148: step, f, dg - the step, function, and derivative at the current step */
150: stx = 0.0;
151: fx = finit;
152: dgx = dginit;
153: sty = 0.0;
154: fy = finit;
155: dgy = dginit;
157: ls->step=ls->initstep;
158: for (i=0; i< ls->max_funcs; i++) {
159: /* Set min and max steps to correspond to the interval of uncertainty */
160: if (mt->bracket) {
161: ls->stepmin = PetscMin(stx,sty);
162: ls->stepmax = PetscMax(stx,sty);
163: } else {
164: ls->stepmin = stx;
165: ls->stepmax = ls->step + xtrapf * (ls->step - stx);
166: }
168: /* Force the step to be within the bounds */
169: ls->step = PetscMax(ls->step,ls->stepmin);
170: ls->step = PetscMin(ls->step,ls->stepmax);
172: /* If an unusual termination is to occur, then let step be the lowest
173: point obtained thus far */
174: if ((stx!=0) && (((mt->bracket) && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) ||
175: ((ls->nfeval+ls->nfgeval) >= ls->max_funcs - 1) || (mt->infoc == 0))) {
176: ls->step = stx;
177: }
179: VecCopy(x,mt->work);
180: VecAXPY(mt->work,ls->step,s); /* W = X + step*S */
182: if (ls->bounded) {
183: VecMedian(ls->lower, mt->work, ls->upper, mt->work);
184: }
185: if (ls->usegts) {
186: TaoLineSearchComputeObjectiveAndGTS(ls,mt->work,f,&dg);
187: g_computed=PETSC_FALSE;
188: } else {
189: TaoLineSearchComputeObjectiveAndGradient(ls,mt->work,f,g);
190: g_computed=PETSC_TRUE;
191: if (ls->bounded) {
192: VecDot(g,x,&dg);
193: VecDot(g,mt->work,&dg2);
194: dg = (dg2 - dg)/ls->step;
195: } else {
196: VecDot(g,s,&dg);
197: }
198: }
200: if (0 == i) {
201: ls->f_fullstep=*f;
202: }
204: if (PetscIsInfOrNanReal(*f) || PetscIsInfOrNanReal(dg)) {
205: /* User provided compute function generated Not-a-Number, assume
206: domain violation and set function value and directional
207: derivative to infinity. */
208: *f = PETSC_INFINITY;
209: dg = PETSC_INFINITY;
210: }
212: ftest1 = finit + ls->step * dgtest;
213: if (ls->bounded) {
214: ftest2 = finit + ls->step * dgtest * ls->ftol;
215: }
216: /* Convergence testing */
217: if (((*f - ftest1 <= 1.0e-10 * PetscAbsReal(finit)) && (PetscAbsReal(dg) + ls->gtol*dginit <= 0.0))) {
218: PetscInfo(ls, "Line search success: Sufficient decrease and directional deriv conditions hold\n");
219: ls->reason = TAOLINESEARCH_SUCCESS;
220: break;
221: }
223: /* Check Armijo if beyond the first breakpoint */
224: if (ls->bounded && (*f <= ftest2) && (ls->step >= bstepmin2)) {
225: PetscInfo(ls,"Line search success: Sufficient decrease.\n");
226: ls->reason = TAOLINESEARCH_SUCCESS;
227: break;
228: }
230: /* Checks for bad cases */
231: if (((mt->bracket) && (ls->step <= ls->stepmin||ls->step >= ls->stepmax)) || (!mt->infoc)) {
232: PetscInfo(ls,"Rounding errors may prevent further progress. May not be a step satisfying\n");
233: PetscInfo(ls,"sufficient decrease and curvature conditions. Tolerances may be too small.\n");
234: ls->reason = TAOLINESEARCH_HALTED_OTHER;
235: break;
236: }
237: if ((ls->step == ls->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
238: PetscInfo1(ls,"Step is at the upper bound, stepmax (%g)\n",(double)ls->stepmax);
239: ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
240: break;
241: }
242: if ((ls->step == ls->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
243: PetscInfo1(ls,"Step is at the lower bound, stepmin (%g)\n",(double)ls->stepmin);
244: ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
245: break;
246: }
247: if ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol*ls->stepmax)){
248: PetscInfo1(ls,"Relative width of interval of uncertainty is at most rtol (%g)\n",(double)ls->rtol);
249: ls->reason = TAOLINESEARCH_HALTED_RTOL;
250: break;
251: }
253: /* In the first stage, we seek a step for which the modified function
254: has a nonpositive value and nonnegative derivative */
255: if ((stage1) && (*f <= ftest1) && (dg >= dginit * PetscMin(ls->ftol, ls->gtol))) {
256: stage1 = 0;
257: }
259: /* A modified function is used to predict the step only if we
260: have not obtained a step for which the modified function has a
261: nonpositive function value and nonnegative derivative, and if a
262: lower function value has been obtained but the decrease is not
263: sufficient */
265: if ((stage1) && (*f <= fx) && (*f > ftest1)) {
266: fm = *f - ls->step * dgtest; /* Define modified function */
267: fxm = fx - stx * dgtest; /* and derivatives */
268: fym = fy - sty * dgtest;
269: dgm = dg - dgtest;
270: dgxm = dgx - dgtest;
271: dgym = dgy - dgtest;
273: /* if (dgxm * (ls->step - stx) >= 0.0) */
274: /* Update the interval of uncertainty and compute the new step */
275: Tao_mcstep(ls,&stx,&fxm,&dgxm,&sty,&fym,&dgym,&ls->step,&fm,&dgm);
277: fx = fxm + stx * dgtest; /* Reset the function and */
278: fy = fym + sty * dgtest; /* gradient values */
279: dgx = dgxm + dgtest;
280: dgy = dgym + dgtest;
281: } else {
282: /* Update the interval of uncertainty and compute the new step */
283: Tao_mcstep(ls,&stx,&fx,&dgx,&sty,&fy,&dgy,&ls->step,f,&dg);
284: }
286: /* Force a sufficient decrease in the interval of uncertainty */
287: if (mt->bracket) {
288: if (PetscAbsReal(sty - stx) >= 0.66 * width1) ls->step = stx + 0.5*(sty - stx);
289: width1 = width;
290: width = PetscAbsReal(sty - stx);
291: }
292: }
293: if ((ls->nfeval+ls->nfgeval) > ls->max_funcs) {
294: PetscInfo2(ls,"Number of line search function evals (%D) > maximum (%D)\n",(ls->nfeval+ls->nfgeval),ls->max_funcs);
295: ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
296: }
298: /* Finish computations */
299: PetscInfo2(ls,"%D function evals in line search, step = %g\n",(ls->nfeval+ls->nfgeval),(double)ls->step);
301: /* Set new solution vector and compute gradient if needed */
302: VecCopy(mt->work,x);
303: if (!g_computed) {
304: TaoLineSearchComputeGradient(ls,mt->work,g);
305: }
306: return(0);
307: }
311: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_MT(TaoLineSearch ls)312: {
313: PetscErrorCode ierr;
314: TaoLineSearch_MT *ctx;
318: PetscNewLog(ls,&ctx);
319: ctx->bracket=0;
320: ctx->infoc=1;
321: ls->data = (void*)ctx;
322: ls->initstep = 1.0;
323: ls->ops->setup=0;
324: ls->ops->reset=0;
325: ls->ops->apply=TaoLineSearchApply_MT;
326: ls->ops->view =TaoLineSearchView_MT;
327: ls->ops->destroy=TaoLineSearchDestroy_MT;
328: ls->ops->setfromoptions=TaoLineSearchSetFromOptions_MT;
329: return(0);
330: }
332: /*
333: The subroutine mcstep is taken from the work of Jorge Nocedal.
334: this is a variant of More' and Thuente's routine.
336: subroutine mcstep
338: the purpose of mcstep is to compute a safeguarded step for
339: a linesearch and to update an interval of uncertainty for
340: a minimizer of the function.
342: the parameter stx contains the step with the least function
343: value. the parameter stp contains the current step. it is
344: assumed that the derivative at stx is negative in the
345: direction of the step. if bracket is set true then a
346: minimizer has been bracketed in an interval of uncertainty
347: with endpoints stx and sty.
349: the subroutine statement is
351: subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
352: stpmin,stpmax,info)
354: where
356: stx, fx, and dx are variables which specify the step,
357: the function, and the derivative at the best step obtained
358: so far. The derivative must be negative in the direction
359: of the step, that is, dx and stp-stx must have opposite
360: signs. On output these parameters are updated appropriately.
362: sty, fy, and dy are variables which specify the step,
363: the function, and the derivative at the other endpoint of
364: the interval of uncertainty. On output these parameters are
365: updated appropriately.
367: stp, fp, and dp are variables which specify the step,
368: the function, and the derivative at the current step.
369: If bracket is set true then on input stp must be
370: between stx and sty. On output stp is set to the new step.
372: bracket is a logical variable which specifies if a minimizer
373: has been bracketed. If the minimizer has not been bracketed
374: then on input bracket must be set false. If the minimizer
375: is bracketed then on output bracket is set true.
377: stpmin and stpmax are input variables which specify lower
378: and upper bounds for the step.
380: info is an integer output variable set as follows:
381: if info = 1,2,3,4,5, then the step has been computed
382: according to one of the five cases below. otherwise
383: info = 0, and this indicates improper input parameters.
385: subprograms called
387: fortran-supplied ... abs,max,min,sqrt
389: argonne national laboratory. minpack project. june 1983
390: jorge j. more', david j. thuente
392: */
396: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,PetscReal *stx,PetscReal *fx,PetscReal *dx,PetscReal *sty,PetscReal *fy,PetscReal *dy,PetscReal *stp,PetscReal *fp,PetscReal *dp)397: {
398: TaoLineSearch_MT *mtP = (TaoLineSearch_MT *) ls->data;
399: PetscReal gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
400: PetscInt bound;
403: /* Check the input parameters for errors */
404: mtP->infoc = 0;
405: if (mtP->bracket && (*stp <= PetscMin(*stx,*sty) || (*stp >= PetscMax(*stx,*sty)))) SETERRQ(PETSC_COMM_SELF,1,"bad stp in bracket");
406: if (*dx * (*stp-*stx) >= 0.0) SETERRQ(PETSC_COMM_SELF,1,"dx * (stp-stx) >= 0.0");
407: if (ls->stepmax < ls->stepmin) SETERRQ(PETSC_COMM_SELF,1,"stepmax > stepmin");
409: /* Determine if the derivatives have opposite sign */
410: sgnd = *dp * (*dx / PetscAbsReal(*dx));
412: if (*fp > *fx) {
413: /* Case 1: a higher function value.
414: The minimum is bracketed. If the cubic step is closer
415: to stx than the quadratic step, the cubic step is taken,
416: else the average of the cubic and quadratic steps is taken. */
418: mtP->infoc = 1;
419: bound = 1;
420: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
421: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
422: s = PetscMax(s,PetscAbsReal(*dp));
423: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
424: if (*stp < *stx) gamma1 = -gamma1;
425: /* Can p be 0? Check */
426: p = (gamma1 - *dx) + theta;
427: q = ((gamma1 - *dx) + gamma1) + *dp;
428: r = p/q;
429: stpc = *stx + r*(*stp - *stx);
430: stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);
432: if (PetscAbsReal(stpc-*stx) < PetscAbsReal(stpq-*stx)) {
433: stpf = stpc;
434: } else {
435: stpf = stpc + 0.5*(stpq - stpc);
436: }
437: mtP->bracket = 1;
438: } else if (sgnd < 0.0) {
439: /* Case 2: A lower function value and derivatives of
440: opposite sign. The minimum is bracketed. If the cubic
441: step is closer to stx than the quadratic (secant) step,
442: the cubic step is taken, else the quadratic step is taken. */
444: mtP->infoc = 2;
445: bound = 0;
446: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
447: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
448: s = PetscMax(s,PetscAbsReal(*dp));
449: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
450: if (*stp > *stx) gamma1 = -gamma1;
451: p = (gamma1 - *dp) + theta;
452: q = ((gamma1 - *dp) + gamma1) + *dx;
453: r = p/q;
454: stpc = *stp + r*(*stx - *stp);
455: stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);
457: if (PetscAbsReal(stpc-*stp) > PetscAbsReal(stpq-*stp)) {
458: stpf = stpc;
459: } else {
460: stpf = stpq;
461: }
462: mtP->bracket = 1;
463: } else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
464: /* Case 3: A lower function value, derivatives of the
465: same sign, and the magnitude of the derivative decreases.
466: The cubic step is only used if the cubic tends to infinity
467: in the direction of the step or if the minimum of the cubic
468: is beyond stp. Otherwise the cubic step is defined to be
469: either stepmin or stepmax. The quadratic (secant) step is also
470: computed and if the minimum is bracketed then the step
471: closest to stx is taken, else the step farthest away is taken. */
473: mtP->infoc = 3;
474: bound = 1;
475: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
476: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
477: s = PetscMax(s,PetscAbsReal(*dp));
479: /* The case gamma1 = 0 only arises if the cubic does not tend
480: to infinity in the direction of the step. */
481: gamma1 = s*PetscSqrtScalar(PetscMax(0.0,PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s)));
482: if (*stp > *stx) gamma1 = -gamma1;
483: p = (gamma1 - *dp) + theta;
484: q = (gamma1 + (*dx - *dp)) + gamma1;
485: r = p/q;
486: if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r*(*stx - *stp);
487: else if (*stp > *stx) stpc = ls->stepmax;
488: else stpc = ls->stepmin;
489: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
491: if (mtP->bracket) {
492: if (PetscAbsReal(*stp-stpc) < PetscAbsReal(*stp-stpq)) {
493: stpf = stpc;
494: } else {
495: stpf = stpq;
496: }
497: } else {
498: if (PetscAbsReal(*stp-stpc) > PetscAbsReal(*stp-stpq)) {
499: stpf = stpc;
500: } else {
501: stpf = stpq;
502: }
503: }
504: } else {
505: /* Case 4: A lower function value, derivatives of the
506: same sign, and the magnitude of the derivative does
507: not decrease. If the minimum is not bracketed, the step
508: is either stpmin or stpmax, else the cubic step is taken. */
510: mtP->infoc = 4;
511: bound = 0;
512: if (mtP->bracket) {
513: theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
514: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dy));
515: s = PetscMax(s,PetscAbsReal(*dp));
516: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dy/s)*(*dp/s));
517: if (*stp > *sty) gamma1 = -gamma1;
518: p = (gamma1 - *dp) + theta;
519: q = ((gamma1 - *dp) + gamma1) + *dy;
520: r = p/q;
521: stpc = *stp + r*(*sty - *stp);
522: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
524: stpf = stpc;
525: } else if (*stp > *stx) {
526: stpf = ls->stepmax;
527: } else {
528: stpf = ls->stepmin;
529: }
530: }
532: /* Update the interval of uncertainty. This update does not
533: depend on the new step or the case analysis above. */
535: if (*fp > *fx) {
536: *sty = *stp;
537: *fy = *fp;
538: *dy = *dp;
539: } else {
540: if (sgnd < 0.0) {
541: *sty = *stx;
542: *fy = *fx;
543: *dy = *dx;
544: }
545: *stx = *stp;
546: *fx = *fp;
547: *dx = *dp;
548: }
550: /* Compute the new step and safeguard it. */
551: stpf = PetscMin(ls->stepmax,stpf);
552: stpf = PetscMax(ls->stepmin,stpf);
553: *stp = stpf;
554: if (mtP->bracket && bound) {
555: if (*sty > *stx) {
556: *stp = PetscMin(*stx+0.66*(*sty-*stx),*stp);
557: } else {
558: *stp = PetscMax(*stx+0.66*(*sty-*stx),*stp);
559: }
560: }
561: return(0);
562: }