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