Actual source code: boundmorethuente.c
1: /*$Id: morethuente.c,v 1.1 2002/08/11 01:15:06 benson Exp $*/
3: #include "morethuente.h" /*I "tao_solver.h
4: " I*/
6: static int TaoSetOptions_BoundLineSearch(TAO_SOLVER,void*);
7: static int TaoApply_BoundLineSearch(TAO_SOLVER,TaoVec*,TaoVec*,TaoVec*,TaoVec*,
8: double*,double*,int*,void*);
10: /*------------------------------------------------------------*/
13: static int TaoSetOptions_BoundLineSearch(TAO_SOLVER tao, void*linectx)
14: {
15: TAO_LINESEARCH *ctx = (TAO_LINESEARCH *)linectx;
16: int info;
18: TaoFunctionBegin;
19: info = TaoOptionsHead("Modified More-Thuente line line search options for bound constrained minimization");CHKERRQ(info);
20: info = TaoOptionsTail();CHKERRQ(info);
21: info = TaoSetOptions_LineSearch(tao,ctx);
22: TaoFunctionReturn(0);
23: }
24: /*------------------------------------------------------------*/
26: /* ---------------------------------------------------------- */
29: /* @
30: TaoApply_BoundLineSearch - This routine performs a line search algorithm.
32: Input Parameters:
33: + tao - TAO_SOLVER context
34: . X - current iterate (on output X contains new iterate, X + step*S)
35: . S - search direction
36: . f - objective function evaluated at X
37: . G - gradient evaluated at X
38: . W - work vector
39: - step - initial estimate of step length
41: Output parameters:
42: + f - objective function evaluated at new iterate, X + step*S
43: . G - gradient evaluated at new iterate, X + step*S
44: . X - new iterate
45: . gnorm - 2-norm of G
46: - step - final step length
49: Info is set to one of:
50: + 0 - the line search succeeds; the sufficient decrease
51: condition and the directional derivative condition hold
53: negative number if an input parameter is invalid
54: . -1 - step < 0
55: . -2 - ftol < 0
56: . -3 - rtol < 0
57: . -4 - gtol < 0
58: . -5 - stepmin < 0
59: . -6 - stepmax < stepmin
60: - -7 - maxfev < 0
62: positive number > 1 if the line search otherwise terminates
63: + 2 - Relative width of the interval of uncertainty is
64: at most rtol.
65: . 3 - Maximum number of function evaluations (maxfev) has
66: been reached.
67: . 4 - Step is at the lower bound, stepmin.
68: . 5 - Step is at the upper bound, stepmax.
69: . 6 - Rounding errors may prevent further progress.
70: There may not be a step that satisfies the
71: sufficient decrease and curvature conditions.
72: Tolerances may be too small.
73: - 7 - Search direction is not a descent direction.
75: Notes:
76: This algorithm is is a modification of the algorithm by More' and Thuente.
77: The modifications concern bounds. This algorithm steps in the direction
78: passed into this routine. This point get projected back into the feasible set.
79: In the context of bound constrained
80: optimization, there may not be a point in the piecewise linear
81: path that satisfies the Wolfe conditions. When the active set
82: is changing, decrease in the objective function may be sufficient
83: to terminate this line search.
85: Note: Much of this coded is identical to the More' Thuente line search.
86: Variations to the code are commented.
88: Notes:
89: This routine is used within the following TAO boiund constrained minimization
90: solvers: Newton linesearch (tao_tron) and limited memory variable metric
91: (tao_blmvm).
93: Level: advanced
95: .keywords: TAO_SOLVER, linesearch
96: @ */
97: static int TaoApply_BoundLineSearch(TAO_SOLVER tao,TaoVec* X,TaoVec* G,TaoVec* S,TaoVec* W,double *f,
98: double *step,int *info2,void*ctx)
99: {
100: TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
101: TaoVec *XL,*XU;
102: double zero = 0.0, two = 2.0, p5 = 0.5, p66 = 0.66, xtrapf = 4.0;
103: double finit, width, width1, dginit,fm, fxm, fym, dgm, dgxm, dgym;
104: double dgx, dgy, dg, fx, fy, stx, sty, dgtest, ftest1=0.0,ftest2=0.0;
105: double bstepmin1,bstepmin2,bstepmax;
106: double dg1,dg2;
107: int info, i, stage1;
109: #if defined(PETSC_USE_COMPLEX)
110: PetscScalar cdginit, cdg, cstep = 0.0;
111: #endif
113: TaoFunctionBegin;
114: /* neP->stepmin - lower bound for step */
115: /* neP->stepmax - upper bound for step */
116: /* neP->rtol - relative tolerance for an acceptable step */
117: /* neP->ftol - tolerance for sufficient decrease condition */
118: /* neP->gtol - tolerance for curvature condition */
119: /* neP->nfev - number of function evaluations */
120: /* neP->maxfev - maximum number of function evaluations */
122: /* Check input parameters for errors */
123: if (*step < zero) {
124: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: step (%g) < 0\n",*step)); CHKERRQ(info);
125: *info2 = -1; TaoFunctionReturn(0);
126: } else if (neP->ftol < zero) {
127: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: ftol (%g) < 0\n",neP->ftol)); CHKERRQ(info);
128: *info2 = -2; TaoFunctionReturn(0);
129: } else if (neP->rtol < zero) {
130: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: rtol (%g) < 0\n",neP->rtol)); CHKERRQ(info);
131: *info2 = -3; TaoFunctionReturn(0);
132: } else if (neP->gtol < zero) {
133: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: gtol (%g) < 0\n",neP->gtol)); CHKERRQ(info);
134: *info2 = -4; TaoFunctionReturn(0);
135: } else if (neP->stepmin < zero) {
136: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: stepmin (%g) < 0\n",neP->stepmin)); CHKERRQ(info);
137: *info2 = -5; TaoFunctionReturn(0);
138: } else if (neP->stepmax < neP->stepmin) {
139: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: stepmax (%g) < stepmin (%g)\n",
140: neP->stepmax,neP->stepmin)); CHKERRQ(info);
141: *info2 = -6; TaoFunctionReturn(0);
142: } else if (neP->maxfev < zero) {
143: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: maxfev (%d) < 0\n",neP->maxfev)); CHKERRQ(info);
144: *info2 = -7; TaoFunctionReturn(0);
145: }
147: /* Next 6 lines added by S. Benson to accomodate bounds on the variables */
148: /* Compute step length needed to make all variables equal a bound */
149: /* Compute the smallest steplength that will make one nonbinding variable equal the bound */
150: info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
151: info = S->Negate(); CHKERRQ(info);
152: info = S->BoundGradientProjection(S,XL,X,XU); CHKERRQ(info);
153: info = S->Negate(); CHKERRQ(info);
154: info = X->StepBoundInfo(XL,XU,S,&bstepmin1,&bstepmin2,&bstepmax); CHKERRQ(info);
155: neP->stepmax=TaoMin(bstepmax,1.0e+15);
157: /* Check that search direction is a descent direction */
159: #if defined(PETSC_USE_COMPLEX)
160: info = G->Dot(S,&cdginit);CHKERRQ(info); dginit = TaoReal(cdginit);
161: #else
162: info = G->Dot(S,&dginit);CHKERRQ(info);
163: #endif
166: if (dginit >= zero) {
167: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Search direction not a descent direction\n")); CHKERRQ(info);
168: *info2 = 7; TaoFunctionReturn(0);
169: }
171: /* Initialization */
172: neP->bracket = 0;
173: *info2 = 0;
174: stage1 = 1;
175: finit = *f;
176: dgtest = neP->ftol * dginit;
177: width = neP->stepmax - neP->stepmin;
178: width1 = width * two;
179: info = W->CopyFrom(X);CHKERRQ(info);
180: /* Variable dictionary:
181: stx, fx, dgx - the step, function, and derivative at the best step
182: sty, fy, dgy - the step, function, and derivative at the other endpoint
183: of the interval of uncertainty
184: step, f, dg - the step, function, and derivative at the current step */
186: stx = zero;
187: fx = finit;
188: dgx = dginit;
189: sty = zero;
190: fy = finit;
191: dgy = dginit;
192:
193: neP->nfev = 0;
194: for (i=0; i< neP->maxfev; i++) {
195: /* Set min and max steps to correspond to the interval of uncertainty */
196: if (neP->bracket) {
197: neP->stepmin = TaoMin(stx,sty);
198: neP->stepmax = TaoMax(stx,sty);
199: } else {
200: neP->stepmin = stx;
201: neP->stepmax = *step + xtrapf * (*step - stx);
202: }
204: /* Force the step to be within the bounds */
205: *step = TaoMax(*step,neP->stepmin);
206: *step = TaoMin(*step,neP->stepmax);
208: /* If an unusual termination is to occur, then let step be the lowest
209: point obtained thus far */
210: if (((neP->bracket) && (*step <= neP->stepmin || *step >= neP->stepmax)) ||
211: ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol * neP->stepmax)) ||
212: (neP->nfev >= neP->maxfev - 1) || (neP->infoc == 0))
213: *step = stx;
214:
215: info = PetscLogInfo((tao,"TaoApply_LineSearch: Try step size: %6.4e\n",*step)); CHKERRQ(info);
216: #if defined(PETSC_USE_COMPLEX)
217: cstep = *step;
218: info = X->Waxpby(cstep,S,1.0,W);CHKERRQ(info);
219: #else
220: info = X->Waxpby(*step,S,1.0,W);CHKERRQ(info); /* X = W + step*S */
221: #endif
222: /* This if statement added by S. Benson. -- Project the solution, if necessary, to keep feasible */
223: if (*step >= bstepmin1){
224: info=X->Median(XL,X,XU);CHKERRQ(info);
225: }
226: info = TaoComputeMeritFunctionGradient(tao,X,f,G);CHKERRQ(info);
227: neP->nfev++;
228: #if defined(PETSC_USE_COMPLEX)
229: info = G->Dot(S,&cdg);CHKERRQ(info); dg = TaoReal(cdg);
230: return 1;
231: #else
232: // info = G->Dot(S,&dg);CHKERRQ(info); /* dg = G^T S */
233: /* Next 3 lines added by S. Benson. to compute an estimate of decrease in the actual direction and curvature */
234: info = G->Dot(W,&dg1);CHKERRQ(info); /* dg = G^T S */
235: info = G->Dot(X,&dg2);CHKERRQ(info); /* dg = G^T S */
236: dg=dg2-dg1;
237: #endif
238: ftest1 = finit + (*step) * dgtest;
240: /* This test added by S. Benson. -- Be satisfied with Armijo condition if the Projection X->Median(XL,X,XU) changed X. */
241: ftest2 = finit + (*step) * dgtest * neP->ftol;
242:
244: /* Convergence testing */
245: if (((neP->bracket) && (*step <= neP->stepmin||*step >= neP->stepmax)) || (!neP->infoc)) {
246: info = PetscLogInfo((tao,"TaoApply_LineSearch:Rounding errors may prevent further progress. May not be a step satisfying\n")); CHKERRQ(info);
247: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:sufficient decrease and curvature conditions. Tolerances may be too small.\n")); CHKERRQ(info);
248: *info2 = 6; break;
249: }
250: if ((*step == neP->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
251: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Step is at the upper bound, stepmax (%g)\n",neP->stepmax)); CHKERRQ(info);
252: *info2 = 5; break;
253: }
254: if ((*step == neP->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
255: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Step is at the lower bound, stepmin (%g)\n",neP->stepmin)); CHKERRQ(info);
256: *info2 = 4; break;
257: }
258: if (neP->nfev >= neP->maxfev) {
259: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Number of line search function evals (%d) > maximum (%d)\n",neP->nfev,neP->maxfev)); CHKERRQ(info);
260: *info2 = 3; break;
261: }
262: if ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol*neP->stepmax)){
263: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Relative width of interval of uncertainty is at most rtol (%g)\n",neP->rtol)); CHKERRQ(info);
264: *info2 = 2; break;
265: }
266: if ((*f <= ftest1) && (TaoAbsDouble(dg) <= neP->gtol*(-dginit))) {
267: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search success: Sufficient decrease and directional deriv conditions hold, %6.4e, %6.4e\n",TaoAbsDouble(dg), -dginit )); CHKERRQ(info);
268: *info2 = 0; break;
269: }
271: /* This test added by S. Benson. -- Be satisfied with Armijo condition if the Projection X->Median(XL,X,XU) changed X. */
272: if ( (*f <= ftest2) && ((*step) < bstepmin2) ) {
273: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search success: Sufficient decrease new bounds apply\n")); CHKERRQ(info);
274: *info2 = 0; break;
275: }
277: /* In the first stage, we seek a step for which the modified function
278: has a nonpositive value and nonnegative derivative */
279: if ((stage1) && (*f <= ftest1) && (dg >= dginit * TaoMin(neP->ftol, neP->gtol)))
280: stage1 = 0;
282: /* A modified function is used to predict the step only if we
283: have not obtained a step for which the modified function has a
284: nonpositive function value and nonnegative derivative, and if a
285: lower function value has been obtained but the decrease is not
286: sufficient */
288: if ((stage1) && (*f <= fx) && (*f > ftest1)) {
289: fm = *f - *step * dgtest; /* Define modified function */
290: fxm = fx - stx * dgtest; /* and derivatives */
291: fym = fy - sty * dgtest;
292: dgm = dg - dgtest;
293: dgxm = dgx - dgtest;
294: dgym = dgy - dgtest;
296: /* Update the interval of uncertainty and compute the new step */
297: info = TaoStep_LineSearch(tao,&stx,&fxm,&dgxm,&sty,&fym,&dgym,step,&fm,&dgm);CHKERRQ(info);
299: fx = fxm + stx * dgtest; /* Reset the function and */
300: fy = fym + sty * dgtest; /* gradient values */
301: dgx = dgxm + dgtest;
302: dgy = dgym + dgtest;
303: } else {
304: /* Update the interval of uncertainty and compute the new step */
305: info = TaoStep_LineSearch(tao,&stx,&fx,&dgx,&sty,&fy,&dgy,step,f,&dg);CHKERRQ(info);
306: }
308: /* Force a sufficient decrease in the interval of uncertainty */
309: if (neP->bracket) {
310: if (TaoAbsDouble(sty - stx) >= p66 * width1) *step = stx + p5*(sty - stx);
311: width1 = width;
312: width = TaoAbsDouble(sty - stx);
313: }
314: }
316: /* Finish computations */
317: info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:%d function evals in line search, step = %10.4e\n",neP->nfev,*step)); CHKERRQ(info);
318: /*
319: info = G->Norm2(dginitt);CHKERRQ(info);
320: */
321: TaoFunctionReturn(0);
322: }
324: /* ---------------------------------------------------------- */
328: /*@C
329: TaoCreateMoreThuenteBoundLineSearch - Create a line search
331: Input Parameters:
332: + tao - TAO_SOLVER context
333: . fftol - the sufficient descent parameter , greater than 0.
334: - ggtol - the curvature tolerance, greater than 0, less than 1.
337: Note:
338: If either fftol or ggtol is 0, default parameters will be used.
340: Note:
341: This algorithm is is a modification of the algorithm by More' and Thuente.
342: The modifications concern bounds. This algorithm steps in the direction
343: passed into this routine. This point get projected back into the feasible set.
344: It tries to satisfy the Wolfe conditions, but in the context of bound constrained
345: optimization, there may not be a point in the piecewise linear
346: path that satisfies the Wolfe conditions. When the active set
347: is changing, decrease in the objective function may be sufficient
348: to terminate this line search.
350: Note:
351: This routine is used within the following TAO bound constrained
352: minimization solvers: Newton trust region (tao_tron) and limited memory variable
353: metric (tao_blmvm).
355: Level: developer
357: .keywords: TAO_SOLVER, linesearch
358: @*/
359: int TaoCreateMoreThuenteBoundLineSearch(TAO_SOLVER tao, double fftol, double ggtol)
360: {
361: int info;
362: TAO_LINESEARCH *neP;
364: TaoFunctionBegin;
366: info = TaoNew(TAO_LINESEARCH,&neP); CHKERRQ(info);
367: info = PetscLogObjectMemory(tao,sizeof(TAO_LINESEARCH)); CHKERRQ(info);
369: if (fftol<=0) fftol=0.001;
370: if (ggtol<=0) ggtol=0.9;
372: neP->ftol = fftol;
373: neP->rtol = 1.0e-10;
374: neP->gtol = ggtol;
375: neP->stepmin = 1.0e-30;
376: neP->stepmax = 1.0e+20;
377: neP->nfev = 0;
378: neP->bracket = 0;
379: neP->infoc = 1;
380: neP->maxfev = 30;
382: info = TaoSetLineSearch(tao,0,
383: TaoSetOptions_BoundLineSearch,
384: TaoApply_BoundLineSearch,
385: TaoView_LineSearch,
386: TaoDestroy_LineSearch,
387: (void *) neP);CHKERRQ(info);
389: TaoFunctionReturn(0);
390: }
395: /*
398: */