Actual source code: boundproj.c

  2: #include "boundproj.h"    /*I "tao_solver.h" I*/

  6: static int TaoDestroy_LineSearch(TAO_SOLVER tao,void *linectx)
  7: {
  8:   int  info;
  9:   TAO_LINESEARCH2 *ctx = (TAO_LINESEARCH2 *)linectx;

 11:   TaoFunctionBegin;
 12:   if (ctx->setupcalled==1){
 13:     info=TaoVecDestroy(ctx->W2);CHKERRQ(info);
 14:   }
 15:   info = TaoFree(ctx);CHKERRQ(info);
 16:   TaoFunctionReturn(0);
 17: }

 21: static int TaoSetOptions_LineSearch(TAO_SOLVER tao,void *linectx)
 22: {
 23:   TAO_LINESEARCH2 *ctx = (TAO_LINESEARCH2 *)linectx;
 24:   int            info;

 26:   TaoFunctionBegin;
 27:   info = TaoOptionsHead("Projected line search options");CHKERRQ(info);
 28:   info = TaoOptionInt("-tao_ls_maxfev","max function evals in line search","",ctx->maxfev,&ctx->maxfev,0);CHKERRQ(info);
 29:   info = TaoOptionDouble("-tao_ls_ftol","tol for sufficient decrease","",ctx->ftol,&ctx->ftol,0);CHKERRQ(info);
 30:   info = TaoOptionDouble("-tao_ls_stepmin","lower bound for step","",ctx->stepmin,&ctx->stepmin,0);CHKERRQ(info);
 31:   info = TaoOptionDouble("-tao_ls_stepmax","upper bound for step","",ctx->stepmax,&ctx->stepmax,0);CHKERRQ(info);
 32:   info = TaoOptionsTail();CHKERRQ(info);

 34:   TaoFunctionReturn(0);
 35: }

 39: static int TaoView_LineSearch(TAO_SOLVER tao,void *ctx)
 40: {
 41:   TAO_LINESEARCH2 *ls = (TAO_LINESEARCH2 *)ctx;
 42:   int            info;

 44:   TaoFunctionBegin;
 45:   info = TaoPrintInt(tao,"  Line search: maxf=%d,",ls->maxfev);CHKERRQ(info);
 46:   info = TaoPrintDouble(tao," ftol=%g\n",ls->ftol);CHKERRQ(info);
 47:   TaoFunctionReturn(0);
 48: }

 52: /* @ TaoApply_LineSearch - This routine performs a line search. It
 53:    backtracks until the Armijo conditions are satisfied.

 55:    Input Parameters:
 56: +  tao - TAO_SOLVER context
 57: .  X - current iterate (on output X contains new iterate, X + step*S)
 58: .  S - search direction
 59: .  f - objective function evaluated at X
 60: .  G - gradient evaluated at X
 61: .  W - work vector
 62: .  gdx - inner product of gradient and the direction of the first linear manifold being searched
 63: -  step - initial estimate of step length

 65:    Output parameters:
 66: +  f - objective function evaluated at new iterate, X + step*S
 67: .  G - gradient evaluated at new iterate, X + step*S
 68: .  X - new iterate
 69: -  step - final step length

 71:    Info is set to one of:
 72: .   0 - the line search succeeds; the sufficient decrease
 73:    condition and the directional derivative condition hold

 75:    negative number if an input parameter is invalid
 76: .   -1 -  step < 0 
 77: .   -2 -  ftol < 0 
 78: -   -7 -  maxfev < 0

 80:    positive number > 1 if the line search otherwise terminates
 81: +    2 -  Relative width of the interval of uncertainty is 
 82:          at most rtol.
 83: .    3 -  Maximum number of function evaluations (maxfev) has 
 84:          been reached.
 85: .    4 -  Step is at the lower bound, stepmin.
 86: .    5 -  Step is at the upper bound, stepmax.
 87: .    6 -  Rounding errors may prevent further progress. 
 88:          There may not be a step that satisfies the 
 89:          sufficient decrease and curvature conditions.  
 90:          Tolerances may be too small.
 91: +    7 -  Search direction is not a descent direction.

 93:    Notes:
 94:    This routine is used within the TAO_NLS method.
 95: @ */
 96: static int TaoApply_LineSearch(TAO_SOLVER tao,TaoVec* X,TaoVec* G,TaoVec* S,TaoVec* Gold,double *f,double *f_full,
 97:                         double *step,TaoInt *info2,void*ctx)
 98: {
 99:   TAO_LINESEARCH2 *neP = (TAO_LINESEARCH2 *) ctx;
100:   int       info;
101:   TaoInt i;
102:   double zero=0.0;
103:   double finit,gdx,dginit,actred,prered,rho,d1=0,d2=0,d3=0;
104:   TaoVec* XL, *XU, *Xold=neP->W2;
105:   TaoTruth flag;

107:   TaoFunctionBegin;
108:   /* neP->stepmin - lower bound for step */
109:   /* neP->stepmax - upper bound for step */
110:   /* neP->rtol           - relative tolerance for an acceptable step */
111:   /* neP->ftol           - tolerance for sufficient decrease condition */
112:   /* neP->gtol           - tolerance for curvature condition */
113:   /* neP->nfev           - number of function evaluations */
114:   /* neP->maxfev  - maximum number of function evaluations */

116:   /* Check input parameters for errors */

118:   *info2 = 0;
119:     /* After 2 failures, Check that search direction is a descent direction */
120:     /* This test is not really sufficient.  */
121:   if (neP->setupcalled){
122:     info=X->Compatible(neP->W2,&flag); CHKERRQ(info);
123:     if (flag==TAO_FALSE){
124:       info=TaoVecDestroy(neP->W2); CHKERRQ(info);neP->W2=0;
125:       neP->setupcalled=0;
126:     }
127:   }
128:   if (neP->setupcalled==0){
129:     info=X->Clone(&neP->W2); CHKERRQ(info);
130:     Xold=neP->W2;
131:     neP->setupcalled=1;
132:   }
133:   info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
134:   info = Gold->ScaleCopyFrom(-1.0,S); CHKERRQ(info);
135:   info = Gold->BoundGradientProjection(S,XL,X,XU); CHKERRQ(info);
136:   info = Gold->Dot(G,&dginit);CHKERRQ(info);  /* dginit = G^T S */
137:   dginit*=-1;
138:   gdx=dginit;
139:   if (dginit >= zero) {
140:     info = G->CopyFrom(Gold); CHKERRQ(info);
141:     info = X->CopyFrom(Xold); CHKERRQ(info);
142:     info = PetscInfo(tao,"TaoApply_LineSearch2:Search direction not a descent direction\n"); CHKERRQ(info);
143:     *info2 = 7; TaoFunctionReturn(0);
144:   }

146:   info = Xold->CopyFrom(X); CHKERRQ(info);
147:   info = Gold->CopyFrom(G); CHKERRQ(info);
148:   info=TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
149:   if (*step < zero) {
150:     info = PetscInfo1(tao,"TaoApply_LineSearch:Line search error: step (%g) < 0\n",*step); CHKERRQ(info);
151:     *info2 = -1; TaoFunctionReturn(0);
152:   } else if (neP->ftol < zero) {
153:     info = PetscInfo1(tao,"TaoApply_LineSearch:Line search error: ftol (%g) < 0\n",neP->ftol); CHKERRQ(info);
154:     *info2 = -2; TaoFunctionReturn(0);
155:   } else if (neP->maxfev < zero) {
156:     info = PetscInfo1(tao,"TaoApply_LineSearch:Line search error: maxfev (%d) < 0\n",neP->maxfev); CHKERRQ(info);
157:     *info2 = -7; TaoFunctionReturn(0);
158:   }

160:   /* Initialization */
161:   neP->nfev = 0;
162:   finit = *f;
163:   tao->new_search=TAO_TRUE;
164:   for (i=0; i< neP->maxfev; i++) {
165:     
166:     /* Force the step to be within the bounds */
167:     *step = TaoMax(*step,neP->stepmin);
168:     *step = TaoMin(*step,neP->stepmax);
169:     
170:     if (0==i){
171:       info = X->Axpy(*step,S);CHKERRQ(info);
172:       info = X->Median(XL,X,XU);CHKERRQ(info);
173:       tao->current_step=*step;
174:       info = TaoComputeMeritFunctionGradient(tao,X,f,G);
175:       tao->new_search=TAO_FALSE;
176:       *f_full = *f;

178:       neP->nfev++;
179:       if (info==0 && *f==*f){ /* Successful function evaluation */
180:         actred = *f - finit;
181:         rho = actred/( (*step) * (-TaoAbsScalar(gdx)) );
182:         if (actred < 0 && rho > neP->ftol){
183:           break;
184:         } else{
185:           info=Xold->StepBoundInfo(XL,XU,S,&d3,&d2,&d1);CHKERRQ(info);
186:           info = PetscInfo1(tao,"stepmax = %10f\n",d1); CHKERRQ(info);

188:           *step = TaoMin(*step,d1);
189:           info = Gold->Dot(X,&d1); CHKERRQ(info);
190:           info = Gold->Dot(Xold,&prered); CHKERRQ(info);
191:           prered=d1-prered;
192:           rho = actred/(-TaoAbsScalar(prered));
193:           if (actred < 0 && rho > neP->ftol){
194:             break;
195:           } else{
196:                 *step = (*step)/2;
197:           }
198:         }
199:       } else { /* Function could not be evaluated at this point */
200:         *step = (*step)*0.7;
201:       }

203:     } else {
204:       info = X->Waxpby(*step,S,1.0,Xold);CHKERRQ(info);
205:       info = X->Median(XL,X,XU);CHKERRQ(info);
206:       
207:       info = G->Waxpby(-1,Xold,1.0,X);CHKERRQ(info);
208:       info = G->Dot(Gold,&prered); CHKERRQ(info);
209:       tao->current_step=*step;
210:       info = TaoComputeMeritFunctionGradient(tao,X,f,G); CHKERRQ(info);
211:       tao->new_search=TAO_FALSE;
212:       neP->nfev++;
213:       if (info==0 && *f==*f){ /* Successful function evaluation */
214:         actred = *f - finit;
215:         rho = actred/(-TaoAbsScalar(prered));
216:         /* 
217:            If sufficient progress has been obtained, accept the
218:            point.   Prered could be positive or negative.  
219:            Otherwise, backtrack. 
220:         */

222:         if (actred < 0 && rho > neP->ftol){
223:           info = PetscInfo4(tao,"TaoApply_LineSearch: steplength: %g, actual reduction: %g (hopefully positive), Predicted reduction: %g, rho: %g\n",*step,-actred,prered,rho); CHKERRQ(info);
224:           break;
225:         } else {
226:           info = PetscInfo4(tao,"TaoApply_LineSearch: steplength: %g, actual reduction: %g (hopefully positive), Predicted reduction: %g, rho: %g\n",*step,-actred,prered,rho); CHKERRQ(info);
227:           if (i<5){
228:             *step = (*step)/2;
229:           } else {
230:             *step = (*step)/2;
231:           }
232:         }
233:       } else { /* Function could not be evaluated at this point */
234:           info = PetscInfo(tao,"TaoApply_LineSearch: Problem in function evaluation\n"); CHKERRQ(info);
235:         *step = (*step)*0.7;
236:       }
237:     }

239:   }
240:   /* Convergence testing */
241:   
242:   if (*step <= neP->stepmin || *step >= neP->stepmax) {
243:     *info2 = 6;
244:     info = PetscInfo(tao,"TaoApply_LineSearch:Rounding errors may prevent further progress.  May not be a step satisfying\n"); CHKERRQ(info);
245:     info = PetscInfo(tao,"TaoApply_LineSearch:sufficient decrease and curvature conditions. Tolerances may be too small.\n"); CHKERRQ(info);
246:   }
247:   if (*step == neP->stepmax) {
248:     info = PetscInfo1(tao,"TaoApply_LineSearch:Step is at the upper bound, stepmax (%g)\n",neP->stepmax); CHKERRQ(info);
249:     *info2 = 5;
250:   }
251:   if (*step == neP->stepmin) {
252:     info = PetscInfo1(tao,"TaoApply_LineSearch:Step is at the lower bound, stepmin (%g)\n",neP->stepmin); CHKERRQ(info);
253:     *info2 = 4;
254:   }
255:   if (neP->nfev >= neP->maxfev) {
256:     info = PetscInfo2(tao,"TaoApply_LineSearch:Number of line search function evals (%d) > maximum (%d)\n",neP->nfev,neP->maxfev); CHKERRQ(info);
257:     *info2 = 3;
258:   }
259:   if ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol*neP->stepmax)){
260:     info = PetscInfo1(tao,"TaoApply_LineSearch:Relative width of interval of uncertainty is at most rtol (%g)\n",neP->rtol); CHKERRQ(info);
261:     *info2 = 2;
262:   }
263:   /*
264:   if ((*f <= ftest1) && (TaoAbsDouble(dg) <= neP->gtol*(-dginit))) {
265:     info = PetscLogInfo((tao,"TaoApply_LineSearch:Line search success: Sufficient decrease and directional deriv conditions hold\n")); CHKERRQ(info);
266:     *info2 = 1;
267:   }
268:   */
269:   
270:   /* Finish computations */
271:   /*
272:   info = PetscLogInfo((tao,"TaoApply_LineSearch:%d function evals in line search, step = %10.4f\n",neP->nfev,*step)); CHKERRQ(info);
273:   */
274:   TaoFunctionReturn(0);
275: }

