Actual source code: tao_util.c

  1: /*$Id: tao_util.c 1.84 05/05/10 15:21:38-05:00 sarich@zorak.(none) $*/

  3: #include "src/tao_impl.h"       /*I   "tao_solver.h"   I*/


  8: /*@C
  9:    TaoVecViewMonitor - Monitors progress of the TAO_SOLVER solvers by calling 
 10:    VecView() for the approximate solution at each iteration.

 12:    Collective on TAO_SOLVER

 14:    Input Parameters:
 15: +  tao - the TAO_SOLVER context
 16: -  dummy - either a viewer or TAO_NULL

 18:    Level: advanced

 20: .keywords: TAO_SOLVER, vector, monitor, view

 22: .seealso: TaoSetMonitor(), TaoDefaultMonitor(), VecView()
 23: @*/
 24: int TaoVecViewMonitor(TAO_SOLVER tao,void *dummy)
 25: {
 26:   int    info;
 27:   TaoVec *xx;

 29:   TaoFunctionBegin;
 30:   info = TaoGetSolution(tao,&xx);CHKERRQ(info);
 31:   info = xx->View();CHKERRQ(info);
 32:   TaoFunctionReturn(0);
 33: }

 37: /*@C
 38:    TaoVecViewMonitorUpdate - Monitors progress of the TAO_SOLVER solvers by calling 
 39:    VecView() for the UPDATE to the solution at each iteration.

 41:    Collective on TAO_SOLVER

 43:    Input Parameters:
 44: +  tao - the TAO_SOLVER context
 45: -  dummy - either a viewer or TAO_NULL

 47:    Level: advanced

 49: .keywords: TAO_SOLVER, vector, monitor, view

 51: .seealso: TaoSetMonitor(), TaoDefaultMonitor(), VecView()
 52: @*/
 53: int TaoVecViewMonitorUpdate(TAO_SOLVER tao,void *dummy)
 54: {
 55:   int    info;
 56:   TaoVec *xx;

 58:   TaoFunctionBegin;
 59:   info = TaoGetStepDirectionVector(tao,&xx);CHKERRQ(info);
 60:   if (xx){  info = xx->View();CHKERRQ(info); }
 61:   TaoFunctionReturn(0);
 62: }

 66: /*@C
 67:    TaoDefaultMonitor - Default routine for monitoring progress of the 
 68:    TAO_SOLVER solvers (default).

 70:    Collective on TAO_SOLVER

 72:    Input Parameters:
 73: +  tao - the TAO_SOLVER context
 74: -  dummy - unused context

 76:    Notes:
 77:    The routine prints the function value and gradient norm at each iteration.

 79:    Level: advanced

 81: .keywords: TAO_SOLVER, default, monitor, norm

 83: .seealso: TaoSetMonitor(), TaoVecViewMonitor()
 84: @*/
 85: int TaoDefaultMonitor(TAO_SOLVER tao,void *dummy)
 86: {
 87:   int info;
 88:   int its;
 89:   double fct,gnorm;

 91:   TaoFunctionBegin;
 92:   its=tao->iter;
 93:   fct=tao->fc;
 94:   gnorm=tao->norm;
 95:   info=TaoPrintInt(tao,"iter = %d,",its); CHKERRQ(info);
 96:   info=TaoPrintDouble(tao," Function value: %12.10e,",fct); CHKERRQ(info);
 97:   info=TaoPrintDouble(tao,"  Residual: %12.10e \n",gnorm);CHKERRQ(info);
 98:   TaoFunctionReturn(0);
 99: }

