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