Actual source code: nlmvm.c

      1: /*$Id: s.nls.c 1.106 04/03/03 13:56:23-06:00 benson@rockies.mcs.anl.gov $*/

      3: #include "nlmvm.h"             /*I "tao_solver.h" I*/

      5: static int TaoDestroy_NLMVM(TAO_SOLVER tao, void *solver);
      6: static int TaoView_NLMVM(TAO_SOLVER tao,void* solver);
      7: static int TaoSolve_NLMVM(TAO_SOLVER tao,void *solver);
      8: static int TaoSetOptions_NLMVM(TAO_SOLVER tao, void *solver);
      9: static int TaoSetUp_NLMVM(TAO_SOLVER tao, void *solver);

     11: /*
     12:    Implements Newton's Method with a line search approach
     13:    for solving unconstrained minimization problems.

     15:    Note:
     16:    The line search algorithm is taken from More and Thuente,
     17:    "Line search algorithms with guaranteed sufficient decrease",
     18:    Argonne National Laboratory, Technical Report MCS-P330-1092.
     19: */

     21: class TaoH00Mat: public TaoMat{
     22:  protected:
     23:   TAO_SOLVER tao;
     24:  public:
     25:   TaoH00Mat(TAO_SOLVER);
     26:   ~TaoH00Mat();
     27:   int Solve(TaoVec*, TaoVec*, TaoTruth*);
     28: };

     32: TaoH00Mat::TaoH00Mat(TAO_SOLVER ttao){
     33:   this->tao=ttao;
     34:   return;
     35: }

     39: TaoH00Mat::~TaoH00Mat(){
     40:   return;
     41: }

     45: int TaoH00Mat::Solve(TaoVec* rhs, TaoVec* dx, TaoTruth *success){
     46:   int info;
     47:   TaoMat *HH;
     48:   info=TaoGetHessian(this->tao,&HH);CHKERRQ(info);
     49:   info=TaoLinearSolve(this->tao,HH,rhs,dx,success);CHKERRQ(info);
     50:   return 0;
     51: }

     55: static int TaoSolve_NLMVM(TAO_SOLVER tao,void *solver)
     56: {
     57:   TAO_NLMVM            *neP = (TAO_NLMVM *) solver;
     58:   int                ii,iter=0, info, line=0;
     59:   double             r0,f,dxnorm,gdx, gnorm, two = 2.0,step=0.0;
     60:   TaoVec             *gg, *xx;
     61:   TaoVec             *rrhs=neP->RHS, *ss=neP->S, *ww=neP->W;
     62:   TaoMat             *HH;
     63:   TaoTerminateReason reason;
     64:   TaoTruth           success;
     65:   double             gamma,gamma_factor;
     66:   TaoH00Mat *HH0;
     67:   TaoFunctionBegin;

     69:   info=TaoGetSolution(tao,&xx);CHKERRQ(info);
     70:   info=TaoGetGradient(tao,&gg);CHKERRQ(info);
     71:   info=TaoGetHessian(tao,&HH);CHKERRQ(info);

     73:   HH0=new TaoH00Mat(tao);
     74:   info = neP->M->SetH0( HH0 );CHKERRQ(info);

     76:   info = TaoComputeMeritFunctionGradient(tao,xx,&f,gg);CHKERRQ(info);
     77:   info = gg->Norm2(&gnorm);CHKERRQ(info);         /* gnorm = || gg || */

     79:   gamma       = neP->gamma;
     80:   gamma_factor= neP->gamma_factor;
     81:
     82:   while (1) {
     83:
     84:     info = TaoMonitor(tao,iter++,f,gnorm,0.0,step,&reason);CHKERRQ(info);
     85:     if (reason!=TAO_CONTINUE_ITERATING) break;

     87:     gamma_factor /= two;
     88:     gamma = gamma_factor*(gnorm);

     90:     info = rrhs->Axpby(-1.0,gg,0.0); CHKERRQ(info);

     92:     info = TaoComputeHessian(tao,xx,HH);CHKERRQ(info);

     94:     info = neP->M->reset(); CHKERRQ(info);
     95:     info = neP->M->update(xx,gg); CHKERRQ(info);

     97:     success = TAO_FALSE;
     98:     while (success==TAO_FALSE) {
     99:       info = TaoPreLinearSolve(tao,HH);CHKERRQ(info);
    100:       info = TaoLinearSolve(tao,HH,rrhs,ss,&success);CHKERRQ(info);

    102:       info = ss->Dot(gg,&gdx); CHKERRQ(info);

    104:       if (gdx>=0 /*|| success==TAO_FALSE*/) { /* Modify diagonal of Hessian */
    105:         gamma_factor *= two;
    106:         gamma = gamma_factor*(gnorm);
    107: #if !defined(PETSC_USE_COMPLEX)
    108:         info = PetscLogInfo((tao,"TaoSolve_NLMVM:  modify diagonal (assume same nonzero structure), gamma_factor=%g, gamma=%g\n",
    109:                              gamma_factor,gamma)); CHKERRQ(info);
    110: #else
    111:         info = PetscLogInfo((tao,"TaoSolve_NLMVM:  modify diagonal (assume same nonzero structure), gamma_factor=%g, gamma=%g\n",
    112:                              gamma_factor,PetscReal(gamma))); CHKERRQ(info);
    113: #endif
    114:         info = HH->ShiftDiagonal(gamma);CHKERRQ(info);
    115:         success = TAO_FALSE;

    117:         /* We currently assume that all diagonal elements were allocated in
    118:            original matrix, so that nonzero pattern is same ... should fix this */
    119:       } else {
    120:         //        gamma_factor /= two;
    121:         success = TAO_TRUE;
    122:       }
    123:     }

    125:     /* Line search */
    126:     step=1.0;
    127:     info = TaoLineSearchApply(tao,xx,gg,ss,ww,&f,&step,&line);
    128:     CHKERRQ(info);
    129:     info = gg->Norm2(&gnorm);CHKERRQ(info);

    131:     for (ii=0;ii<neP->lmvmits;ii++){
    132:       /* Test for convergence */
    133:       info=TaoMonitor(tao,iter,f,gnorm,0,step,&reason);CHKERRQ(info);
    134:       if (reason!=TAO_CONTINUE_ITERATING) break;
    135:
    136:       /* Compute step direction */
    137:       info = neP->M->update(xx,gg); CHKERRQ(info);
    138:       info = neP->M->Solve(gg,ss,&success); CHKERRQ(info);
    139:       info = ss->Negate(); CHKERRQ(info);
    140:       info = ss->Dot(gg,&gdx); CHKERRQ(info);
    141:       if (gdx>=0){
    142:         info = ss->ScaleCopyFrom(-1.0,gg); CHKERRQ(info);
    143:         gdx = -gnorm*gnorm;
    144:       }
    145:
    146:       /* Line Search */
    147:       info = ss->Norm2(&dxnorm);  CHKERRQ(info);
    148:       info = TaoGetInitialTrustRegionRadius(tao,&r0); CHKERRQ(info);
    149:       step = TaoMin(1.0,r0/dxnorm);
    150:       info = TaoLineSearchApply(tao,xx,gg,ss,ww,&f,&step,&line);
    151:       CHKERRQ(info);
    152:       info = gg->Norm2(&gnorm);CHKERRQ(info);
    153:     }
    154:
    155:     //    neP->gamma=gamma;
    156:     //.    neP->gamma_factor=gamma_factor;

    158:   }

    160:   delete HH0;
    161:   TaoFunctionReturn(0);
    162: }
    163: /* ---------------------------------------------------------- */
    166: static int TaoSetUp_NLMVM(TAO_SOLVER tao, void *solver)
    167: {
    168:   int        info,kspmaxit,lm,n;
    169:   TAO_NLMVM *ctx = (TAO_NLMVM *)solver;
    170:   TaoVec *xx;
    171:   TaoMat *HH;
    172:   TaoLinearSolver *tls;

    174:   TaoFunctionBegin;
    175:   info = TaoGetSolution(tao,&xx);CHKERRQ(info);

    177:   info=xx->Clone(&ctx->RHS);CHKERRQ(info);
    178:   info=xx->Clone(&ctx->S);CHKERRQ(info);
    179:   info=xx->Clone(&ctx->W);CHKERRQ(info);
    180:   info=xx->Clone(&ctx->G);CHKERRQ(info);

    182:   info = TaoSetLagrangianGradientVector(tao,ctx->G);CHKERRQ(info);
    183:   info = TaoSetStepDirectionVector(tao,ctx->S);CHKERRQ(info);

    185:   info = TaoLMVMGetSize(tao,&lm); CHKERRQ(info);
    186:   ctx->M = new TaoLMVMMat(xx,TaoMin(lm,ctx->lmvmits));

    188:   info = TaoGetHessian(tao,&HH); CHKERRQ(info);
    189:   info = TaoCreateLinearSolver(tao,HH,200,0); CHKERRQ(info);
    190:   info = TaoGetLinearSolver(tao,&tls); CHKERRQ(info);
    191:   info = xx->GetDimension(&n);CHKERRQ(info);
    192:   kspmaxit = ctx->max_kspiter_factor * ((int) sqrt((double)n));
    193:   info = tls->SetTolerances(1e-5,1e-40,1e30,kspmaxit);CHKERRQ(info);

    195:   info = TaoCheckFGH(tao);CHKERRQ(info);

    197:   TaoFunctionReturn(0);
    198: }
    199: /*------------------------------------------------------------*/
    202: static int TaoDestroy_NLMVM(TAO_SOLVER tao, void *solver)
    203: {
    204:   TAO_NLMVM *ctx = (TAO_NLMVM *)solver;
    205:   int     info;

    207:   TaoFunctionBegin;

    209:   info=TaoVecDestroy(ctx->RHS);CHKERRQ(info);
    210:   info=TaoVecDestroy(ctx->S);CHKERRQ(info);
    211:   info=TaoVecDestroy(ctx->W);CHKERRQ(info);
    212:   info=TaoVecDestroy(ctx->G);CHKERRQ(info);
    213:   info=TaoMatDestroy(ctx->M);CHKERRQ(info);
    214:   info = TaoDestroyLinearSolver(tao);CHKERRQ(info);

    216:   TaoFunctionReturn(0);
    217: }
    218: /*------------------------------------------------------------*/
    221: static int TaoSetOptions_NLMVM(TAO_SOLVER tao, void *solver)
    222: {
    223:   TAO_NLMVM    *ctx = (TAO_NLMVM *)solver;
    224:   int        info;
    225:   TaoTruth   flg;

    227:   TaoFunctionBegin;
    228:   info = TaoOptionsHead("Newton line search method for unconstrained optimization");CHKERRQ(info);

    230:   info = TaoOptionDouble("-tao_nls_gamma_factor","damping parameter","",ctx->gamma_factor,&ctx->gamma_factor,&flg);CHKERRQ(info);

    232:   info = TaoOptionInt("-tao_nls_maxkspf","max ksp iterates","",ctx->max_kspiter_factor,&ctx->max_kspiter_factor,&flg);CHKERRQ(info);

    234:   info = TaoOptionInt("-tao_nlmvmits","LMVM updates per iteration","",ctx->lmvmits,&ctx->lmvmits,&flg);CHKERRQ(info);

    236:   info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);
    237:   /*  info = TaoSetLinearSolverOptions(tao);CHKERRQ(info); */

    239:   info = TaoOptionsTail();CHKERRQ(info);

    241:   TaoFunctionReturn(0);
    242: }

    245: /*------------------------------------------------------------*/
    248: static int TaoView_NLMVM(TAO_SOLVER tao,void* solver)
    249: {
    250:   TAO_NLMVM *ctx = (TAO_NLMVM *)solver;
    251:   int     info;

    253:   TaoFunctionBegin;

    255:   info = TaoPrintDouble(tao,"  gamma_f=%g,",ctx->gamma_factor);CHKERRQ(info);
    256:   info = TaoPrintInt(tao," maxkspf=%d\n",ctx->max_kspiter_factor);CHKERRQ(info);

    258:   info = TaoPrintStatement(tao,"  Linear Solver for symmetric matrix: \n");CHKERRQ(info);
    259:   info = TaoLineSearchView(tao);CHKERRQ(info);

    261:   TaoFunctionReturn(0);
    262: }

    264: /* ---------------------------------------------------------- */
    268: int TaoCreate_NLMVM(TAO_SOLVER tao)
    269: {
    270:   TAO_NLMVM         *neP;
    271:   int             info;

    273:   TaoFunctionBegin;

    275:   info = TaoNew(TAO_NLMVM,&neP); CHKERRQ(info);
    276:   info = PetscLogObjectMemory(tao,sizeof(TAO_NLMVM)); CHKERRQ(info);

    278:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_NLMVM,(void*)neP); CHKERRQ(info);
    279:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_NLMVM,TaoDestroy_NLMVM); CHKERRQ(info);
    280:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_NLMVM); CHKERRQ(info);
    281:   info=TaoSetTaoViewRoutine(tao,TaoView_NLMVM); CHKERRQ(info);

    283:   info = TaoSetMaximumIterates(tao,50); CHKERRQ(info);
    284:   info = TaoSetTolerances(tao,1e-16,1e-16,0,0); CHKERRQ(info);
    285:   info = TaoLMVMSetSize(tao,5); CHKERRQ(info);
    286:   info = TaoSetTrustRegionRadius(tao,1.0e16); CHKERRQ(info);

    288:   neP->gamma                        = 0.0;
    289:   neP->gamma_factor                = 0.01;
    290:   neP->max_kspiter_factor       = 5;
    291:   neP->lmvmits                  = 2;

    293:   info = TaoCreateMoreThuenteLineSearch(tao,0,0); CHKERRQ(info);
    294:   /*
    295:   info = PetscObjectComposeFunctionDynamic((PetscObject)tao,
    296:                                     "TaoLineSearchGetDampingParameter_C",
    297:                                     "TaoLineSearchGetDampingParameter_NLMVM",
    298:                                     (void*)TaoLineSearchGetDampingParameter_NLMVM);
    299:   CHKERRQ(info);
    300:   */
    301:   TaoFunctionReturn(0);
    302: }

    308: int TaoLineSearchGetDampingParameter_NLMVM(TAO_SOLVER tao,double *damp)
    309: {
    310:   TAO_NLMVM *neP;
    311:   int info;

    313:   TaoFunctionBegin;
    314:   info = TaoGetSolverContext(tao,"tao_nlmvm",(void**)&neP); CHKERRQ(info);
    315:   if (neP){
    316:     *damp = neP->gamma;
    317:   }
    318:   TaoFunctionReturn(0);
    319: }