Actual source code: nls.c
1: /*$Id: nls.c 1.111 05/05/10 15:21:39-05:00 sarich@zorak.(none) $*/
3: #include "nls.h" /*I "tao_solver.h" I*/
5: static int TaoDestroy_NLS(TAO_SOLVER tao, void *solver);
6: static int TaoView_NLS(TAO_SOLVER tao,void* solver);
7: static int TaoSolve_NLS(TAO_SOLVER tao,void *solver);
8: static int TaoSetOptions_NLS(TAO_SOLVER tao, void *solver);
9: static int TaoSetUp_NLS(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: */
23: static int TaoSolve_NLS(TAO_SOLVER tao,void *solver)
24: {
25: TAO_NLS *neP = (TAO_NLS *) solver;
26: int iter=0, info, line=0;
27: double f, gdx, gnorm, two = 2.0,step=0.0;
28: TaoVec *gg, *xx;
29: TaoVec *rrhs=neP->RHS, *ss=neP->S, *ww=neP->W;
30: TaoMat *HH;
31: TaoTerminateReason reason;
32: TaoTruth success;
33: double gamma,gamma_factor;
35: TaoFunctionBegin;
37: info=TaoGetSolution(tao,&xx);CHKERRQ(info);
38: info=TaoGetGradient(tao,&gg);CHKERRQ(info);
39: info=TaoGetHessian(tao,&HH);CHKERRQ(info);
41: info = TaoComputeMeritFunctionGradient(tao,xx,&f,gg);CHKERRQ(info);
42: info = gg->Norm2(&gnorm);CHKERRQ(info); /* gnorm = || gg || */
44: gamma = neP->gamma;
45: gamma_factor= neP->gamma_factor;
46:
47: while (1) {
48:
49: info = TaoMonitor(tao,iter++,f,gnorm,0.0,step,&reason);CHKERRQ(info);
50: if (reason!=TAO_CONTINUE_ITERATING) break;
52: gamma_factor /= two;
53: gamma = gamma_factor*(gnorm);
55: info = rrhs->Axpby(-1.0,gg,0.0); CHKERRQ(info);
57: info = TaoComputeHessian(tao,xx,HH);CHKERRQ(info);
59: success = TAO_FALSE;
60: while (success==TAO_FALSE) {
61: info = TaoPreLinearSolve(tao,HH);CHKERRQ(info);
62: info = TaoLinearSolve(tao,HH,rrhs,ss,&success);CHKERRQ(info);
64: info = ss->Dot(gg,&gdx); CHKERRQ(info);
66: if (gdx>=0 /*|| success==TAO_FALSE*/) { /* Modify diagonal of Hessian */
67: gamma_factor *= two;
68: gamma = gamma_factor*(gnorm);
69: #if !defined(PETSC_USE_COMPLEX)
70: info = PetscLogInfo((tao,"TaoSolve_NLS: modify diagonal (assume same nonzero structure), gamma_factor=%g, gamma=%g\n", gamma_factor,gamma));CHKERRQ(info);
71: #else
72: info = PetscLogInfo((tao,"TaoSolve_NLS: modify diagonal (assume same nonzero structure), gamma_factor=%g, gamma=%g\n",
73: gamma_factor,PetscReal(gamma))); CHKERRQ(info);
74: #endif
75: info = HH->ShiftDiagonal(gamma);CHKERRQ(info);
76: success = TAO_FALSE;
78: /* We currently assume that all diagonal elements were allocated in
79: original matrix, so that nonzero pattern is same ... should fix this */
80: } else {
81: // gamma_factor /= two;
82: success = TAO_TRUE;
83: }
84: }
86: /* Line search */
87: step=1.0;
88: info = TaoLineSearchApply(tao,xx,gg,ss,ww,&f,&step,&line);
89: CHKERRQ(info);
90: info = gg->Norm2(&gnorm);CHKERRQ(info);
91: // neP->gamma=gamma;
92: //. neP->gamma_factor=gamma_factor;
94: }
97: TaoFunctionReturn(0);
98: }
99: /* ---------------------------------------------------------- */
102: static int TaoSetUp_NLS(TAO_SOLVER tao, void *solver)
103: {
104: int info,kspmaxit,n;
105: TAO_NLS *ctx = (TAO_NLS *)solver;
106: TaoVec *xx;
107: TaoMat *HH;
108: TaoLinearSolver *tls;
110: TaoFunctionBegin;
111: info = TaoGetSolution(tao,&xx);CHKERRQ(info);
113: info=xx->Clone(&ctx->RHS);CHKERRQ(info);
114: info=xx->Clone(&ctx->S);CHKERRQ(info);
115: info=xx->Clone(&ctx->W);CHKERRQ(info);
116: info=xx->Clone(&ctx->G);CHKERRQ(info);
118: info = TaoSetLagrangianGradientVector(tao,ctx->G);CHKERRQ(info);
119: info = TaoSetStepDirectionVector(tao,ctx->S);CHKERRQ(info);
121: info = TaoGetHessian(tao,&HH); CHKERRQ(info);
122: info = TaoCreateLinearSolver(tao,HH,200,0); CHKERRQ(info);
123: info = TaoGetLinearSolver(tao,&tls); CHKERRQ(info);
124: info = xx->GetDimension(&n);CHKERRQ(info);
125: kspmaxit = ctx->max_kspiter_factor * ((int) sqrt((double)n));
126: info = tls->SetTolerances(1e-5,1e-40,1e30,kspmaxit);CHKERRQ(info);
128: info = TaoCheckFGH(tao);CHKERRQ(info);
130: TaoFunctionReturn(0);
131: }
132: /*------------------------------------------------------------*/
135: static int TaoDestroy_NLS(TAO_SOLVER tao, void *solver)
136: {
137: TAO_NLS *ctx = (TAO_NLS *)solver;
138: int info;
140: TaoFunctionBegin;
142: info=TaoVecDestroy(ctx->RHS);CHKERRQ(info);
143: info=TaoVecDestroy(ctx->S);CHKERRQ(info);
144: info=TaoVecDestroy(ctx->W);CHKERRQ(info);
145: info=TaoVecDestroy(ctx->G);CHKERRQ(info);
146: info = TaoDestroyLinearSolver(tao);CHKERRQ(info);
148: info = TaoSetLagrangianGradientVector(tao,0);CHKERRQ(info);
149: info = TaoSetStepDirectionVector(tao,0);CHKERRQ(info);
151: TaoFunctionReturn(0);
152: }
153: /*------------------------------------------------------------*/
156: static int TaoSetOptions_NLS(TAO_SOLVER tao, void *solver)
157: {
158: TAO_NLS *ctx = (TAO_NLS *)solver;
159: int info;
160: TaoTruth flg;
162: TaoFunctionBegin;
163: info = TaoOptionsHead("Newton line search method for unconstrained optimization");CHKERRQ(info);
165: info = TaoOptionDouble("-tao_nls_gamma_factor","damping parameter","",ctx->gamma_factor,&ctx->gamma_factor,&flg);CHKERRQ(info);
167: info = TaoOptionInt("-tao_nls_maxkspf","max ksp iterates","",ctx->max_kspiter_factor,&ctx->max_kspiter_factor,&flg);CHKERRQ(info);
169: info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);
170: /* info = TaoSetLinearSolverOptions(tao);CHKERRQ(info); */
172: info = TaoOptionsTail();CHKERRQ(info);
174: TaoFunctionReturn(0);
175: }
178: /*------------------------------------------------------------*/
181: static int TaoView_NLS(TAO_SOLVER tao,void* solver)
182: {
183: TAO_NLS *ctx = (TAO_NLS *)solver;
184: int info;
186: TaoFunctionBegin;
188: info = TaoPrintDouble(tao," gamma_f=%g,",ctx->gamma_factor);CHKERRQ(info);
189: info = TaoPrintInt(tao," maxkspf=%d\n",ctx->max_kspiter_factor);CHKERRQ(info);
191: info = TaoPrintStatement(tao," Linear Solver for symmetric matrix: \n");CHKERRQ(info);
192: info = TaoLineSearchView(tao);CHKERRQ(info);
194: TaoFunctionReturn(0);
195: }
197: /* ---------------------------------------------------------- */
201: int TaoCreate_NLS(TAO_SOLVER tao)
202: {
203: TAO_NLS *neP;
204: int info;
206: TaoFunctionBegin;
208: info = TaoNew(TAO_NLS,&neP); CHKERRQ(info);
209: info = PetscLogObjectMemory(tao,sizeof(TAO_NLS)); CHKERRQ(info);
211: info=TaoSetTaoSolveRoutine(tao,TaoSolve_NLS,(void*)neP); CHKERRQ(info);
212: info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_NLS,TaoDestroy_NLS); CHKERRQ(info);
213: info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_NLS); CHKERRQ(info);
214: info=TaoSetTaoViewRoutine(tao,TaoView_NLS); CHKERRQ(info);
216: info = TaoSetMaximumIterates(tao,50); CHKERRQ(info);
217: info = TaoSetTolerances(tao,1e-16,1e-16,0,0); CHKERRQ(info);
219: neP->gamma = 0.0;
220: neP->gamma_factor = 0.01;
221: neP->max_kspiter_factor = 5;
223: info = TaoCreateMoreThuenteLineSearch(tao,0,0); CHKERRQ(info);
224: /*
225: info = PetscObjectComposeFunctionDynamic((PetscObject)tao,
226: "TaoLineSearchGetDampingParameter_C",
227: "TaoLineSearchGetDampingParameter_NLS",
228: (void*)TaoLineSearchGetDampingParameter_NLS);
229: CHKERRQ(info);
230: */
231: TaoFunctionReturn(0);
232: }
238: int TaoLineSearchGetDampingParameter_NLS(TAO_SOLVER tao,double *damp)
239: {
240: TAO_NLS *neP;
241: int info;
243: TaoFunctionBegin;
244: info = TaoGetSolverContext(tao,"tao_nls",(void**)&neP); CHKERRQ(info);
245: if (neP){
246: *damp = neP->gamma;
247: }
248: TaoFunctionReturn(0);
249: }