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);
12: static PetscErrorCode TaoLineSearchDestroy_MT(TaoLineSearch ls) 13: {
14: PetscErrorCode ierr;
15: TaoLineSearch_MT *mt;
19: mt = (TaoLineSearch_MT*)(ls->data);
20: if (mt->x) {
21: PetscObjectDereference((PetscObject)mt->x);
22: }
23: VecDestroy(&mt->work);
24: PetscFree(ls->data);
25: return(0);
26: }
28: static PetscErrorCode TaoLineSearchSetFromOptions_MT(PetscOptionItems *PetscOptionsObject,TaoLineSearch ls) 29: {
32: return(0);
33: }
35: static PetscErrorCode TaoLineSearchMonitor_MT(TaoLineSearch ls) 36: {
37: TaoLineSearch_MT *mt = (TaoLineSearch_MT*)ls->data;
38: PetscErrorCode ierr;
39: 41: PetscViewerASCIIPrintf(ls->viewer, "stx: %g, fx: %g, dgx: %g\n", (double)mt->stx, (double)mt->fx, (double)mt->dgx);
42: PetscViewerASCIIPrintf(ls->viewer, "sty: %g, fy: %g, dgy: %g\n", (double)mt->sty, (double)mt->fy, (double)mt->dgy);
43: return(0);
44: }
46: static PetscErrorCode TaoLineSearchApply_MT(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s) 47: {
48: PetscErrorCode ierr;
49: TaoLineSearch_MT *mt;
51: PetscReal xtrapf = 4.0;
52: PetscReal finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
53: PetscReal dgx, dgy, dg, dg2, fx, fy, stx, sty, dgtest;
54: PetscReal ftest1=0.0, ftest2=0.0;
55: PetscInt i, stage1,n1,n2,nn1,nn2;
56: PetscReal bstepmin1, bstepmin2, bstepmax;
57: PetscBool g_computed=PETSC_FALSE; /* to prevent extra gradient computation */
65: 66: TaoLineSearchMonitor(ls, 0, *f, 0.0);
68: /* comm,type,size checks are done in interface TaoLineSearchApply */
69: mt = (TaoLineSearch_MT*)(ls->data);
70: ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
72: /* Check work vector */
73: if (!mt->work) {
74: VecDuplicate(x,&mt->work);
75: mt->x = x;
76: PetscObjectReference((PetscObject)mt->x);
77: } else if (x != mt->x) {
78: VecDestroy(&mt->work);
79: VecDuplicate(x,&mt->work);
80: PetscObjectDereference((PetscObject)mt->x);
81: mt->x = x;
82: PetscObjectReference((PetscObject)mt->x);
83: }
85: if (ls->bounded) {
86: /* Compute step length needed to make all variables equal a bound */
87: /* Compute the smallest steplength that will make one nonbinding variable
88: equal the bound */
89: VecGetLocalSize(ls->upper,&n1);
90: VecGetLocalSize(mt->x, &n2);
91: VecGetSize(ls->upper,&nn1);
92: VecGetSize(mt->x,&nn2);
93: if (n1 != n2 || nn1 != nn2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Variable vector not compatible with bounds vector");
94: VecScale(s,-1.0);
95: VecBoundGradientProjection(s,x,ls->lower,ls->upper,s);
96: VecScale(s,-1.0);
97: VecStepBoundInfo(x,s,ls->lower,ls->upper,&bstepmin1,&bstepmin2,&bstepmax);
98: ls->stepmax = PetscMin(bstepmax,1.0e15);
99: }
101: VecDot(g,s,&dginit);
102: if (PetscIsInfOrNanReal(dginit)) {
103: PetscInfo1(ls,"Initial Line Search step * g is Inf or Nan (%g)\n",(double)dginit);
104: ls->reason=TAOLINESEARCH_FAILED_INFORNAN;
105: return(0);
106: }
107: if (dginit >= 0.0) {
108: PetscInfo1(ls,"Initial Line Search step * g is not descent direction (%g)\n",(double)dginit);
109: ls->reason = TAOLINESEARCH_FAILED_ASCENT;
110: return(0);
111: }
113: /* Initialization */
114: mt->bracket = 0;
115: stage1 = 1;
116: finit = *f;
117: dgtest = ls->ftol * dginit;
118: width = ls->stepmax - ls->stepmin;
119: width1 = width * 2.0;
120: VecCopy(x,mt->work);
121: /* Variable dictionary:
122: stx, fx, dgx - the step, function, and derivative at the best step
123: sty, fy, dgy - the step, function, and derivative at the other endpoint
124: of the interval of uncertainty
125: step, f, dg - the step, function, and derivative at the current step */
127: stx = 0.0;
128: fx = finit;
129: dgx = dginit;
130: sty = 0.0;
131: fy = finit;
132: dgy = dginit;
134: ls->step=ls->initstep;
135: for (i=0; i< ls->max_funcs; i++) {
136: /* Set min and max steps to correspond to the interval of uncertainty */
137: if (mt->bracket) {
138: ls->stepmin = PetscMin(stx,sty);
139: ls->stepmax = PetscMax(stx,sty);
140: } else {
141: ls->stepmin = stx;
142: ls->stepmax = ls->step + xtrapf * (ls->step - stx);
143: }
145: /* Force the step to be within the bounds */
146: ls->step = PetscMax(ls->step,ls->stepmin);
147: ls->step = PetscMin(ls->step,ls->stepmax);
149: /* If an unusual termination is to occur, then let step be the lowest
150: point obtained thus far */
151: if ((stx!=0) && (((mt->bracket) && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) ||
152: ((ls->nfeval+ls->nfgeval) >= ls->max_funcs - 1) || (mt->infoc == 0))) {
153: ls->step = stx;
154: }
156: VecCopy(x,mt->work);
157: VecAXPY(mt->work,ls->step,s); /* W = X + step*S */
159: if (ls->bounded) {
160: VecMedian(ls->lower, mt->work, ls->upper, mt->work);
161: }
162: if (ls->usegts) {
163: TaoLineSearchComputeObjectiveAndGTS(ls,mt->work,f,&dg);
164: g_computed=PETSC_FALSE;
165: } else {
166: TaoLineSearchComputeObjectiveAndGradient(ls,mt->work,f,g);
167: g_computed=PETSC_TRUE;
168: if (ls->bounded) {
169: VecDot(g,x,&dg);
170: VecDot(g,mt->work,&dg2);
171: dg = (dg2 - dg)/ls->step;
172: } else {
173: VecDot(g,s,&dg);
174: }
175: }
176: 177: /* update bracketing parameters in the MT context for printouts in monitor */
178: mt->stx = stx;
179: mt->fx = fx;
180: mt->dgx = dgx;
181: mt->sty = sty;
182: mt->fy = fy;
183: mt->dgy = dgy;
184: TaoLineSearchMonitor(ls, i+1, *f, ls->step);
186: if (0 == i) {
187: ls->f_fullstep=*f;
188: }
190: if (PetscIsInfOrNanReal(*f) || PetscIsInfOrNanReal(dg)) {
191: /* User provided compute function generated Not-a-Number, assume
192: domain violation and set function value and directional
193: derivative to infinity. */
194: *f = PETSC_INFINITY;
195: dg = PETSC_INFINITY;
196: }
198: ftest1 = finit + ls->step * dgtest;
199: if (ls->bounded) {
200: ftest2 = finit + ls->step * dgtest * ls->ftol;
201: }
202: /* Convergence testing */
203: if (((*f - ftest1 <= 1.0e-10 * PetscAbsReal(finit)) && (PetscAbsReal(dg) + ls->gtol*dginit <= 0.0))) {
204: PetscInfo(ls, "Line search success: Sufficient decrease and directional deriv conditions hold\n");
205: ls->reason = TAOLINESEARCH_SUCCESS;
206: break;
207: }
209: /* Check Armijo if beyond the first breakpoint */
210: if (ls->bounded && (*f <= ftest2) && (ls->step >= bstepmin2)) {
211: PetscInfo(ls,"Line search success: Sufficient decrease.\n");
212: ls->reason = TAOLINESEARCH_SUCCESS;
213: break;
214: }
216: /* Checks for bad cases */
217: if (((mt->bracket) && (ls->step <= ls->stepmin||ls->step >= ls->stepmax)) || (!mt->infoc)) {
218: PetscInfo(ls,"Rounding errors may prevent further progress. May not be a step satisfying\n");
219: PetscInfo(ls,"sufficient decrease and curvature conditions. Tolerances may be too small.\n");
220: ls->reason = TAOLINESEARCH_HALTED_OTHER;
221: break;
222: }
223: if ((ls->step == ls->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
224: PetscInfo1(ls,"Step is at the upper bound, stepmax (%g)\n",(double)ls->stepmax);
225: ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
226: break;
227: }
228: if ((ls->step == ls->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
229: PetscInfo1(ls,"Step is at the lower bound, stepmin (%g)\n",(double)ls->stepmin);
230: ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
231: break;
232: }
233: if ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol*ls->stepmax)){
234: PetscInfo1(ls,"Relative width of interval of uncertainty is at most rtol (%g)\n",(double)ls->rtol);
235: ls->reason = TAOLINESEARCH_HALTED_RTOL;
236: break;
237: }
239: /* In the first stage, we seek a step for which the modified function
240: has a nonpositive value and nonnegative derivative */
241: if ((stage1) && (*f <= ftest1) && (dg >= dginit * PetscMin(ls->ftol, ls->gtol))) {
242: stage1 = 0;
243: }
245: /* A modified function is used to predict the step only if we
246: have not obtained a step for which the modified function has a
247: nonpositive function value and nonnegative derivative, and if a
248: lower function value has been obtained but the decrease is not
249: sufficient */
251: if ((stage1) && (*f <= fx) && (*f > ftest1)) {
252: fm = *f - ls->step * dgtest; /* Define modified function */
253: fxm = fx - stx * dgtest; /* and derivatives */
254: fym = fy - sty * dgtest;
255: dgm = dg - dgtest;
256: dgxm = dgx - dgtest;
257: dgym = dgy - dgtest;
259: /* if (dgxm * (ls->step - stx) >= 0.0) */
260: /* Update the interval of uncertainty and compute the new step */
261: Tao_mcstep(ls,&stx,&fxm,&dgxm,&sty,&fym,&dgym,&ls->step,&fm,&dgm);
263: fx = fxm + stx * dgtest; /* Reset the function and */
264: fy = fym + sty * dgtest; /* gradient values */
265: dgx = dgxm + dgtest;
266: dgy = dgym + dgtest;
267: } else {
268: /* Update the interval of uncertainty and compute the new step */
269: Tao_mcstep(ls,&stx,&fx,&dgx,&sty,&fy,&dgy,&ls->step,f,&dg);
270: }
272: /* Force a sufficient decrease in the interval of uncertainty */
273: if (mt->bracket) {
274: if (PetscAbsReal(sty - stx) >= 0.66 * width1) ls->step = stx + 0.5*(sty - stx);
275: width1 = width;
276: width = PetscAbsReal(sty - stx);
277: }
278: }
279: if ((ls->nfeval+ls->nfgeval) > ls->max_funcs) {
280: PetscInfo2(ls,"Number of line search function evals (%D) > maximum (%D)\n",(ls->nfeval+ls->nfgeval),ls->max_funcs);
281: ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
282: }
284: /* Finish computations */
285: PetscInfo2(ls,"%D function evals in line search, step = %g\n",(ls->nfeval+ls->nfgeval),(double)ls->step);
287: /* Set new solution vector and compute gradient if needed */
288: VecCopy(mt->work,x);
289: if (!g_computed) {
290: TaoLineSearchComputeGradient(ls,mt->work,g);
291: }
292: return(0);
293: }
295: /*MC
296: TAOLINESEARCHMT - Line-search type with cubic interpolation that satisfies both the sufficient decrease and
297: curvature conditions. This method can take step lengths greater than 1.
299: More-Thuente line-search can be selected with "-tao_ls_type more-thuente".
301: References:
302: . 1. - JORGE J. MORE AND DAVID J. THUENTE, LINE SEARCH ALGORITHMS WITH GUARANTEED SUFFICIENT DECREASE.
303: ACM Trans. Math. Software 20, no. 3 (1994): 286-307.
305: Level: developer
307: .seealso: TaoLineSearchCreate(), TaoLineSearchSetType(), TaoLineSearchApply()
309: .keywords: Tao, linesearch
310: M*/
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->destroy=TaoLineSearchDestroy_MT;
327: ls->ops->setfromoptions=TaoLineSearchSetFromOptions_MT;
328: ls->ops->monitor=TaoLineSearchMonitor_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: */
394: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,PetscReal *stx,PetscReal *fx,PetscReal *dx,PetscReal *sty,PetscReal *fy,PetscReal *dy,PetscReal *stp,PetscReal *fp,PetscReal *dp)395: {
396: TaoLineSearch_MT *mtP = (TaoLineSearch_MT *) ls->data;
397: PetscReal gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
398: PetscInt bound;
401: /* Check the input parameters for errors */
402: mtP->infoc = 0;
403: if (mtP->bracket && (*stp <= PetscMin(*stx,*sty) || (*stp >= PetscMax(*stx,*sty)))) SETERRQ(PETSC_COMM_SELF,1,"bad stp in bracket");
404: if (*dx * (*stp-*stx) >= 0.0) SETERRQ(PETSC_COMM_SELF,1,"dx * (stp-stx) >= 0.0");
405: if (ls->stepmax < ls->stepmin) SETERRQ(PETSC_COMM_SELF,1,"stepmax > stepmin");
407: /* Determine if the derivatives have opposite sign */
408: sgnd = *dp * (*dx / PetscAbsReal(*dx));
410: if (*fp > *fx) {
411: /* Case 1: a higher function value.
412: The minimum is bracketed. If the cubic step is closer
413: to stx than the quadratic step, the cubic step is taken,
414: else the average of the cubic and quadratic steps is taken. */
416: mtP->infoc = 1;
417: bound = 1;
418: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
419: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
420: s = PetscMax(s,PetscAbsReal(*dp));
421: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
422: if (*stp < *stx) gamma1 = -gamma1;
423: /* Can p be 0? Check */
424: p = (gamma1 - *dx) + theta;
425: q = ((gamma1 - *dx) + gamma1) + *dp;
426: r = p/q;
427: stpc = *stx + r*(*stp - *stx);
428: stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);
430: if (PetscAbsReal(stpc-*stx) < PetscAbsReal(stpq-*stx)) {
431: stpf = stpc;
432: } else {
433: stpf = stpc + 0.5*(stpq - stpc);
434: }
435: mtP->bracket = 1;
436: } else if (sgnd < 0.0) {
437: /* Case 2: A lower function value and derivatives of
438: opposite sign. The minimum is bracketed. If the cubic
439: step is closer to stx than the quadratic (secant) step,
440: the cubic step is taken, else the quadratic step is taken. */
442: mtP->infoc = 2;
443: bound = 0;
444: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
445: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
446: s = PetscMax(s,PetscAbsReal(*dp));
447: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
448: if (*stp > *stx) gamma1 = -gamma1;
449: p = (gamma1 - *dp) + theta;
450: q = ((gamma1 - *dp) + gamma1) + *dx;
451: r = p/q;
452: stpc = *stp + r*(*stx - *stp);
453: stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);
455: if (PetscAbsReal(stpc-*stp) > PetscAbsReal(stpq-*stp)) {
456: stpf = stpc;
457: } else {
458: stpf = stpq;
459: }
460: mtP->bracket = 1;
461: } else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
462: /* Case 3: A lower function value, derivatives of the
463: same sign, and the magnitude of the derivative decreases.
464: The cubic step is only used if the cubic tends to infinity
465: in the direction of the step or if the minimum of the cubic
466: is beyond stp. Otherwise the cubic step is defined to be
467: either stepmin or stepmax. The quadratic (secant) step is also
468: computed and if the minimum is bracketed then the step
469: closest to stx is taken, else the step farthest away is taken. */
471: mtP->infoc = 3;
472: bound = 1;
473: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
474: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
475: s = PetscMax(s,PetscAbsReal(*dp));
477: /* The case gamma1 = 0 only arises if the cubic does not tend
478: to infinity in the direction of the step. */
479: gamma1 = s*PetscSqrtScalar(PetscMax(0.0,PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s)));
480: if (*stp > *stx) gamma1 = -gamma1;
481: p = (gamma1 - *dp) + theta;
482: q = (gamma1 + (*dx - *dp)) + gamma1;
483: r = p/q;
484: if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r*(*stx - *stp);
485: else if (*stp > *stx) stpc = ls->stepmax;
486: else stpc = ls->stepmin;
487: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
489: if (mtP->bracket) {
490: if (PetscAbsReal(*stp-stpc) < PetscAbsReal(*stp-stpq)) {
491: stpf = stpc;
492: } else {
493: stpf = stpq;
494: }
495: } else {
496: if (PetscAbsReal(*stp-stpc) > PetscAbsReal(*stp-stpq)) {
497: stpf = stpc;
498: } else {
499: stpf = stpq;
500: }
501: }
502: } else {
503: /* Case 4: A lower function value, derivatives of the
504: same sign, and the magnitude of the derivative does
505: not decrease. If the minimum is not bracketed, the step
506: is either stpmin or stpmax, else the cubic step is taken. */
508: mtP->infoc = 4;
509: bound = 0;
510: if (mtP->bracket) {
511: theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
512: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dy));
513: s = PetscMax(s,PetscAbsReal(*dp));
514: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dy/s)*(*dp/s));
515: if (*stp > *sty) gamma1 = -gamma1;
516: p = (gamma1 - *dp) + theta;
517: q = ((gamma1 - *dp) + gamma1) + *dy;
518: r = p/q;
519: stpc = *stp + r*(*sty - *stp);
520: stpf = stpc;
521: } else if (*stp > *stx) {
522: stpf = ls->stepmax;
523: } else {
524: stpf = ls->stepmin;
525: }
526: }
528: /* Update the interval of uncertainty. This update does not
529: depend on the new step or the case analysis above. */
531: if (*fp > *fx) {
532: *sty = *stp;
533: *fy = *fp;
534: *dy = *dp;
535: } else {
536: if (sgnd < 0.0) {
537: *sty = *stx;
538: *fy = *fx;
539: *dy = *dx;
540: }
541: *stx = *stp;
542: *fx = *fp;
543: *dx = *dp;
544: }
546: /* Compute the new step and safeguard it. */
547: stpf = PetscMin(ls->stepmax,stpf);
548: stpf = PetscMax(ls->stepmin,stpf);
549: *stp = stpf;
550: if (mtP->bracket && bound) {
551: if (*sty > *stx) {
552: *stp = PetscMin(*stx+0.66*(*sty-*stx),*stp);
553: } else {
554: *stp = PetscMax(*stx+0.66*(*sty-*stx),*stp);
555: }
556: }
557: return(0);
558: }