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