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