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