Actual source code: ssils.c

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

  3: /*------------------------------------------------------------*/
  6: static int TaoSolve_SSILS(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 = TaoSetMeritFunction(tao, Tao_SSLS_Function, Tao_SSLS_FunctionGradient,
 39:                              TAO_NULL, TAO_NULL, ssls); CHKERRQ(info);

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

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

 48:   while (1) {
 49:     info=PetscLogInfo((tao, "TaoSolve_SSILS: %d: merit: %5.4e, ndpsi: %5.4e\n",
 50:                        iter, ssls->merit, ndpsi)); CHKERRQ(info);
 51:     // Check the termination criteria
 52:     info = TaoMonitor(tao,iter++,ssls->merit,ndpsi,0.0,t,&reason);
 53:            CHKERRQ(info);
 54:     if (reason!=TAO_CONTINUE_ITERATING) break;

 56:     // Calculate direction.  (Really negative of newton direction.  Therefore,
 57:     // rest of the code uses -d.)
 58:     info = TaoPreLinearSolve(tao, J); CHKERRQ(info);
 59:     info = TaoLinearSolve(tao, J, ff, d, &flag); CHKERRQ(info);
 60:     info = d->Norm2(&normd); CHKERRQ(info);
 61:     info = d->Dot(dpsi, &innerd); CHKERRQ(info);

 63:     // Make sure that we have a descent direction
 64:     if (innerd <= delta*pow(normd, rho)) {
 65:       info = PetscLogInfo((tao, "TaoSolve_SSILS: %d: newton direction not descent\n", iter)); CHKERRQ(info);
 66:       info = d->CopyFrom(dpsi); CHKERRQ(info);
 67:       info = d->Dot(dpsi, &innerd); CHKERRQ(info);
 68:     }

 70:     info = d->Negate(); CHKERRQ(info);
 71:     innerd = -innerd;

 73:     t = 1;
 74:     info = TaoLineSearchApply(tao, x, dpsi, d, w,
 75:                               &psi, &t, &tao->lsflag); CHKERRQ(info);
 76:     info = dpsi->Norm2(&ndpsi);
 77:   }

 79:   TaoFunctionReturn(0);
 80: }

 82: /* ---------------------------------------------------------- */
 86: int TaoCreate_SSILS(TAO_SOLVER tao)
 87: {
 88:   TAO_SSLS *ssls;
 89:   int        info;

 91:   TaoFunctionBegin;

 93:   info = TaoNew(TAO_SSLS,&ssls); CHKERRQ(info);
 94:   info = PetscLogObjectMemory(tao, sizeof(TAO_SSLS)); CHKERRQ(info);

 96:   ssls->delta = 1e-10;
 97:   ssls->rho = 2.1;

 99:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_SSILS,(void*)ssls); CHKERRQ(info);
100:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_SSLS,TaoSetDown_SSLS); CHKERRQ(info);
101:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_SSLS); CHKERRQ(info);
102:   info=TaoSetTaoViewRoutine(tao,TaoView_SSLS); CHKERRQ(info);

104:   info = TaoCreateArmijoLineSearch(tao); CHKERRQ(info);

106:   info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
107:   info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);

109:   info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
110:   info = TaoSetGradientTolerances(tao,1.0e-16,0.0,0.0); CHKERRQ(info);
111:   info = TaoSetFunctionLowerBound(tao,1.0e-8); CHKERRQ(info);

113:   TaoFunctionReturn(0);
114: }