103: /*
104:      Default (short) TAO_SOLVER Monitor, same as TaoDefaultMonitor() except
105:   it prints fewer digits of the residual as the residual gets smaller.
106:   This is because the later digits are meaningless and are often 
107:   different on different machines; by using this routine different 
108:   machines will usually generate the same output.
109: */
110: int TaoDefaultSMonitor(TAO_SOLVER tao,void *dummy)
111: {
112:   int info;
113:   int its;
114:   double  fct,gnorm;

116:   TaoFunctionBegin;
117:   its=tao->iter;
118:   fct=tao->fc;
119:   gnorm=tao->norm;

121:   if (gnorm > 1.e-6) {
122:     info=TaoPrintInt(tao,"iter = %d,",its); CHKERRQ(info);
123:     info=TaoPrintDouble(tao," Function value %g,",fct); CHKERRQ(info);
124:     info=TaoPrintDouble(tao," Residual: %7.6f \n",gnorm);CHKERRQ(info);
125:   } else if (gnorm > 1.e-11) {
126:     info=TaoPrintInt(tao,"iter = %d,",its); CHKERRQ(info);
127:     info=TaoPrintDouble(tao," Function value %g,",fct); CHKERRQ(info);
128:     info=TaoPrintStatement(tao," Residual: < 1.0e-6 \n");
129:   } else {
130:     info=TaoPrintInt(tao,"iter = %d,",its); CHKERRQ(info);
131:     info=TaoPrintDouble(tao," Function value %g,",fct); CHKERRQ(info);
132:     info=TaoPrintStatement(tao," Residual: < 1.0e-11 \n");
133:   }
134:   TaoFunctionReturn(0);
135: }


138: /* ---------------------------------------------------------- */
141: /*@ 
142:    TaoConverged_MaxIts - Determines whether the solver has reached maximum number
143:    of iterations. 

145:    Collective on TAO_SOLVER

147:    Input Parameters:
148: +  tao - the TAO_SOLVER context
149: -  dummy - unused dummy context

151:    Level: developer

153: .seealso: TaoSetTolerances(),TaoGetTerminationReason(),TaoSetTerminationReason()
154: @*/
155: int TaoConverged_MaxIts(TAO_SOLVER tao,void *dummy)
156: {
157:   int info,its, maxits=tao->max_its;
158:   TaoTerminateReason reason=TAO_CONTINUE_ITERATING;

160:   TaoFunctionBegin;
161:   info = TaoGetSolutionStatus(tao,&its,0,0,0,0,0);CHKERRQ(info);
162:   if (its >= maxits) {
163:     info = PetscLogInfo((tao,"TaoConverged_Default: Exceeded maximum number of iterations: %d > %d\n",its,maxits)); CHKERRQ(info);
164:     reason = TAO_DIVERGED_MAXITS;
165:     info = TaoSetTerminationReason(tao,reason); CHKERRQ(info);
166:   }
167:   TaoFunctionReturn(0);
168: }

