Actual source code: lmvm.c

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

      3: #include "lmvm.h"             /*I "tao_solver.h" I*/
      4: #include "src/matrix/lmvmmat.h"             /*I "tao_solver.h" I*/

      6: static int TaoDestroy_LMVM(TAO_SOLVER, void*);

      8: /* ---  TESTING ------------------------------------------------------- */
     11: int TaoInitializeLMVMmatrix(TAO_SOLVER tao, TaoVec *HV){
     12:   int info;
     13:   TAO_LMVM *lmvm;
     14:
     15:   TaoFunctionBegin;
     16:   info = TaoGetSolverContext(tao,"tao_lmvm",(void**)&lmvm); CHKERRQ(info);
     17:   if (lmvm && lmvm->M){
     18:     info = lmvm->M->InitialApproximation( HV );CHKERRQ(info);
     19:   }
     20:   TaoFunctionReturn(0);
     21: }

     25: int TaoSetUp_LMVM(TAO_SOLVER tao, void* solver)
     26: {
     27:   TAO_LMVM *lmvm = (TAO_LMVM *)solver;
     28:   TaoVec  *xx;
     29:   int    info, lm;

     31:   TaoFunctionBegin;

     33:   info = TaoGetSolution(tao,&xx);CHKERRQ(info);

     35:   info = TaoLMVMGetSize(tao,&lm); CHKERRQ(info);
     36:   lmvm->M = new TaoLMVMMat(xx,lm);

     38:   info = xx->Clone(&lmvm->G); CHKERRQ(info);
     39:   info = xx->Clone(&lmvm->Work); CHKERRQ(info);
     40:   info = xx->Clone(&lmvm->DX); CHKERRQ(info);

     42:   info = TaoSetLagrangianGradientVector(tao,lmvm->G);CHKERRQ(info);
     43:   info = TaoSetStepDirectionVector(tao,lmvm->DX);CHKERRQ(info);

     45:   TaoFunctionReturn(0);
     46: }

     48: /* ---------------------------------------------------------- */
     51: static int TaoDestroy_LMVM(TAO_SOLVER tao, void*solver)
     52: {
     53:   int info;
     54:   TAO_LMVM *lmvm = (TAO_LMVM *)solver;

     56:   TaoFunctionBegin;

     58:   info=TaoVecDestroy(lmvm->G);CHKERRQ(info);
     59:   info=TaoVecDestroy(lmvm->DX);CHKERRQ(info);
     60:   info=TaoVecDestroy(lmvm->Work);CHKERRQ(info);
     61:   info=TaoMatDestroy(lmvm->M);CHKERRQ(info);

     63:   info = TaoSetLagrangianGradientVector(tao,0);CHKERRQ(info);
     64:   info = TaoSetStepDirectionVector(tao,0);CHKERRQ(info);
     65:   lmvm->DX=0;lmvm->Work=0;lmvm->M=0;

     67:   TaoFunctionReturn(0);
     68: }

     70: /*------------------------------------------------------------*/
     73: static int TaoSolve_LMVM(TAO_SOLVER tao, void *solver)
     74: {
     75:   TAO_LMVM *lmvm = (TAO_LMVM *)solver;
     76:   TaoVec       *xx, *gg=lmvm->G;         /* gradient vector */
     77:   TaoVec       *dxdx=lmvm->DX, *ww=lmvm->Work;
     78:   int          info,lsflag=0,iter=0;
     79:   double       f, gdx, gnorm, dxnorm, r0, step=0;
     80:   TaoTerminateReason reason;
     81:   TaoTruth     success;

     83:   TaoFunctionBegin;

     85:   info = TaoGetSolution(tao,&xx);CHKERRQ(info);

     87:   info = TaoComputeMeritFunctionGradient(tao,xx,&f,gg);CHKERRQ(info);
     88:   info = gg->Norm2(&gnorm);  CHKERRQ(info);

     90:   /* Enter loop */
     91:   while (1){

     93:     /* Test for convergence */
     94:     info=TaoMonitor(tao,iter++,f,gnorm,0,step,&reason);CHKERRQ(info);
     95:     if (reason!=TAO_CONTINUE_ITERATING) break;

     97:     /* Compute step direction */
     98:     info = lmvm->M->update(xx,gg); CHKERRQ(info);
     99:     info = lmvm->M->Solve(gg,dxdx,&success); CHKERRQ(info);
    100:     info = dxdx->Negate(); CHKERRQ(info);
    101:     info = dxdx->Dot(gg,&gdx); CHKERRQ(info);
    102:     if (gdx>=0){
    103:       info = dxdx->ScaleCopyFrom(-1.0,gg); CHKERRQ(info);
    104:       gdx = -gnorm*gnorm;
    105:     }

    107:     /* Line Search */
    108:     info = dxdx->Norm2(&dxnorm);  CHKERRQ(info);
    109:     info = TaoGetInitialTrustRegionRadius(tao,&r0); CHKERRQ(info);
    110:     step = TaoMin(1.0,r0/dxnorm);
    111:     if (step < 1.0){
    112:       info = PetscLogInfo((tao,"Restrict step size to %g so the change in solution <= %g\n",step,r0)); CHKERRQ(info);
    113:     }
    114:     info = TaoLineSearchApply(tao,xx,gg,dxdx,ww,&f,&step,&lsflag);
    115:     CHKERRQ(info);
    116:     info = gg->Norm2(&gnorm);CHKERRQ(info);

    118:   }

    120:   TaoFunctionReturn(0);

    122: }

    125: /*------------------------------------------------------------*/
    128: static int TaoSetOptions_LMVM(TAO_SOLVER tao, void *solver)
    129: {
    130:   int        info;

    132:   TaoFunctionBegin;
    133:   info = TaoOptionsHead("Limited Memory Variable Metric method for unconstrained optimization");CHKERRQ(info);
    134:   info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);
    135:   info = TaoOptionsTail();CHKERRQ(info);

    137:   TaoFunctionReturn(0);
    138: }

    140: /* ---------------------------------------------------------- */
    143: int TaoLMVMGetX0(TAO_SOLVER tao,TaoVec *x0)
    144: {
    145:   TAO_LMVM *lmvm;
    146:   int info;
    147:   TaoFunctionBegin;
    148:   TaoValidHeaderSpecific(tao,TAO_COOKIE,1);
    149:   info=TaoGetSolverContext(tao,"tao_lmvm",(void**)&lmvm);CHKERRQ(info);
    150:   if (lmvm && lmvm->M){
    151:     info=lmvm->M->GetX0(x0);CHKERRQ(info);
    152:   }
    153:   TaoFunctionReturn(0);
    154: }

    156: /*------------------------------------------------------------*/
    159: static int TaoView_LMVM(TAO_SOLVER tao,void *solver)
    160: {
    161:   int        info;

    163:   TaoFunctionBegin;
    164:   info = TaoLineSearchView(tao);CHKERRQ(info);
    165:   TaoFunctionReturn(0);
    166: }

    168: /* ---------------------------------------------------------- */

    173: int TaoCreate_LMVM(TAO_SOLVER tao)
    174: {
    175:   TAO_LMVM *lmvm;
    176:   int info;

    178:   TaoFunctionBegin;

    180:   info = TaoNew(TAO_LMVM,&lmvm); CHKERRQ(info);
    181:   info = PetscLogObjectMemory(tao,sizeof(TAO_LMVM)); CHKERRQ(info);

    183:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_LMVM,(void*)lmvm); CHKERRQ(info);
    184:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_LMVM,TaoDestroy_LMVM); CHKERRQ(info);
    185:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_LMVM); CHKERRQ(info);
    186:   info=TaoSetTaoViewRoutine(tao,TaoView_LMVM); CHKERRQ(info);

    188:   info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
    189:   info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);
    190:   info = TaoSetTolerances(tao,1e-4,1e-4,0,0); CHKERRQ(info);
    191:   info = TaoLMVMSetSize(tao,5); CHKERRQ(info);
    192:   info = TaoSetTrustRegionRadius(tao,1.0e16); CHKERRQ(info);
    193:   info = TaoCreateMoreThuenteLineSearch(tao,0,0.99); CHKERRQ(info);
    194:   lmvm->DX=0;lmvm->Work=0;lmvm->M=0;

    196:   TaoFunctionReturn(0);
    197: }