Actual source code: blmvm.c

      1: /*$Id: blmvm.c 1.87 05/05/09 11:09:17-05:00 sarich@zorak.(none) $*/

      3: #include "blmvm.h"      /*I "tao_solver.h" I*/
      4: static int TaoSetDown_BLMVM(TAO_SOLVER, void*);

      6: /* ---------------------------------------------------------- */

     10: int TaoSetUp_BLMVM(TAO_SOLVER tao, void*solver){

     12:   TAO_BLMVM *blmvm = (TAO_BLMVM *)solver;
     13:   TaoVec*  X;
     14:   int info, lm;

     16:   TaoFunctionBegin;

     18:   info = TaoGetSolution(tao,&X);CHKERRQ(info);

     20:   info = X->Clone(&blmvm->DX); CHKERRQ(info);
     21:   info = X->Clone(&blmvm->GP); CHKERRQ(info);
     22:   info = X->Clone(&blmvm->G); CHKERRQ(info);
     23:   info = X->Clone(&blmvm->XL); CHKERRQ(info);
     24:   info = X->Clone(&blmvm->XU); CHKERRQ(info);

     26:   info = TaoSetLagrangianGradientVector(tao,blmvm->GP);CHKERRQ(info);
     27:   info = TaoSetStepDirectionVector(tao,blmvm->DX);CHKERRQ(info);
     28:   info = TaoSetVariableBounds(tao,blmvm->XL,blmvm->XU);CHKERRQ(info);
     29:
     30:   info = TaoLMVMGetSize(tao,&lm); CHKERRQ(info);
     31:   blmvm->M = new TaoLMVMMat(X,lm);

     33:   TaoFunctionReturn(0);
     34: }

     36: /* ---------------------------------------------------------- */
     39: static int TaoSetDown_BLMVM(TAO_SOLVER tao, void*solver)
     40: {
     41:   int info;
     42:   TAO_BLMVM *blmvm = (TAO_BLMVM *)solver;

     44:   TaoFunctionBegin;

     46:   info=TaoVecDestroy(blmvm->DX);CHKERRQ(info);blmvm->DX=0;
     47:   info=TaoVecDestroy(blmvm->GP);CHKERRQ(info);blmvm->GP=0;
     48:   info=TaoVecDestroy(blmvm->G);CHKERRQ(info);
     49:   info=TaoVecDestroy(blmvm->XL);CHKERRQ(info);
     50:   info=TaoVecDestroy(blmvm->XU);CHKERRQ(info);
     51:   info=TaoMatDestroy(blmvm->M);CHKERRQ(info);blmvm->M=0;
     52:
     53:   TaoFunctionReturn(0);
     54: }

     56: /*------------------------------------------------------------*/
     59: static int TaoSolve_BLMVM(TAO_SOLVER tao, void *solver){

     61:   TAO_BLMVM *blmvm = (TAO_BLMVM *) solver;
     62:   TaoVec  *X,*G=blmvm->G,*GP=blmvm->GP;         /* solution TaoVec */
     63:   TaoVec  *DX=blmvm->DX;
     64:   TaoVec  *XL=blmvm->XL,*XU=blmvm->XU;
     65:   int       info,lsflag,iter=0;
     66:   double    f,gnorm,step=0,gdx;
     67:   TaoTerminateReason reason;
     68:   TaoTruth success;

     70:   TaoFunctionBegin;

     72:   info = TaoGetSolution(tao,&X); CHKERRQ(info);

     74:   info = TaoEvaluateVariableBounds(tao,XL,XU); CHKERRQ(info);
     75:   info = X->Median(XL,X,XU); CHKERRQ(info);

     77:   info = TaoComputeMeritFunctionGradient(tao,X,&f,G); CHKERRQ(info);
     78:   blmvm->pgits=0;

     80:   /* Enter loop */
     81:   while (1){

     83:     info = GP->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);
     84:     info = GP->Norm2(&gnorm); CHKERRQ(info);

     86:     /* Test for convergence */
     87:     info = TaoMonitor(tao,iter++,f,gnorm,0.0,step,&reason); CHKERRQ(info);
     88:     if (reason!=TAO_CONTINUE_ITERATING) break;

     90:     info = blmvm->M->update(X,GP); CHKERRQ(info);
     91:     info = blmvm->M->Solve(G,DX,&success); CHKERRQ(info);
     92:     info = GP->BoundGradientProjection(DX,XL,X,XU); CHKERRQ(info);

     94:     info = GP->Dot(G,&gdx); CHKERRQ(info);
     95:     if (gdx<=0){
     96:       info = DX->ScaleCopyFrom(-1.0,G); CHKERRQ(info);
     97:       gdx=-gnorm*gnorm;
     98:       blmvm->pgits++;
     99:     } else {
    100:       info = DX->Negate(); CHKERRQ(info);
    101:       gdx*= -1;
    102:     }

    104:     /* Line Search */
    105:     if (iter==1){
    106:       info = TaoGetInitialTrustRegionRadius(tao,&step); CHKERRQ(info);
    107:       if (step <= 0 ) step=1.0;
    108:     } else step=1.0;

    110:     info = TaoLineSearchApply(tao,X,G,DX,GP,&f,&step,&lsflag);
    111:     CHKERRQ(info)

    113:   }

    115:   TaoFunctionReturn(0);
    116: }

    119: /*------------------------------------------------------------*/
    122: static int TaoSetOptions_BLMVM(TAO_SOLVER tao, void*solver)
    123: {
    124:   int        info;

    126:   TaoFunctionBegin;
    127:   info = TaoOptionsHead("Variable metric method for bound constrained optimization");CHKERRQ(info);
    128:   info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);
    129:   info = TaoOptionsTail();CHKERRQ(info);
    130:   TaoFunctionReturn(0);
    131: }

    134: /*------------------------------------------------------------*/
    137: static int TaoView_BLMVM(TAO_SOLVER tao,void*solver)
    138: {
    139:   TAO_BLMVM *blmvm = (TAO_BLMVM *) solver;
    140:   int        info;

    142:   TaoFunctionBegin;
    143:
    144:   info=TaoPrintInt(tao,"  Projected gradient iterates: %d\n",blmvm->pgits); CHKERRQ(info);
    145:   info=TaoPrintInt(tao,"  Rejected matrix updates: %d\n",blmvm->M->countrejects());
    146:   CHKERRQ(info);
    147:
    148:   info = TaoLineSearchView(tao);CHKERRQ(info);

    150:   TaoFunctionReturn(0);
    151: }

    155: static int TaoGetDualVariables_BLMVM(TAO_SOLVER tao, TaoVec *DXL, TaoVec* DXU, void* solver)
    156: {
    157:   TAO_BLMVM *blmvm = (TAO_BLMVM *) solver;
    158:   TaoVec  *G=blmvm->G,*GP=blmvm->GP;
    159:   int       info;

    161:   TaoFunctionBegin;
    162:
    163:   info = DXL->Waxpby(-1.0,G,1.0,GP); CHKERRQ(info);
    164:   info = DXU->SetToZero(); CHKERRQ(info);
    165:   info = DXL->PointwiseMaximum(DXL,DXU); CHKERRQ(info);

    167:   info = DXU->Waxpby(-1.0,GP,1.0,G); CHKERRQ(info);
    168:   info = DXU->Axpy(1.0,DXL); CHKERRQ(info);

    170:   TaoFunctionReturn(0);
    171: }

    174: /* ---------------------------------------------------------- */
    178: int TaoCreate_BLMVM(TAO_SOLVER tao)
    179: {
    180:   TAO_BLMVM *blmvm;
    181:   int       info;

    183:   TaoFunctionBegin;

    185:   info = TaoNew(TAO_BLMVM,&blmvm); CHKERRQ(info);
    186:   info = PetscLogObjectMemory(tao,sizeof(TAO_BLMVM)); CHKERRQ(info);

    188:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_BLMVM,(void*)blmvm); CHKERRQ(info);
    189:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_BLMVM,TaoSetDown_BLMVM); CHKERRQ(info);
    190:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_BLMVM); CHKERRQ(info);
    191:   info=TaoSetTaoViewRoutine(tao,TaoView_BLMVM); CHKERRQ(info);
    192:   info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_BLMVM); CHKERRQ(info);

    194:   info = TaoSetMaximumIterates(tao,10000); CHKERRQ(info);
    195:   info = TaoSetMaximumFunctionEvaluations(tao,10000); CHKERRQ(info);
    196:   info = TaoSetTolerances(tao,1e-4,1e-4,0,0); CHKERRQ(info);
    197:   info = TaoLMVMSetSize(tao,5); CHKERRQ(info);
    198:   info = TaoSetTrustRegionRadius(tao,0.01);CHKERRQ(info);

    200:   //  info = TaoCreateProjectedLineSearch(tao); CHKERRQ(info);
    201:   info = TaoCreateMoreThuenteBoundLineSearch(tao,0,0.99); CHKERRQ(info);

    203:   blmvm->DX=0;blmvm->GP=0;blmvm->M=0;
    204:   TaoFunctionReturn(0);
    205: }