170: /* ---------------------------------------------------------- */
173: /*@ 
174:    TaoConverged_Default - Determines whether the solver should continue iterating
175:    or terminate. 

177:    Collective on TAO_SOLVER

179:    Input Parameters:
180: +  tao - the TAO_SOLVER context
181: -  dummy - unused dummy context

183:    Output Parameter:
184: .  reason - for terminating

186:    Notes:
187:    This routine checks the residual in the optimality conditions, the 
188:    relative residual in the optimity conditions, the number of function
189:    evaluations, and the function value to test convergence.  Some
190:    solvers may use different convergence routines.

192:    Level: developer

194: .seealso: TaoSetTolerances(),TaoGetTerminationReason(),TaoSetTerminationReason()
195: @*/
196: int TaoConverged_Default(TAO_SOLVER tao,void *dummy)
197: {
198:   int info,its;
199:   int nfuncs=tao->nfuncs+tao->ngrads, max_funcs=tao->max_funcs;
200:   double gnorm, gnorm0=tao->norm0;
201:   double f, trtol=tao->trtol,trradius=tao->step;
202:   double gatol,grtol,gttol,fatol,frtol,catol,crtol;
203:   double fmin=tao->fmin, cnorm, cnorm0=tao->cnorm0;
204:   double gnorm2;
205:   TaoTerminateReason reason=TAO_CONTINUE_ITERATING;

207:   TaoFunctionBegin;
208:   info = TaoGetSolutionStatus(tao,&its,&f,&gnorm,&cnorm,&trradius,&reason);
209:   info = TaoGetTolerances(tao,&fatol,&frtol,&catol,&crtol);CHKERRQ(info);
210:   info = TaoGetGradientTolerances(tao,&gatol,&grtol,&gttol);CHKERRQ(info);
211:   gnorm2=gnorm*gnorm;

213:   if (f != f ) {
214:     info = PetscLogInfo((tao,"TaoConverged_Default: Failed to converged, function is NaN\n")); CHKERRQ(info);
215:     reason = TAO_DIVERGED_NAN;
216:   } else if (f <= fmin && cnorm <=catol) {
217:     info = PetscLogInfo((tao,"TaoConverged_Default: Converged due to function value %g < minimum function value %g\n", f,fmin)); CHKERRQ(info);
218:     reason = TAO_CONVERGED_MINF;
219:   } else if (gnorm2 <= fatol && cnorm <=catol) {
220:     info = PetscLogInfo((tao,"TaoConverged_Default: Converged due to residual norm %g < %g\n",gnorm2,fatol)); CHKERRQ(info);
221:     reason = TAO_CONVERGED_ATOL;
222:   } else if (gnorm2 / TaoAbsScalar(f+1.0e-10)<= frtol && cnorm/TaoMax(cnorm0,1.0) <= crtol) {
223:     info = PetscLogInfo((tao,"TaoConverged_Default: Converged due to relative residual norm %g < %g\n",gnorm2/TaoAbsScalar(f+1.0e-10),frtol)); CHKERRQ(info);
224:     reason = TAO_CONVERGED_RTOL;
225:   } else if (gnorm<= gatol && cnorm <=catol) {
226:     info = PetscLogInfo((tao,"TaoConverged_Default: Converged due to residual norm %g < %g\n",gnorm,gatol)); CHKERRQ(info);
227:     reason = TAO_CONVERGED_ATOL;
228:   } else if ( f!=0 && TaoAbsScalar(gnorm/f) <= grtol && cnorm <= crtol) {
229:     info = PetscLogInfo((tao,"TaoConverged_Default: Converged due to residual norm %g < |%g| %g\n",gnorm,f,grtol)); CHKERRQ(info);
230:     reason = TAO_CONVERGED_ATOL;
231:   } else if (gnorm/gnorm0 <= gttol && cnorm <= crtol) {
232:     info = PetscLogInfo((tao,"TaoConverged_Default: Converged due to relative residual norm %g < %g\n",gnorm/gnorm0,gttol)); CHKERRQ(info);
233:     reason = TAO_CONVERGED_RTOL;
234:   } else if (nfuncs > max_funcs){
235:     info = PetscLogInfo((tao,"TaoConverged_Default: Exceeded maximum number of function evaluations: %d > %d\n", nfuncs,max_funcs)); CHKERRQ(info);
236:     reason = TAO_DIVERGED_MAXFCN;
237:   } else if ( tao->lsflag != 0 ){
238:     info = PetscLogInfo((tao,"TaoConverged_Default: Tao Line Search failure.\n")); CHKERRQ(info);
239:     reason = TAO_DIVERGED_LS_FAILURE;
240:   } else if (trradius < trtol && its > 0){
241:     info = PetscLogInfo((tao,"TaoConverged_Default: Trust region/step size too small: %g < %g\n", trradius,trtol)); CHKERRQ(info);
242:     reason = TAO_CONVERGED_TRTOL;
243:   } else {
244:     reason = TAO_CONTINUE_ITERATING;
245:   }
246:   info = TaoSetTerminationReason(tao,reason); CHKERRQ(info);
247:   info=TaoConverged_MaxIts(tao,0); CHKERRQ(info);
248:   TaoFunctionReturn(0);
249: }

