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: }
116: /* Initialization */
117: mt->bracket = 0;
118: stage1 = 1;
119: finit = *f;
120: dgtest = ls->ftol * dginit;
121: width = ls->stepmax - ls->stepmin;
122: width1 = width * 2.0;
123: VecCopy(x,mt->work);
124: /* Variable dictionary:
125: stx, fx, dgx - the step, function, and derivative at the best step
126: sty, fy, dgy - the step, function, and derivative at the other endpoint
127: of the interval of uncertainty
128: step, f, dg - the step, function, and derivative at the current step */
130: stx = 0.0;
131: fx = finit;
132: dgx = dginit;
133: sty = 0.0;
134: fy = finit;
135: dgy = dginit;
137: ls->step=ls->initstep;
138: for (i=0; i< ls->max_funcs; i++) {
139: /* Set min and max steps to correspond to the interval of uncertainty */
140: if (mt->bracket) {
141: ls->stepmin = PetscMin(stx,sty);
142: ls->stepmax = PetscMax(stx,sty);
143: } else {
144: ls->stepmin = stx;
145: ls->stepmax = ls->step + xtrapf * (ls->step - stx);
146: }
148: /* Force the step to be within the bounds */
149: ls->step = PetscMax(ls->step,ls->stepmin);
150: ls->step = PetscMin(ls->step,ls->stepmax);
152: /* If an unusual termination is to occur, then let step be the lowest
153: point obtained thus far */
154: if ((stx!=0) && (((mt->bracket) && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) ||
155: ((ls->nfeval+ls->nfgeval) >= ls->max_funcs - 1) || (mt->infoc == 0))) {
156: ls->step = stx;
157: }
159: VecCopy(x,mt->work);
160: VecAXPY(mt->work,ls->step,s); /* W = X + step*S */
162: if (ls->bounded) {
163: VecMedian(ls->lower, mt->work, ls->upper, mt->work);
164: }
165: if (ls->usegts) {
166: TaoLineSearchComputeObjectiveAndGTS(ls,mt->work,f,&dg);
167: g_computed=PETSC_FALSE;
168: } else {
169: TaoLineSearchComputeObjectiveAndGradient(ls,mt->work,f,g);
170: g_computed=PETSC_TRUE;
171: if (ls->bounded) {
172: VecDot(g,x,&dg);
173: VecDot(g,mt->work,&dg2);
174: dg = (dg2 - dg)/ls->step;
175: } else {
176: VecDot(g,s,&dg);
177: }
178: }
180: if (0 == i) {
181: ls->f_fullstep=*f;
182: }
184: if (PetscIsInfOrNanReal(*f) || PetscIsInfOrNanReal(dg)) {
185: /* User provided compute function generated Not-a-Number, assume
186: domain violation and set function value and directional
187: derivative to infinity. */
188: *f = PETSC_INFINITY;
189: dg = PETSC_INFINITY;
190: }
192: ftest1 = finit + ls->step * dgtest;
193: if (ls->bounded) {
194: ftest2 = finit + ls->step * dgtest * ls->ftol;
195: }
196: /* Convergence testing */
197: if (((*f - ftest1 <= 1.0e-10 * PetscAbsReal(finit)) && (PetscAbsReal(dg) + ls->gtol*dginit <= 0.0))) {
198: PetscInfo(ls, "Line search success: Sufficient decrease and directional deriv conditions hold\n");
199: ls->reason = TAOLINESEARCH_SUCCESS;
200: break;
201: }
203: /* Check Armijo if beyond the first breakpoint */
204: if (ls->bounded && (*f <= ftest2) && (ls->step >= bstepmin2)) {
205: PetscInfo(ls,"Line search success: Sufficient decrease.\n");
206: ls->reason = TAOLINESEARCH_SUCCESS;
207: break;
208: }
210: /* Checks for bad cases */
211: if (((mt->bracket) && (ls->step <= ls->stepmin||ls->step >= ls->stepmax)) || (!mt->infoc)) {
212: PetscInfo(ls,"Rounding errors may prevent further progress. May not be a step satisfying\n");
213: PetscInfo(ls,"sufficient decrease and curvature conditions. Tolerances may be too small.\n");
214: ls->reason = TAOLINESEARCH_HALTED_OTHER;
215: break;
216: }
217: if ((ls->step == ls->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
218: PetscInfo1(ls,"Step is at the upper bound, stepmax (%g)\n",(double)ls->stepmax);
219: ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
220: break;
221: }
222: if ((ls->step == ls->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
223: PetscInfo1(ls,"Step is at the lower bound, stepmin (%g)\n",(double)ls->stepmin);
224: ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
225: break;
226: }
227: if ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol*ls->stepmax)){
228: PetscInfo1(ls,"Relative width of interval of uncertainty is at most rtol (%g)\n",(double)ls->rtol);
229: ls->reason = TAOLINESEARCH_HALTED_RTOL;
230: break;
231: }
233: /* In the first stage, we seek a step for which the modified function
234: has a nonpositive value and nonnegative derivative */
235: if ((stage1) && (*f <= ftest1) && (dg >= dginit * PetscMin(ls->ftol, ls->gtol))) {
236: stage1 = 0;
237: }
239: /* A modified function is used to predict the step only if we
240: have not obtained a step for which the modified function has a
241: nonpositive function value and nonnegative derivative, and if a
242: lower function value has been obtained but the decrease is not
243: sufficient */
245: if ((stage1) && (*f <= fx) && (*f > ftest1)) {
246: fm = *f - ls->step * dgtest; /* Define modified function */
247: fxm = fx - stx * dgtest; /* and derivatives */
248: fym = fy - sty * dgtest;
249: dgm = dg - dgtest;
250: dgxm = dgx - dgtest;
251: dgym = dgy - dgtest;
253: /* if (dgxm * (ls->step - stx) >= 0.0) */
254: /* Update the interval of uncertainty and compute the new step */
255: Tao_mcstep(ls,&stx,&fxm,&dgxm,&sty,&fym,&dgym,&ls->step,&fm,&dgm);
257: fx = fxm + stx * dgtest; /* Reset the function and */
258: fy = fym + sty * dgtest; /* gradient values */
259: dgx = dgxm + dgtest;
260: dgy = dgym + dgtest;
261: } else {
262: /* Update the interval of uncertainty and compute the new step */
263: Tao_mcstep(ls,&stx,&fx,&dgx,&sty,&fy,&dgy,&ls->step,f,&dg);
264: }
266: /* Force a sufficient decrease in the interval of uncertainty */
267: if (mt->bracket) {
268: if (PetscAbsReal(sty - stx) >= 0.66 * width1) ls->step = stx + 0.5*(sty - stx);
269: width1 = width;
270: width = PetscAbsReal(sty - stx);
271: }
272: }
273: if ((ls->nfeval+ls->nfgeval) > ls->max_funcs) {
274: PetscInfo2(ls,"Number of line search function evals (%D) > maximum (%D)\n",(ls->nfeval+ls->nfgeval),ls->max_funcs);
275: ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
276: }
278: /* Finish computations */
279: PetscInfo2(ls,"%D function evals in line search, step = %g\n",(ls->nfeval+ls->nfgeval),(double)ls->step);
281: /* Set new solution vector and compute gradient if needed */
282: VecCopy(mt->work,x);
283: if (!g_computed) {
284: TaoLineSearchComputeGradient(ls,mt->work,g);
285: }
286: return(0);
287: }
289: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_MT(TaoLineSearch ls)290: {
291: PetscErrorCode ierr;
292: TaoLineSearch_MT *ctx;
296: PetscNewLog(ls,&ctx);
297: ctx->bracket=0;
298: ctx->infoc=1;
299: ls->data = (void*)ctx;
300: ls->initstep = 1.0;
301: ls->ops->setup=0;
302: ls->ops->reset=0;
303: ls->ops->apply=TaoLineSearchApply_MT;
304: ls->ops->destroy=TaoLineSearchDestroy_MT;
305: ls->ops->setfromoptions=TaoLineSearchSetFromOptions_MT;
306: return(0);
307: }
309: /*
310: The subroutine mcstep is taken from the work of Jorge Nocedal.
311: this is a variant of More' and Thuente's routine.
313: subroutine mcstep
315: the purpose of mcstep is to compute a safeguarded step for
316: a linesearch and to update an interval of uncertainty for
317: a minimizer of the function.
319: the parameter stx contains the step with the least function
320: value. the parameter stp contains the current step. it is
321: assumed that the derivative at stx is negative in the
322: direction of the step. if bracket is set true then a
323: minimizer has been bracketed in an interval of uncertainty
324: with endpoints stx and sty.
326: the subroutine statement is
328: subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
329: stpmin,stpmax,info)
331: where
333: stx, fx, and dx are variables which specify the step,
334: the function, and the derivative at the best step obtained
335: so far. The derivative must be negative in the direction
336: of the step, that is, dx and stp-stx must have opposite
337: signs. On output these parameters are updated appropriately.
339: sty, fy, and dy are variables which specify the step,
340: the function, and the derivative at the other endpoint of
341: the interval of uncertainty. On output these parameters are
342: updated appropriately.
344: stp, fp, and dp are variables which specify the step,
345: the function, and the derivative at the current step.
346: If bracket is set true then on input stp must be
347: between stx and sty. On output stp is set to the new step.
349: bracket is a logical variable which specifies if a minimizer
350: has been bracketed. If the minimizer has not been bracketed
351: then on input bracket must be set false. If the minimizer
352: is bracketed then on output bracket is set true.
354: stpmin and stpmax are input variables which specify lower
355: and upper bounds for the step.
357: info is an integer output variable set as follows:
358: if info = 1,2,3,4,5, then the step has been computed
359: according to one of the five cases below. otherwise
360: info = 0, and this indicates improper input parameters.
362: subprograms called
364: fortran-supplied ... abs,max,min,sqrt
366: argonne national laboratory. minpack project. june 1983
367: jorge j. more', david j. thuente
369: */
371: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,PetscReal *stx,PetscReal *fx,PetscReal *dx,PetscReal *sty,PetscReal *fy,PetscReal *dy,PetscReal *stp,PetscReal *fp,PetscReal *dp)372: {
373: TaoLineSearch_MT *mtP = (TaoLineSearch_MT *) ls->data;
374: PetscReal gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
375: PetscInt bound;
378: /* Check the input parameters for errors */
379: mtP->infoc = 0;
380: if (mtP->bracket && (*stp <= PetscMin(*stx,*sty) || (*stp >= PetscMax(*stx,*sty)))) SETERRQ(PETSC_COMM_SELF,1,"bad stp in bracket");
381: if (*dx * (*stp-*stx) >= 0.0) SETERRQ(PETSC_COMM_SELF,1,"dx * (stp-stx) >= 0.0");
382: if (ls->stepmax < ls->stepmin) SETERRQ(PETSC_COMM_SELF,1,"stepmax > stepmin");
384: /* Determine if the derivatives have opposite sign */
385: sgnd = *dp * (*dx / PetscAbsReal(*dx));
387: if (*fp > *fx) {
388: /* Case 1: a higher function value.
389: The minimum is bracketed. If the cubic step is closer
390: to stx than the quadratic step, the cubic step is taken,
391: else the average of the cubic and quadratic steps is taken. */
393: mtP->infoc = 1;
394: bound = 1;
395: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
396: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
397: s = PetscMax(s,PetscAbsReal(*dp));
398: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
399: if (*stp < *stx) gamma1 = -gamma1;
400: /* Can p be 0? Check */
401: p = (gamma1 - *dx) + theta;
402: q = ((gamma1 - *dx) + gamma1) + *dp;
403: r = p/q;
404: stpc = *stx + r*(*stp - *stx);
405: stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);
407: if (PetscAbsReal(stpc-*stx) < PetscAbsReal(stpq-*stx)) {
408: stpf = stpc;
409: } else {
410: stpf = stpc + 0.5*(stpq - stpc);
411: }
412: mtP->bracket = 1;
413: } else if (sgnd < 0.0) {
414: /* Case 2: A lower function value and derivatives of
415: opposite sign. The minimum is bracketed. If the cubic
416: step is closer to stx than the quadratic (secant) step,
417: the cubic step is taken, else the quadratic step is taken. */
419: mtP->infoc = 2;
420: bound = 0;
421: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
422: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
423: s = PetscMax(s,PetscAbsReal(*dp));
424: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
425: if (*stp > *stx) gamma1 = -gamma1;
426: p = (gamma1 - *dp) + theta;
427: q = ((gamma1 - *dp) + gamma1) + *dx;
428: r = p/q;
429: stpc = *stp + r*(*stx - *stp);
430: stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);
432: if (PetscAbsReal(stpc-*stp) > PetscAbsReal(stpq-*stp)) {
433: stpf = stpc;
434: } else {
435: stpf = stpq;
436: }
437: mtP->bracket = 1;
438: } else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
439: /* Case 3: A lower function value, derivatives of the
440: same sign, and the magnitude of the derivative decreases.
441: The cubic step is only used if the cubic tends to infinity
442: in the direction of the step or if the minimum of the cubic
443: is beyond stp. Otherwise the cubic step is defined to be
444: either stepmin or stepmax. The quadratic (secant) step is also
445: computed and if the minimum is bracketed then the step
446: closest to stx is taken, else the step farthest away is taken. */
448: mtP->infoc = 3;
449: bound = 1;
450: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
451: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
452: s = PetscMax(s,PetscAbsReal(*dp));
454: /* The case gamma1 = 0 only arises if the cubic does not tend
455: to infinity in the direction of the step. */
456: gamma1 = s*PetscSqrtScalar(PetscMax(0.0,PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s)));
457: if (*stp > *stx) gamma1 = -gamma1;
458: p = (gamma1 - *dp) + theta;
459: q = (gamma1 + (*dx - *dp)) + gamma1;
460: r = p/q;
461: if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r*(*stx - *stp);
462: else if (*stp > *stx) stpc = ls->stepmax;
463: else stpc = ls->stepmin;
464: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
466: if (mtP->bracket) {
467: if (PetscAbsReal(*stp-stpc) < PetscAbsReal(*stp-stpq)) {
468: stpf = stpc;
469: } else {
470: stpf = stpq;
471: }
472: } else {
473: if (PetscAbsReal(*stp-stpc) > PetscAbsReal(*stp-stpq)) {
474: stpf = stpc;
475: } else {
476: stpf = stpq;
477: }
478: }
479: } else {
480: /* Case 4: A lower function value, derivatives of the
481: same sign, and the magnitude of the derivative does
482: not decrease. If the minimum is not bracketed, the step
483: is either stpmin or stpmax, else the cubic step is taken. */
485: mtP->infoc = 4;
486: bound = 0;
487: if (mtP->bracket) {
488: theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
489: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dy));
490: s = PetscMax(s,PetscAbsReal(*dp));
491: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dy/s)*(*dp/s));
492: if (*stp > *sty) gamma1 = -gamma1;
493: p = (gamma1 - *dp) + theta;
494: q = ((gamma1 - *dp) + gamma1) + *dy;
495: r = p/q;
496: stpc = *stp + r*(*sty - *stp);
497: stpf = stpc;
498: } else if (*stp > *stx) {
499: stpf = ls->stepmax;
500: } else {
501: stpf = ls->stepmin;
502: }
503: }
505: /* Update the interval of uncertainty. This update does not
506: depend on the new step or the case analysis above. */
508: if (*fp > *fx) {
509: *sty = *stp;
510: *fy = *fp;
511: *dy = *dp;
512: } else {
513: if (sgnd < 0.0) {
514: *sty = *stx;
515: *fy = *fx;
516: *dy = *dx;
517: }
518: *stx = *stp;
519: *fx = *fp;
520: *dx = *dp;
521: }
523: /* Compute the new step and safeguard it. */
524: stpf = PetscMin(ls->stepmax,stpf);
525: stpf = PetscMax(ls->stepmin,stpf);
526: *stp = stpf;
527: if (mtP->bracket && bound) {
528: if (*sty > *stx) {
529: *stp = PetscMin(*stx+0.66*(*sty-*stx),*stp);
530: } else {
531: *stp = PetscMax(*stx+0.66*(*sty-*stx),*stp);
532: }
533: }
534: return(0);
535: }