Actual source code: ssfls.c

      1: #include "src/complementarity/impls/ssls/ssls.h"

      3: /*------------------------------------------------------------*/
      6: static int TaoSolve_SSFLS(TAO_SOLVER tao, void *solver)
      7: {
      8:   TAO_SSLS *ssls = (TAO_SSLS *)solver;
      9:   //  TaoLinearSolver *lsolver;
     10:   TaoVec *x, *l, *u, *ff, *dpsi, *d, *w;
     11:   TaoMat *J;
     12:   double psi, ndpsi, normd, innerd, t=0;
     13:   double delta, rho;
     14:   int iter=0, info;
     15:   TaoTerminateReason reason;
     16:   TaoTruth flag;

     18:   TaoFunctionBegin;

     20:   // Assume that Setup has been called!
     21:   // Set the structure for the Jacobian and create a linear solver.
     22:
     23:   delta = ssls->delta;
     24:   rho = ssls->rho;

     26:   info = TaoGetSolution(tao, &x); CHKERRQ(info);
     27:   l=ssls->xl;
     28:   u=ssls->xu;
     29:   info = TaoEvaluateVariableBounds(tao,l,u); CHKERRQ(info);
     30:   info = x->Median(l,x,u); CHKERRQ(info);
     31:   info = TaoGetJacobian(tao, &J); CHKERRQ(info);

     33:   ff = ssls->ff;
     34:   dpsi = ssls->dpsi;
     35:   d = ssls->d;
     36:   w = ssls->w;

     38:   info = x->PointwiseMaximum(x, l); CHKERRQ(info);
     39:   info = x->PointwiseMinimum(x, u); CHKERRQ(info);
     40:   info = TaoSetMeritFunction(tao, Tao_SSLS_Function, Tao_SSLS_FunctionGradient,
     41:                              TAO_NULL, TAO_NULL, ssls); CHKERRQ(info);

     43:   info = J->SetStructure(TAOMAT_UNSYMMETRIC); CHKERRQ(info);

     45:   // Calculate the function value and fischer function value at the
     46:   // current iterate
     47:   info = TaoComputeMeritFunctionGradient(tao, x, &psi, dpsi); CHKERRQ(info);
     48:   info = dpsi->Norm2(&ndpsi);

     50:   while (1) {
     51:     info=PetscLogInfo((tao, "TaoSolve_SSFLS: %d: merit: %5.4e, ndpsi: %5.4e\n",
     52:                        iter, ssls->merit, ndpsi));CHKERRQ(info);

     54:     // Check the termination criteria
     55:     info = TaoMonitor(tao,iter++,ssls->merit,ndpsi,0.0,t,&reason);
     56:            CHKERRQ(info);
     57:     if (reason!=TAO_CONTINUE_ITERATING) break;

     59:     // Calculate direction.  (Really negative of newton direction.  Therefore,
     60:     // rest of the code uses -d.)
     61:     info = TaoPreLinearSolve(tao, J); CHKERRQ(info);
     62:     info = TaoLinearSolve(tao, J, ff, d, &flag); CHKERRQ(info);
     63:
     64:     info = w->CopyFrom(d); CHKERRQ(info);
     65:     info = w->Negate(); CHKERRQ(info);
     66:     info = w->BoundGradientProjection(w,l, x, u);

     68:     info = w->Norm2(&normd); CHKERRQ(info);
     69:     info = w->Dot(dpsi, &innerd); CHKERRQ(info);

     71:     // Make sure that we have a descent direction
     72:     if (innerd >= -delta*pow(normd, rho)) {
     73:       info = PetscLogInfo((tao, "TaoSolve_SSFLS: %d: newton direction not descent\n", iter)); CHKERRQ(info);
     74:       info = d->CopyFrom(dpsi); CHKERRQ(info);
     75:       info = w->Dot(dpsi, &innerd); CHKERRQ(info);
     76:     }
     77:     info = d->Negate(); CHKERRQ(info);
     78:     innerd = -innerd;

     80:     t = 1;
     81:     info = TaoLineSearchApply(tao, x, dpsi, d, w,
     82:                               &psi, &t, &tao->lsflag); CHKERRQ(info);
     83:     info = dpsi->Norm2(&ndpsi);
     84:   }
     85:   TaoFunctionReturn(0);
     86: }

     88: /* ---------------------------------------------------------- */
     92: int TaoCreate_SSFLS(TAO_SOLVER tao)
     93: {
     94:   TAO_SSLS *ssls;
     95:   int        info;

     97:   TaoFunctionBegin;

     99:   info = TaoNew(TAO_SSLS,&ssls); CHKERRQ(info);
    100:   info = PetscLogObjectMemory(tao, sizeof(TAO_SSLS)); CHKERRQ(info);

    102:   ssls->delta = 1e-10;
    103:   ssls->rho = 2.1;

    105:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_SSFLS,(void*)ssls); CHKERRQ(info);
    106:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_SSLS,TaoSetDown_SSLS); CHKERRQ(info);
    107:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_SSLS); CHKERRQ(info);
    108:   info=TaoSetTaoViewRoutine(tao,TaoView_SSLS); CHKERRQ(info);

    110:   info = TaoCreateProjectedArmijoLineSearch(tao); CHKERRQ(info);

    112:   info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
    113:   info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);

    115:   info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
    116:   info = TaoSetGradientTolerances(tao,1.0e-16,0.0,0.0); CHKERRQ(info);
    117:   info = TaoSetFunctionLowerBound(tao,1.0e-8); CHKERRQ(info);

    119:   TaoFunctionReturn(0);
    120: }