251: /* ---------------------------------------------------------- */
254: /*@C 
255:    TaoCheckFG - Check if the function and gradient vectors have
256:    been set properly and are compatible. 

258:    Collective on TAO_SOLVER

260:    Input Parameters:
261: .  tao - the TAO_SOLVER context

263:    Level: developer

265: .seealso: TaoCheckFGH()
266: @*/
267: int TaoCheckFG(TAO_SOLVER tao)
268: {
269:   TaoVec *xx,*gg;
270:   TaoTruth flag;
271:   int info;

273:   TaoFunctionBegin;
274:   TaoValidHeaderSpecific(tao,TAO_COOKIE,1);

276:   info=TaoGetSolution(tao,&xx);CHKERRQ(info);
277:   info=TaoGetGradient(tao,&gg);CHKERRQ(info);

279:   info=xx->Compatible(gg,&flag);CHKERRQ(info);
280:   if (flag==TAO_FALSE){
281:     SETERRQ(1,"Gradient and Variable vectors must have identical structure.");
282:   }

284:   TaoFunctionReturn(0);
285: }

289: /*@C 
290:    TaoCheckFGH - Check if the function, gradient vector, and
291:    Hessian matrix have been set properly and are compatible. 

293:    Collective on TAO_SOLVER

295:    Input Parameters:
296: .  tao - the TAO_SOLVER context

298:    Level: developer

300: .seealso: TaoCheckFG()
301: @*/
302: int TaoCheckFGH(TAO_SOLVER tao)
303: {
304:   int info;
305:   TaoMat *HH;
306:   TaoFunctionBegin;
307: 
308:   info=TaoCheckFG(tao);CHKERRQ(info);
309:   info = TaoGetHessian(tao,&HH);CHKERRQ(info);
310:   if (!HH) {
311:     SETERRQ(1,"Must Provide Hessian Matrix");
312:   }

314:   TaoFunctionReturn(0);
315: }

319: /*@C 
320:    TaoCheckConstraints - Check if the nonlinear constraints have
321:    been set and if the data structures are compatible.

323:    Collective on TAO_SOLVER

325:    Input Parameters:
326: .  tao - the TAO_SOLVER context

328:    Level: developer

330: .seealso: TaoCheckFG()
331: @*/
332: int TaoCheckConstraints(TAO_SOLVER tao)
333: {
334:   TaoVec *solu, *cons;
335:   TaoMat *jac;
336:   TaoTruth flag;
337:   int info;

339:   TaoFunctionBegin;
340:   TaoValidHeaderSpecific(tao,TAO_COOKIE,1);

342:   info = TaoGetSolution(tao, &solu); CHKERRQ(info);
343:   info = TaoGetConstraints(tao, &cons); CHKERRQ(info);
344:   info = TaoGetJacobian(tao, &jac); CHKERRQ(info);
345: 
346:   info = jac->Compatible(solu,cons,&flag);CHKERRQ(info);
347:   if (flag == TAO_FALSE){
348:     SETERRQ(1,"Jacobian matrix not consistent with Variable Vector or Constraint Vector");
349:   }

351:   TaoFunctionReturn(0);
352: }
353: 
356: /*@C 
357:    TaoCheckBounds - Check if the variable bounds have been
358:    set and the data structures are compatible.

360:    Collective on TAO_SOLVER

362:    Input Parameters:
363: .  tao - the TAO_SOLVER context

365:    Level: developer

367: .seealso: TaoCheckFG()
368: @*/
369: int TaoCheckBounds(TAO_SOLVER tao)
370: {
371:   int info;
372:   TaoTruth flag;
373:   TaoVec *xx,*xxll,*xxuu;
374:   //  TaoIndexSet *SS;

376:   info=TaoGetSolution(tao,&xx);CHKERRQ(info);
377:   info=TaoGetVariableBounds(tao,&xxll,&xxuu);CHKERRQ(info);

379:   info=xx->Compatible(xxll,&flag);CHKERRQ(info);
380:   if (flag==TAO_FALSE){
381:     SETERRQ(1,"Vector of lower bounds not Compatible with Variable Vector");
382:   }

384:   info=xx->Compatible(xxuu,&flag);CHKERRQ(info);
385:   if (flag==TAO_FALSE){
386:     SETERRQ(1,"Vector of upper bounds not Compatible with Variable vector");
387:   }
388:   /*
389:   info = xx->CreateIndexSet(&SS);CHKERRQ(info);

391:   info = SS->whichLessThan(xxuu,xxll);CHKERRQ(info);
392:   info = SS->GetSize(&n);CHKERRQ(info);
393:   if (n>0)
394:     SETERRQ(1,"Error: Upper bound not greater than lower bound.");
395:   */
396:   /*
397:   info = S->whichLessThanEqualTo(XU,TAO_NINFINITY);CHKERRQ(info);
398:   info = S->GetSize(&n);CHKERRQ(info);
399:   if (n>0)
400:     SETERRQ(1,0,"Upper bound cannot be negative infinity .");

402:   info = S->whichGreaterThanEqualTo(XL,TAO_INFINITY);CHKERRQ(info);
403:   info = S->GetSize(&n);CHKERRQ(info);
404:   if (n>0)
405:     SETERRQ(1,0,"Error: Lower bound cannot be infinity.");
406:   */
407:   /*
408:   info = TaoIndexSetDestroy(SS);CHKERRQ(info);
409:   */
410:   TaoFunctionReturn(0);
411: }


