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: }
36: /* @ TaoApply_LineSearch - This routine takes step length of 1.0.
38: Input Parameters:
39: + tao - Tao context
40: . X - current iterate (on output X contains new iterate, X + step*S)
41: . f - objective function evaluated at X
42: . G - gradient evaluated at X
43: - D - search direction
46: Info is set to 0.
48: @ */
50: static PetscErrorCode TaoLineSearchApply_MT(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s) 51: {
52: PetscErrorCode ierr;
53: TaoLineSearch_MT *mt;
55: PetscReal xtrapf = 4.0;
56: PetscReal finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
57: PetscReal dgx, dgy, dg, dg2, fx, fy, stx, sty, dgtest;
58: PetscReal ftest1=0.0, ftest2=0.0;
59: PetscInt i, stage1,n1,n2,nn1,nn2;
60: PetscReal bstepmin1, bstepmin2, bstepmax;
61: PetscBool g_computed=PETSC_FALSE; /* to prevent extra gradient computation */
70: /* comm,type,size checks are done in interface TaoLineSearchApply */
71: mt = (TaoLineSearch_MT*)(ls->data);
72: ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
74: /* Check work vector */
75: if (!mt->work) {
76: VecDuplicate(x,&mt->work);
77: mt->x = x;
78: PetscObjectReference((PetscObject)mt->x);
79: } else if (x != mt->x) {
80: VecDestroy(&mt->work);
81: VecDuplicate(x,&mt->work);
82: PetscObjectDereference((PetscObject)mt->x);
83: mt->x = x;
84: PetscObjectReference((PetscObject)mt->x);
85: }
87: if (ls->bounded) {
88: /* Compute step length needed to make all variables equal a bound */
89: /* Compute the smallest steplength that will make one nonbinding variable
90: equal the bound */
91: VecGetLocalSize(ls->upper,&n1);
92: VecGetLocalSize(mt->x, &n2);
93: VecGetSize(ls->upper,&nn1);
94: VecGetSize(mt->x,&nn2);
95: if (n1 != n2 || nn1 != nn2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Variable vector not compatible with bounds vector");
96: VecScale(s,-1.0);
97: VecBoundGradientProjection(s,x,ls->lower,ls->upper,s);
98: VecScale(s,-1.0);
99: VecStepBoundInfo(x,s,ls->lower,ls->upper,&bstepmin1,&bstepmin2,&bstepmax);
100: ls->stepmax = PetscMin(bstepmax,1.0e15);
101: }
103: VecDot(g,s,&dginit);
104: if (PetscIsInfOrNanReal(dginit)) {
105: PetscInfo1(ls,"Initial Line Search step * g is Inf or Nan (%g)\n",(double)dginit);
106: ls->reason=TAOLINESEARCH_FAILED_INFORNAN;
107: return(0);
108: }
109: if (dginit >= 0.0) {
110: PetscInfo1(ls,"Initial Line Search step * g is not descent direction (%g)\n",(double)dginit);
111: ls->reason = TAOLINESEARCH_FAILED_ASCENT;
112: return(0);
113: }
115: /* Initialization */
116: mt->bracket = 0;
117: stage1 = 1;
118: finit = *f;
119: dgtest = ls->ftol * dginit;
120: width = ls->stepmax - ls->stepmin;
121: width1 = width * 2.0;
122: VecCopy(x,mt->work);
123: /* Variable dictionary:
124: stx, fx, dgx - the step, function, and derivative at the best step
125: sty, fy, dgy - the step, function, and derivative at the other endpoint
126: of the interval of uncertainty
127: step, f, dg - the step, function, and derivative at the current step */
129: stx = 0.0;
130: fx = finit;
131: dgx = dginit;
132: sty = 0.0;
133: fy = finit;
134: dgy = dginit;
136: ls->step=ls->initstep;
137: for (i=0; i< ls->max_funcs; i++) {
138: /* Set min and max steps to correspond to the interval of uncertainty */
139: if (mt->bracket) {
140: ls->stepmin = PetscMin(stx,sty);
141: ls->stepmax = PetscMax(stx,sty);
142: } else {
143: ls->stepmin = stx;
144: ls->stepmax = ls->step + xtrapf * (ls->step - stx);
145: }
147: /* Force the step to be within the bounds */
148: ls->step = PetscMax(ls->step,ls->stepmin);
149: ls->step = PetscMin(ls->step,ls->stepmax);
151: /* If an unusual termination is to occur, then let step be the lowest
152: point obtained thus far */
153: if ((stx!=0) && (((mt->bracket) && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) ||
154: ((ls->nfeval+ls->nfgeval) >= ls->max_funcs - 1) || (mt->infoc == 0))) {
155: ls->step = stx;
156: }
158: VecCopy(x,mt->work);
159: VecAXPY(mt->work,ls->step,s); /* W = X + step*S */
161: if (ls->bounded) {
162: VecMedian(ls->lower, mt->work, ls->upper, mt->work);
163: }
164: if (ls->usegts) {
165: TaoLineSearchComputeObjectiveAndGTS(ls,mt->work,f,&dg);
166: g_computed=PETSC_FALSE;
167: } else {
168: TaoLineSearchComputeObjectiveAndGradient(ls,mt->work,f,g);
169: g_computed=PETSC_TRUE;
170: if (ls->bounded) {
171: VecDot(g,x,&dg);
172: VecDot(g,mt->work,&dg2);
173: dg = (dg2 - dg)/ls->step;
174: } else {
175: VecDot(g,s,&dg);
176: }
177: }
179: if (0 == i) {
180: ls->f_fullstep=*f;
181: }
183: if (PetscIsInfOrNanReal(*f) || PetscIsInfOrNanReal(dg)) {
184: /* User provided compute function generated Not-a-Number, assume
185: domain violation and set function value and directional
186: derivative to infinity. */
187: *f = PETSC_INFINITY;
188: dg = PETSC_INFINITY;
189: }
191: ftest1 = finit + ls->step * dgtest;
192: if (ls->bounded) {
193: ftest2 = finit + ls->step * dgtest * ls->ftol;
194: }
195: /* Convergence testing */
196: if (((*f - ftest1 <= 1.0e-10 * PetscAbsReal(finit)) && (PetscAbsReal(dg) + ls->gtol*dginit <= 0.0))) {
197: PetscInfo(ls, "Line search success: Sufficient decrease and directional deriv conditions hold\n");
198: ls->reason = TAOLINESEARCH_SUCCESS;
199: break;
200: }
202: /* Check Armijo if beyond the first breakpoint */
203: if (ls->bounded && (*f <= ftest2) && (ls->step >= bstepmin2)) {
204: PetscInfo(ls,"Line search success: Sufficient decrease.\n");
205: ls->reason = TAOLINESEARCH_SUCCESS;
206: break;
207: }
209: /* Checks for bad cases */
210: if (((mt->bracket) && (ls->step <= ls->stepmin||ls->step >= ls->stepmax)) || (!mt->infoc)) {
211: PetscInfo(ls,"Rounding errors may prevent further progress. May not be a step satisfying\n");
212: PetscInfo(ls,"sufficient decrease and curvature conditions. Tolerances may be too small.\n");
213: ls->reason = TAOLINESEARCH_HALTED_OTHER;
214: break;
215: }
216: if ((ls->step == ls->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
217: PetscInfo1(ls,"Step is at the upper bound, stepmax (%g)\n",(double)ls->stepmax);
218: ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
219: break;
220: }
221: if ((ls->step == ls->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
222: PetscInfo1(ls,"Step is at the lower bound, stepmin (%g)\n",(double)ls->stepmin);
223: ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
224: break;
225: }
226: if ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol*ls->stepmax)){
227: PetscInfo1(ls,"Relative width of interval of uncertainty is at most rtol (%g)\n",(double)ls->rtol);
228: ls->reason = TAOLINESEARCH_HALTED_RTOL;
229: break;
230: }
232: /* In the first stage, we seek a step for which the modified function
233: has a nonpositive value and nonnegative derivative */
234: if ((stage1) && (*f <= ftest1) && (dg >= dginit * PetscMin(ls->ftol, ls->gtol))) {
235: stage1 = 0;
236: }
238: /* A modified function is used to predict the step only if we
239: have not obtained a step for which the modified function has a
240: nonpositive function value and nonnegative derivative, and if a
241: lower function value has been obtained but the decrease is not
242: sufficient */
244: if ((stage1) && (*f <= fx) && (*f > ftest1)) {
245: fm = *f - ls->step * dgtest; /* Define modified function */
246: fxm = fx - stx * dgtest; /* and derivatives */
247: fym = fy - sty * dgtest;
248: dgm = dg - dgtest;
249: dgxm = dgx - dgtest;
250: dgym = dgy - dgtest;
252: /* if (dgxm * (ls->step - stx) >= 0.0) */
253: /* Update the interval of uncertainty and compute the new step */
254: Tao_mcstep(ls,&stx,&fxm,&dgxm,&sty,&fym,&dgym,&ls->step,&fm,&dgm);
256: fx = fxm + stx * dgtest; /* Reset the function and */
257: fy = fym + sty * dgtest; /* gradient values */
258: dgx = dgxm + dgtest;
259: dgy = dgym + dgtest;
260: } else {
261: /* Update the interval of uncertainty and compute the new step */
262: Tao_mcstep(ls,&stx,&fx,&dgx,&sty,&fy,&dgy,&ls->step,f,&dg);
263: }
265: /* Force a sufficient decrease in the interval of uncertainty */
266: if (mt->bracket) {
267: if (PetscAbsReal(sty - stx) >= 0.66 * width1) ls->step = stx + 0.5*(sty - stx);
268: width1 = width;
269: width = PetscAbsReal(sty - stx);
270: }
271: }
272: if ((ls->nfeval+ls->nfgeval) > ls->max_funcs) {
273: PetscInfo2(ls,"Number of line search function evals (%D) > maximum (%D)\n",(ls->nfeval+ls->nfgeval),ls->max_funcs);
274: ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
275: }
277: /* Finish computations */
278: PetscInfo2(ls,"%D function evals in line search, step = %g\n",(ls->nfeval+ls->nfgeval),(double)ls->step);
280: /* Set new solution vector and compute gradient if needed */
281: VecCopy(mt->work,x);
282: if (!g_computed) {
283: TaoLineSearchComputeGradient(ls,mt->work,g);
284: }
285: return(0);
286: }
288: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_MT(TaoLineSearch ls)289: {
290: PetscErrorCode ierr;
291: TaoLineSearch_MT *ctx;
295: PetscNewLog(ls,&ctx);
296: ctx->bracket=0;
297: ctx->infoc=1;
298: ls->data = (void*)ctx;
299: ls->initstep = 1.0;
300: ls->ops->setup=0;
301: ls->ops->reset=0;
302: ls->ops->apply=TaoLineSearchApply_MT;
303: ls->ops->destroy=TaoLineSearchDestroy_MT;
304: ls->ops->setfromoptions=TaoLineSearchSetFromOptions_MT;
305: return(0);
306: }
308: /*
309: The subroutine mcstep is taken from the work of Jorge Nocedal.
310: this is a variant of More' and Thuente's routine.
312: subroutine mcstep
314: the purpose of mcstep is to compute a safeguarded step for
315: a linesearch and to update an interval of uncertainty for
316: a minimizer of the function.
318: the parameter stx contains the step with the least function
319: value. the parameter stp contains the current step. it is
320: assumed that the derivative at stx is negative in the
321: direction of the step. if bracket is set true then a
322: minimizer has been bracketed in an interval of uncertainty
323: with endpoints stx and sty.
325: the subroutine statement is
327: subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
328: stpmin,stpmax,info)
330: where
332: stx, fx, and dx are variables which specify the step,
333: the function, and the derivative at the best step obtained
334: so far. The derivative must be negative in the direction
335: of the step, that is, dx and stp-stx must have opposite
336: signs. On output these parameters are updated appropriately.
338: sty, fy, and dy are variables which specify the step,
339: the function, and the derivative at the other endpoint of
340: the interval of uncertainty. On output these parameters are
341: updated appropriately.
343: stp, fp, and dp are variables which specify the step,
344: the function, and the derivative at the current step.
345: If bracket is set true then on input stp must be
346: between stx and sty. On output stp is set to the new step.
348: bracket is a logical variable which specifies if a minimizer
349: has been bracketed. If the minimizer has not been bracketed
350: then on input bracket must be set false. If the minimizer
351: is bracketed then on output bracket is set true.
353: stpmin and stpmax are input variables which specify lower
354: and upper bounds for the step.
356: info is an integer output variable set as follows:
357: if info = 1,2,3,4,5, then the step has been computed
358: according to one of the five cases below. otherwise
359: info = 0, and this indicates improper input parameters.
361: subprograms called
363: fortran-supplied ... abs,max,min,sqrt
365: argonne national laboratory. minpack project. june 1983
366: jorge j. more', david j. thuente
368: */
370: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,PetscReal *stx,PetscReal *fx,PetscReal *dx,PetscReal *sty,PetscReal *fy,PetscReal *dy,PetscReal *stp,PetscReal *fp,PetscReal *dp)371: {
372: TaoLineSearch_MT *mtP = (TaoLineSearch_MT *) ls->data;
373: PetscReal gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
374: PetscInt bound;
377: /* Check the input parameters for errors */
378: mtP->infoc = 0;
379: if (mtP->bracket && (*stp <= PetscMin(*stx,*sty) || (*stp >= PetscMax(*stx,*sty)))) SETERRQ(PETSC_COMM_SELF,1,"bad stp in bracket");
380: if (*dx * (*stp-*stx) >= 0.0) SETERRQ(PETSC_COMM_SELF,1,"dx * (stp-stx) >= 0.0");
381: if (ls->stepmax < ls->stepmin) SETERRQ(PETSC_COMM_SELF,1,"stepmax > stepmin");
383: /* Determine if the derivatives have opposite sign */
384: sgnd = *dp * (*dx / PetscAbsReal(*dx));
386: if (*fp > *fx) {
387: /* Case 1: a higher function value.
388: The minimum is bracketed. If the cubic step is closer
389: to stx than the quadratic step, the cubic step is taken,
390: else the average of the cubic and quadratic steps is taken. */
392: mtP->infoc = 1;
393: bound = 1;
394: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
395: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
396: s = PetscMax(s,PetscAbsReal(*dp));
397: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
398: if (*stp < *stx) gamma1 = -gamma1;
399: /* Can p be 0? Check */
400: p = (gamma1 - *dx) + theta;
401: q = ((gamma1 - *dx) + gamma1) + *dp;
402: r = p/q;
403: stpc = *stx + r*(*stp - *stx);
404: stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);
406: if (PetscAbsReal(stpc-*stx) < PetscAbsReal(stpq-*stx)) {
407: stpf = stpc;
408: } else {
409: stpf = stpc + 0.5*(stpq - stpc);
410: }
411: mtP->bracket = 1;
412: } else if (sgnd < 0.0) {
413: /* Case 2: A lower function value and derivatives of
414: opposite sign. The minimum is bracketed. If the cubic
415: step is closer to stx than the quadratic (secant) step,
416: the cubic step is taken, else the quadratic step is taken. */
418: mtP->infoc = 2;
419: bound = 0;
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: p = (gamma1 - *dp) + theta;
426: q = ((gamma1 - *dp) + gamma1) + *dx;
427: r = p/q;
428: stpc = *stp + r*(*stx - *stp);
429: stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);
431: if (PetscAbsReal(stpc-*stp) > PetscAbsReal(stpq-*stp)) {
432: stpf = stpc;
433: } else {
434: stpf = stpq;
435: }
436: mtP->bracket = 1;
437: } else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
438: /* Case 3: A lower function value, derivatives of the
439: same sign, and the magnitude of the derivative decreases.
440: The cubic step is only used if the cubic tends to infinity
441: in the direction of the step or if the minimum of the cubic
442: is beyond stp. Otherwise the cubic step is defined to be
443: either stepmin or stepmax. The quadratic (secant) step is also
444: computed and if the minimum is bracketed then the step
445: closest to stx is taken, else the step farthest away is taken. */
447: mtP->infoc = 3;
448: bound = 1;
449: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
450: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
451: s = PetscMax(s,PetscAbsReal(*dp));
453: /* The case gamma1 = 0 only arises if the cubic does not tend
454: to infinity in the direction of the step. */
455: gamma1 = s*PetscSqrtScalar(PetscMax(0.0,PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s)));
456: if (*stp > *stx) gamma1 = -gamma1;
457: p = (gamma1 - *dp) + theta;
458: q = (gamma1 + (*dx - *dp)) + gamma1;
459: r = p/q;
460: if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r*(*stx - *stp);
461: else if (*stp > *stx) stpc = ls->stepmax;
462: else stpc = ls->stepmin;
463: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
465: if (mtP->bracket) {
466: if (PetscAbsReal(*stp-stpc) < PetscAbsReal(*stp-stpq)) {
467: stpf = stpc;
468: } else {
469: stpf = stpq;
470: }
471: } else {
472: if (PetscAbsReal(*stp-stpc) > PetscAbsReal(*stp-stpq)) {
473: stpf = stpc;
474: } else {
475: stpf = stpq;
476: }
477: }
478: } else {
479: /* Case 4: A lower function value, derivatives of the
480: same sign, and the magnitude of the derivative does
481: not decrease. If the minimum is not bracketed, the step
482: is either stpmin or stpmax, else the cubic step is taken. */
484: mtP->infoc = 4;
485: bound = 0;
486: if (mtP->bracket) {
487: theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
488: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dy));
489: s = PetscMax(s,PetscAbsReal(*dp));
490: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dy/s)*(*dp/s));
491: if (*stp > *sty) gamma1 = -gamma1;
492: p = (gamma1 - *dp) + theta;
493: q = ((gamma1 - *dp) + gamma1) + *dy;
494: r = p/q;
495: stpc = *stp + r*(*sty - *stp);
496: stpf = stpc;
497: } else if (*stp > *stx) {
498: stpf = ls->stepmax;
499: } else {
500: stpf = ls->stepmin;
501: }
502: }
504: /* Update the interval of uncertainty. This update does not
505: depend on the new step or the case analysis above. */
507: if (*fp > *fx) {
508: *sty = *stp;
509: *fy = *fp;
510: *dy = *dp;
511: } else {
512: if (sgnd < 0.0) {
513: *sty = *stx;
514: *fy = *fx;
515: *dy = *dx;
516: }
517: *stx = *stp;
518: *fx = *fp;
519: *dx = *dp;
520: }
522: /* Compute the new step and safeguard it. */
523: stpf = PetscMin(ls->stepmax,stpf);
524: stpf = PetscMax(ls->stepmin,stpf);
525: *stp = stpf;
526: if (mtP->bracket && bound) {
527: if (*sty > *stx) {
528: *stp = PetscMin(*stx+0.66*(*sty-*stx),*stp);
529: } else {
530: *stp = PetscMax(*stx+0.66*(*sty-*stx),*stp);
531: }
532: }
533: return(0);
534: }