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