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