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: /* @ TaoApply_LineSearch - This routine takes step length of 1.0.
48: Input Parameters:
49: + tao - Tao context
50: . X - current iterate (on output X contains new iterate, X + step*S)
51: . f - objective function evaluated at X
52: . G - gradient evaluated at X
53: - D - search direction
56: Info is set to 0.
58: @ */
60: static PetscErrorCode TaoLineSearchApply_MT(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s) 61: {
62: PetscErrorCode ierr;
63: TaoLineSearch_MT *mt;
65: PetscReal xtrapf = 4.0;
66: PetscReal finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
67: PetscReal dgx, dgy, dg, dg2, fx, fy, stx, sty, dgtest;
68: PetscReal ftest1=0.0, ftest2=0.0;
69: PetscInt i, stage1,n1,n2,nn1,nn2;
70: PetscReal bstepmin1, bstepmin2, bstepmax;
71: PetscBool g_computed=PETSC_FALSE; /* to prevent extra gradient computation */
79: 80: TaoLineSearchMonitor(ls, 0, *f, 0.0);
82: /* comm,type,size checks are done in interface TaoLineSearchApply */
83: mt = (TaoLineSearch_MT*)(ls->data);
84: ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
86: /* Check work vector */
87: if (!mt->work) {
88: VecDuplicate(x,&mt->work);
89: mt->x = x;
90: PetscObjectReference((PetscObject)mt->x);
91: } else if (x != mt->x) {
92: VecDestroy(&mt->work);
93: VecDuplicate(x,&mt->work);
94: PetscObjectDereference((PetscObject)mt->x);
95: mt->x = x;
96: PetscObjectReference((PetscObject)mt->x);
97: }
99: if (ls->bounded) {
100: /* Compute step length needed to make all variables equal a bound */
101: /* Compute the smallest steplength that will make one nonbinding variable
102: equal the bound */
103: VecGetLocalSize(ls->upper,&n1);
104: VecGetLocalSize(mt->x, &n2);
105: VecGetSize(ls->upper,&nn1);
106: VecGetSize(mt->x,&nn2);
107: if (n1 != n2 || nn1 != nn2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Variable vector not compatible with bounds vector");
108: VecScale(s,-1.0);
109: VecBoundGradientProjection(s,x,ls->lower,ls->upper,s);
110: VecScale(s,-1.0);
111: VecStepBoundInfo(x,s,ls->lower,ls->upper,&bstepmin1,&bstepmin2,&bstepmax);
112: ls->stepmax = PetscMin(bstepmax,1.0e15);
113: }
115: VecDot(g,s,&dginit);
116: if (PetscIsInfOrNanReal(dginit)) {
117: PetscInfo1(ls,"Initial Line Search step * g is Inf or Nan (%g)\n",(double)dginit);
118: ls->reason=TAOLINESEARCH_FAILED_INFORNAN;
119: return(0);
120: }
121: if (dginit >= 0.0) {
122: PetscInfo1(ls,"Initial Line Search step * g is not descent direction (%g)\n",(double)dginit);
123: ls->reason = TAOLINESEARCH_FAILED_ASCENT;
124: return(0);
125: }
127: /* Initialization */
128: mt->bracket = 0;
129: stage1 = 1;
130: finit = *f;
131: dgtest = ls->ftol * dginit;
132: width = ls->stepmax - ls->stepmin;
133: width1 = width * 2.0;
134: VecCopy(x,mt->work);
135: /* Variable dictionary:
136: stx, fx, dgx - the step, function, and derivative at the best step
137: sty, fy, dgy - the step, function, and derivative at the other endpoint
138: of the interval of uncertainty
139: step, f, dg - the step, function, and derivative at the current step */
141: stx = 0.0;
142: fx = finit;
143: dgx = dginit;
144: sty = 0.0;
145: fy = finit;
146: dgy = dginit;
148: ls->step=ls->initstep;
149: for (i=0; i< ls->max_funcs; i++) {
150: /* Set min and max steps to correspond to the interval of uncertainty */
151: if (mt->bracket) {
152: ls->stepmin = PetscMin(stx,sty);
153: ls->stepmax = PetscMax(stx,sty);
154: } else {
155: ls->stepmin = stx;
156: ls->stepmax = ls->step + xtrapf * (ls->step - stx);
157: }
159: /* Force the step to be within the bounds */
160: ls->step = PetscMax(ls->step,ls->stepmin);
161: ls->step = PetscMin(ls->step,ls->stepmax);
163: /* If an unusual termination is to occur, then let step be the lowest
164: point obtained thus far */
165: if ((stx!=0) && (((mt->bracket) && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) ||
166: ((ls->nfeval+ls->nfgeval) >= ls->max_funcs - 1) || (mt->infoc == 0))) {
167: ls->step = stx;
168: }
170: VecCopy(x,mt->work);
171: VecAXPY(mt->work,ls->step,s); /* W = X + step*S */
173: if (ls->bounded) {
174: VecMedian(ls->lower, mt->work, ls->upper, mt->work);
175: }
176: if (ls->usegts) {
177: TaoLineSearchComputeObjectiveAndGTS(ls,mt->work,f,&dg);
178: g_computed=PETSC_FALSE;
179: } else {
180: TaoLineSearchComputeObjectiveAndGradient(ls,mt->work,f,g);
181: g_computed=PETSC_TRUE;
182: if (ls->bounded) {
183: VecDot(g,x,&dg);
184: VecDot(g,mt->work,&dg2);
185: dg = (dg2 - dg)/ls->step;
186: } else {
187: VecDot(g,s,&dg);
188: }
189: }
190: 191: /* update bracketing parameters in the MT context for printouts in monitor */
192: mt->stx = stx;
193: mt->fx = fx;
194: mt->dgx = dgx;
195: mt->sty = sty;
196: mt->fy = fy;
197: mt->dgy = dgy;
198: TaoLineSearchMonitor(ls, i+1, *f, ls->step);
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: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_MT(TaoLineSearch ls)310: {
311: PetscErrorCode ierr;
312: TaoLineSearch_MT *ctx;
316: PetscNewLog(ls,&ctx);
317: ctx->bracket=0;
318: ctx->infoc=1;
319: ls->data = (void*)ctx;
320: ls->initstep = 1.0;
321: ls->ops->setup=0;
322: ls->ops->reset=0;
323: ls->ops->apply=TaoLineSearchApply_MT;
324: ls->ops->destroy=TaoLineSearchDestroy_MT;
325: ls->ops->setfromoptions=TaoLineSearchSetFromOptions_MT;
326: ls->ops->monitor=TaoLineSearchMonitor_MT;
327: return(0);
328: }
330: /*
331: The subroutine mcstep is taken from the work of Jorge Nocedal.
332: this is a variant of More' and Thuente's routine.
334: subroutine mcstep
336: the purpose of mcstep is to compute a safeguarded step for
337: a linesearch and to update an interval of uncertainty for
338: a minimizer of the function.
340: the parameter stx contains the step with the least function
341: value. the parameter stp contains the current step. it is
342: assumed that the derivative at stx is negative in the
343: direction of the step. if bracket is set true then a
344: minimizer has been bracketed in an interval of uncertainty
345: with endpoints stx and sty.
347: the subroutine statement is
349: subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
350: stpmin,stpmax,info)
352: where
354: stx, fx, and dx are variables which specify the step,
355: the function, and the derivative at the best step obtained
356: so far. The derivative must be negative in the direction
357: of the step, that is, dx and stp-stx must have opposite
358: signs. On output these parameters are updated appropriately.
360: sty, fy, and dy are variables which specify the step,
361: the function, and the derivative at the other endpoint of
362: the interval of uncertainty. On output these parameters are
363: updated appropriately.
365: stp, fp, and dp are variables which specify the step,
366: the function, and the derivative at the current step.
367: If bracket is set true then on input stp must be
368: between stx and sty. On output stp is set to the new step.
370: bracket is a logical variable which specifies if a minimizer
371: has been bracketed. If the minimizer has not been bracketed
372: then on input bracket must be set false. If the minimizer
373: is bracketed then on output bracket is set true.
375: stpmin and stpmax are input variables which specify lower
376: and upper bounds for the step.
378: info is an integer output variable set as follows:
379: if info = 1,2,3,4,5, then the step has been computed
380: according to one of the five cases below. otherwise
381: info = 0, and this indicates improper input parameters.
383: subprograms called
385: fortran-supplied ... abs,max,min,sqrt
387: argonne national laboratory. minpack project. june 1983
388: jorge j. more', david j. thuente
390: */
392: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,PetscReal *stx,PetscReal *fx,PetscReal *dx,PetscReal *sty,PetscReal *fy,PetscReal *dy,PetscReal *stp,PetscReal *fp,PetscReal *dp)393: {
394: TaoLineSearch_MT *mtP = (TaoLineSearch_MT *) ls->data;
395: PetscReal gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
396: PetscInt bound;
399: /* Check the input parameters for errors */
400: mtP->infoc = 0;
401: if (mtP->bracket && (*stp <= PetscMin(*stx,*sty) || (*stp >= PetscMax(*stx,*sty)))) SETERRQ(PETSC_COMM_SELF,1,"bad stp in bracket");
402: if (*dx * (*stp-*stx) >= 0.0) SETERRQ(PETSC_COMM_SELF,1,"dx * (stp-stx) >= 0.0");
403: if (ls->stepmax < ls->stepmin) SETERRQ(PETSC_COMM_SELF,1,"stepmax > stepmin");
405: /* Determine if the derivatives have opposite sign */
406: sgnd = *dp * (*dx / PetscAbsReal(*dx));
408: if (*fp > *fx) {
409: /* Case 1: a higher function value.
410: The minimum is bracketed. If the cubic step is closer
411: to stx than the quadratic step, the cubic step is taken,
412: else the average of the cubic and quadratic steps is taken. */
414: mtP->infoc = 1;
415: bound = 1;
416: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
417: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
418: s = PetscMax(s,PetscAbsReal(*dp));
419: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
420: if (*stp < *stx) gamma1 = -gamma1;
421: /* Can p be 0? Check */
422: p = (gamma1 - *dx) + theta;
423: q = ((gamma1 - *dx) + gamma1) + *dp;
424: r = p/q;
425: stpc = *stx + r*(*stp - *stx);
426: stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);
428: if (PetscAbsReal(stpc-*stx) < PetscAbsReal(stpq-*stx)) {
429: stpf = stpc;
430: } else {
431: stpf = stpc + 0.5*(stpq - stpc);
432: }
433: mtP->bracket = 1;
434: } else if (sgnd < 0.0) {
435: /* Case 2: A lower function value and derivatives of
436: opposite sign. The minimum is bracketed. If the cubic
437: step is closer to stx than the quadratic (secant) step,
438: the cubic step is taken, else the quadratic step is taken. */
440: mtP->infoc = 2;
441: bound = 0;
442: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
443: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
444: s = PetscMax(s,PetscAbsReal(*dp));
445: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
446: if (*stp > *stx) gamma1 = -gamma1;
447: p = (gamma1 - *dp) + theta;
448: q = ((gamma1 - *dp) + gamma1) + *dx;
449: r = p/q;
450: stpc = *stp + r*(*stx - *stp);
451: stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);
453: if (PetscAbsReal(stpc-*stp) > PetscAbsReal(stpq-*stp)) {
454: stpf = stpc;
455: } else {
456: stpf = stpq;
457: }
458: mtP->bracket = 1;
459: } else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
460: /* Case 3: A lower function value, derivatives of the
461: same sign, and the magnitude of the derivative decreases.
462: The cubic step is only used if the cubic tends to infinity
463: in the direction of the step or if the minimum of the cubic
464: is beyond stp. Otherwise the cubic step is defined to be
465: either stepmin or stepmax. The quadratic (secant) step is also
466: computed and if the minimum is bracketed then the step
467: closest to stx is taken, else the step farthest away is taken. */
469: mtP->infoc = 3;
470: bound = 1;
471: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
472: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
473: s = PetscMax(s,PetscAbsReal(*dp));
475: /* The case gamma1 = 0 only arises if the cubic does not tend
476: to infinity in the direction of the step. */
477: gamma1 = s*PetscSqrtScalar(PetscMax(0.0,PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s)));
478: if (*stp > *stx) gamma1 = -gamma1;
479: p = (gamma1 - *dp) + theta;
480: q = (gamma1 + (*dx - *dp)) + gamma1;
481: r = p/q;
482: if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r*(*stx - *stp);
483: else if (*stp > *stx) stpc = ls->stepmax;
484: else stpc = ls->stepmin;
485: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
487: if (mtP->bracket) {
488: if (PetscAbsReal(*stp-stpc) < PetscAbsReal(*stp-stpq)) {
489: stpf = stpc;
490: } else {
491: stpf = stpq;
492: }
493: } else {
494: if (PetscAbsReal(*stp-stpc) > PetscAbsReal(*stp-stpq)) {
495: stpf = stpc;
496: } else {
497: stpf = stpq;
498: }
499: }
500: } else {
501: /* Case 4: A lower function value, derivatives of the
502: same sign, and the magnitude of the derivative does
503: not decrease. If the minimum is not bracketed, the step
504: is either stpmin or stpmax, else the cubic step is taken. */
506: mtP->infoc = 4;
507: bound = 0;
508: if (mtP->bracket) {
509: theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
510: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dy));
511: s = PetscMax(s,PetscAbsReal(*dp));
512: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dy/s)*(*dp/s));
513: if (*stp > *sty) gamma1 = -gamma1;
514: p = (gamma1 - *dp) + theta;
515: q = ((gamma1 - *dp) + gamma1) + *dy;
516: r = p/q;
517: stpc = *stp + r*(*sty - *stp);
518: stpf = stpc;
519: } else if (*stp > *stx) {
520: stpf = ls->stepmax;
521: } else {
522: stpf = ls->stepmin;
523: }
524: }
526: /* Update the interval of uncertainty. This update does not
527: depend on the new step or the case analysis above. */
529: if (*fp > *fx) {
530: *sty = *stp;
531: *fy = *fp;
532: *dy = *dp;
533: } else {
534: if (sgnd < 0.0) {
535: *sty = *stx;
536: *fy = *fx;
537: *dy = *dx;
538: }
539: *stx = *stp;
540: *fx = *fp;
541: *dx = *dp;
542: }
544: /* Compute the new step and safeguard it. */
545: stpf = PetscMin(ls->stepmax,stpf);
546: stpf = PetscMax(ls->stepmin,stpf);
547: *stp = stpf;
548: if (mtP->bracket && bound) {
549: if (*sty > *stx) {
550: *stp = PetscMin(*stx+0.66*(*sty-*stx),*stp);
551: } else {
552: *stp = PetscMax(*stx+0.66*(*sty-*stx),*stp);
553: }
554: }
555: return(0);
556: }