279: /*@C
280:    TaoCreateProjectedLineSearch - Create a line search

282:    Input Parameters:
283: .  tao - TAO_SOLVER context

285:    Note:
286:    This routine is used within the following TAO bound constrained 
287:    minimization solvers: TRON (tao_tron) and limited memory variable metric
288:    (tao_blmvm). This line search projects points onto the feasible, bound
289:    constrained region.  It only enforces the Armijo descent condition.

291:    Level: developer

293: .keywords: TAO_SOLVER, linesearch
294: @*/
295: int TaoCreateProjectedLineSearch(TAO_SOLVER tao)
296: {
297:   int info;
298:   TAO_LINESEARCH2 *neP;

300:   TaoFunctionBegin;

302:   info = TaoNew(TAO_LINESEARCH2,&neP); CHKERRQ(info);
303:   info = PetscLogObjectMemory(tao,sizeof(TAO_LINESEARCH2)); CHKERRQ(info);

305:   neP->ftol                  = 0.001;
306:   neP->rtol                  = 0.0;
307:   neP->gtol                  = 0.0;
308:   neP->stepmin                  = 1.0e-30;
309:   neP->stepmax                  = 1.0e+20;
310:   neP->nfev                  = 0; 
311:   neP->bracket                  = 0; 
312:   neP->infoc              = 1;
313:   neP->maxfev                  = 30;
314:   neP->W2                 = TAO_NULL;
315:   neP->setupcalled        = 0;

317:   info = TaoSetLineSearch(tao,0,
318:                           TaoSetOptions_LineSearch,
319:                           TaoApply_LineSearch,
320:                           TaoView_LineSearch,
321:                           TaoDestroy_LineSearch,
322:                           (void *) neP);CHKERRQ(info);

324:   TaoFunctionReturn(0);
325: }