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: }