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(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: }
309: EXTERN_C_BEGIN
312: PetscErrorCode TaoLineSearchCreate_MT(TaoLineSearch ls)313: {
314: PetscErrorCode ierr;
315: TaoLineSearch_MT *ctx;
319: PetscNewLog(ls,&ctx);
320: ctx->bracket=0;
321: ctx->infoc=1;
322: ls->data = (void*)ctx;
323: ls->initstep = 1.0;
324: ls->ops->setup=0;
325: ls->ops->reset=0;
326: ls->ops->apply=TaoLineSearchApply_MT;
327: ls->ops->view =TaoLineSearchView_MT;
328: ls->ops->destroy=TaoLineSearchDestroy_MT;
329: ls->ops->setfromoptions=TaoLineSearchSetFromOptions_MT;
330: return(0);
331: }
332: EXTERN_C_END
334: /*
335: The subroutine mcstep is taken from the work of Jorge Nocedal.
336: this is a variant of More' and Thuente's routine.
338: subroutine mcstep
340: the purpose of mcstep is to compute a safeguarded step for
341: a linesearch and to update an interval of uncertainty for
342: a minimizer of the function.
344: the parameter stx contains the step with the least function
345: value. the parameter stp contains the current step. it is
346: assumed that the derivative at stx is negative in the
347: direction of the step. if bracket is set true then a
348: minimizer has been bracketed in an interval of uncertainty
349: with endpoints stx and sty.
351: the subroutine statement is
353: subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
354: stpmin,stpmax,info)
356: where
358: stx, fx, and dx are variables which specify the step,
359: the function, and the derivative at the best step obtained
360: so far. The derivative must be negative in the direction
361: of the step, that is, dx and stp-stx must have opposite
362: signs. On output these parameters are updated appropriately.
364: sty, fy, and dy are variables which specify the step,
365: the function, and the derivative at the other endpoint of
366: the interval of uncertainty. On output these parameters are
367: updated appropriately.
369: stp, fp, and dp are variables which specify the step,
370: the function, and the derivative at the current step.
371: If bracket is set true then on input stp must be
372: between stx and sty. On output stp is set to the new step.
374: bracket is a logical variable which specifies if a minimizer
375: has been bracketed. If the minimizer has not been bracketed
376: then on input bracket must be set false. If the minimizer
377: is bracketed then on output bracket is set true.
379: stpmin and stpmax are input variables which specify lower
380: and upper bounds for the step.
382: info is an integer output variable set as follows:
383: if info = 1,2,3,4,5, then the step has been computed
384: according to one of the five cases below. otherwise
385: info = 0, and this indicates improper input parameters.
387: subprograms called
389: fortran-supplied ... abs,max,min,sqrt
391: argonne national laboratory. minpack project. june 1983
392: jorge j. more', david j. thuente
394: */
398: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,PetscReal *stx,PetscReal *fx,PetscReal *dx,PetscReal *sty,PetscReal *fy,PetscReal *dy,PetscReal *stp,PetscReal *fp,PetscReal *dp)399: {
400: TaoLineSearch_MT *mtP = (TaoLineSearch_MT *) ls->data;
401: PetscReal gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
402: PetscInt bound;
405: /* Check the input parameters for errors */
406: mtP->infoc = 0;
407: if (mtP->bracket && (*stp <= PetscMin(*stx,*sty) || (*stp >= PetscMax(*stx,*sty)))) SETERRQ(PETSC_COMM_SELF,1,"bad stp in bracket");
408: if (*dx * (*stp-*stx) >= 0.0) SETERRQ(PETSC_COMM_SELF,1,"dx * (stp-stx) >= 0.0");
409: if (ls->stepmax < ls->stepmin) SETERRQ(PETSC_COMM_SELF,1,"stepmax > stepmin");
411: /* Determine if the derivatives have opposite sign */
412: sgnd = *dp * (*dx / PetscAbsReal(*dx));
414: if (*fp > *fx) {
415: /* Case 1: a higher function value.
416: The minimum is bracketed. If the cubic step is closer
417: to stx than the quadratic step, the cubic step is taken,
418: else the average of the cubic and quadratic steps is taken. */
420: mtP->infoc = 1;
421: bound = 1;
422: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
423: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
424: s = PetscMax(s,PetscAbsReal(*dp));
425: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
426: if (*stp < *stx) gamma1 = -gamma1;
427: /* Can p be 0? Check */
428: p = (gamma1 - *dx) + theta;
429: q = ((gamma1 - *dx) + gamma1) + *dp;
430: r = p/q;
431: stpc = *stx + r*(*stp - *stx);
432: stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);
434: if (PetscAbsReal(stpc-*stx) < PetscAbsReal(stpq-*stx)) {
435: stpf = stpc;
436: } else {
437: stpf = stpc + 0.5*(stpq - stpc);
438: }
439: mtP->bracket = 1;
440: } else if (sgnd < 0.0) {
441: /* Case 2: A lower function value and derivatives of
442: opposite sign. The minimum is bracketed. If the cubic
443: step is closer to stx than the quadratic (secant) step,
444: the cubic step is taken, else the quadratic step is taken. */
446: mtP->infoc = 2;
447: bound = 0;
448: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
449: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
450: s = PetscMax(s,PetscAbsReal(*dp));
451: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
452: if (*stp > *stx) gamma1 = -gamma1;
453: p = (gamma1 - *dp) + theta;
454: q = ((gamma1 - *dp) + gamma1) + *dx;
455: r = p/q;
456: stpc = *stp + r*(*stx - *stp);
457: stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);
459: if (PetscAbsReal(stpc-*stp) > PetscAbsReal(stpq-*stp)) {
460: stpf = stpc;
461: } else {
462: stpf = stpq;
463: }
464: mtP->bracket = 1;
465: } else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
466: /* Case 3: A lower function value, derivatives of the
467: same sign, and the magnitude of the derivative decreases.
468: The cubic step is only used if the cubic tends to infinity
469: in the direction of the step or if the minimum of the cubic
470: is beyond stp. Otherwise the cubic step is defined to be
471: either stepmin or stepmax. The quadratic (secant) step is also
472: computed and if the minimum is bracketed then the step
473: closest to stx is taken, else the step farthest away is taken. */
475: mtP->infoc = 3;
476: bound = 1;
477: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
478: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
479: s = PetscMax(s,PetscAbsReal(*dp));
481: /* The case gamma1 = 0 only arises if the cubic does not tend
482: to infinity in the direction of the step. */
483: gamma1 = s*PetscSqrtScalar(PetscMax(0.0,PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s)));
484: if (*stp > *stx) gamma1 = -gamma1;
485: p = (gamma1 - *dp) + theta;
486: q = (gamma1 + (*dx - *dp)) + gamma1;
487: r = p/q;
488: if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r*(*stx - *stp);
489: else if (*stp > *stx) stpc = ls->stepmax;
490: else stpc = ls->stepmin;
491: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
493: if (mtP->bracket) {
494: if (PetscAbsReal(*stp-stpc) < PetscAbsReal(*stp-stpq)) {
495: stpf = stpc;
496: } else {
497: stpf = stpq;
498: }
499: } else {
500: if (PetscAbsReal(*stp-stpc) > PetscAbsReal(*stp-stpq)) {
501: stpf = stpc;
502: } else {
503: stpf = stpq;
504: }
505: }
506: } else {
507: /* Case 4: A lower function value, derivatives of the
508: same sign, and the magnitude of the derivative does
509: not decrease. If the minimum is not bracketed, the step
510: is either stpmin or stpmax, else the cubic step is taken. */
512: mtP->infoc = 4;
513: bound = 0;
514: if (mtP->bracket) {
515: theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
516: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dy));
517: s = PetscMax(s,PetscAbsReal(*dp));
518: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dy/s)*(*dp/s));
519: if (*stp > *sty) gamma1 = -gamma1;
520: p = (gamma1 - *dp) + theta;
521: q = ((gamma1 - *dp) + gamma1) + *dy;
522: r = p/q;
523: stpc = *stp + r*(*sty - *stp);
524: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
526: stpf = stpc;
527: } else if (*stp > *stx) {
528: stpf = ls->stepmax;
529: } else {
530: stpf = ls->stepmin;
531: }
532: }
534: /* Update the interval of uncertainty. This update does not
535: depend on the new step or the case analysis above. */
537: if (*fp > *fx) {
538: *sty = *stp;
539: *fy = *fp;
540: *dy = *dp;
541: } else {
542: if (sgnd < 0.0) {
543: *sty = *stx;
544: *fy = *fx;
545: *dy = *dx;
546: }
547: *stx = *stp;
548: *fx = *fp;
549: *dx = *dp;
550: }
552: /* Compute the new step and safeguard it. */
553: stpf = PetscMin(ls->stepmax,stpf);
554: stpf = PetscMax(ls->stepmin,stpf);
555: *stp = stpf;
556: if (mtP->bracket && bound) {
557: if (*sty > *stx) {
558: *stp = PetscMin(*stx+0.66*(*sty-*stx),*stp);
559: } else {
560: *stp = PetscMax(*stx+0.66*(*sty-*stx),*stp);
561: }
562: }
563: return(0);
564: }