414: static int TaoDefaultMeritFunction(TAO_SOLVER,TaoVec*,double*,void*);
415: static int TaoDefaultMeritFunctionGradient(TAO_SOLVER,TaoVec*,double*,TaoVec*,void*);
416: static int TaoDefaultMeritGradient(TAO_SOLVER,TaoVec*,TaoVec*,void*);


421: static int TaoDefaultMeritFunction(TAO_SOLVER tao, TaoVec *xx, double *f, void *ctx)
422: {
423:   int info;

425:   TaoFunctionBegin;
426:   info=TaoComputeFunction(tao,xx,f);CHKERRQ(info);
427:   TaoFunctionReturn(0);
428: }

432: static int TaoDefaultMeritFunctionGradient(TAO_SOLVER tao, TaoVec *xx, double *f, TaoVec *gg, void *ctx)
433: {
434:   int info;

436:   TaoFunctionBegin;
437:   info=TaoComputeFunctionGradient(tao,xx,f,gg);CHKERRQ(info);
438:   TaoFunctionReturn(0);
439: }

443: static int TaoDefaultMeritGradient(TAO_SOLVER tao, TaoVec *xx, TaoVec *gg, void *ctx)
444: {
445:   int info;

447:   TaoFunctionBegin;
448:   info=TaoComputeGradient(tao,xx,gg);CHKERRQ(info);
449:   TaoFunctionReturn(0);
450: }



456: /*@
457:    TaoSetDefaultMeritFunction - Set the merit function equal to the
458:    objective function

460:    Collective on TAO_SOLVER

462:    Input Parameters:
463: .  tao - the TAO_SOLVER context

465:    Level: developer

467: @*/
468: int TaoSetDefaultMeritFunction(TAO_SOLVER tao)
469: {
470:   int info;
471:   TaoFunctionBegin;
472:   info=TaoSetMeritFunction(tao,TaoDefaultMeritFunction,TaoDefaultMeritFunctionGradient,TaoDefaultMeritGradient,0,0);
473:   CHKERRQ(info);
474:   TaoFunctionReturn(0);
475: }
478: /* @C 
479:    TaoSetQuadraticMeritFunction - Set the merit function to a quadratic
480:    objective function 1/2 x^T A x + b^T x + c

482:    Collective on TAO_SOLVER

484:    Input Parameters:
485: +  tao - the TAO_SOLVER context
486: .  A - The matrix term of the quadratic function
487: .  b - The linear term of the quadratic function
488: -  c - The constant term of the quadratic funtion

490:    Level: developer

492: @ */
493: int TaoSetQuadraticMeritFunction(TAO_SOLVER tao, TaoMat *A, TaoVec*b, double c)
494: {
495:   TaoFunctionBegin;
496:   SETERRQ(1,"TAO ERROR: Not implemented");
497:   //  TaoFunctionReturn(0);
498: }