1: #include <../src/tao/complementarity/impls/ssls/ssls.h>

  5: PetscErrorCode TaoSetUp_SSILS(Tao tao)
  6: {
  7:   TAO_SSLS       *ssls = (TAO_SSLS *)tao->data;

 11:   VecDuplicate(tao->solution,&tao->gradient);
 12:   VecDuplicate(tao->solution,&tao->stepdirection);
 13:   VecDuplicate(tao->solution,&ssls->ff);
 14:   VecDuplicate(tao->solution,&ssls->dpsi);
 15:   VecDuplicate(tao->solution,&ssls->da);
 16:   VecDuplicate(tao->solution,&ssls->db);
 17:   VecDuplicate(tao->solution,&ssls->t1);
 18:   VecDuplicate(tao->solution,&ssls->t2);
 19:   return(0);
 20: }

 24: PetscErrorCode TaoDestroy_SSILS(Tao tao)
 25: {
 26:   TAO_SSLS       *ssls = (TAO_SSLS *)tao->data;

 30:   VecDestroy(&ssls->ff);
 31:   VecDestroy(&ssls->dpsi);
 32:   VecDestroy(&ssls->da);
 33:   VecDestroy(&ssls->db);
 34:   VecDestroy(&ssls->t1);
 35:   VecDestroy(&ssls->t2);
 36:   PetscFree(tao->data);
 37:   return(0);
 38: }

 42: static PetscErrorCode TaoSolve_SSILS(Tao tao)
 43: {
 44:   TAO_SSLS                     *ssls = (TAO_SSLS *)tao->data;
 45:   PetscReal                    psi, ndpsi, normd, innerd, t=0;
 46:   PetscReal                    delta, rho;
 47:   PetscInt                     iter=0,kspits;
 48:   TaoConvergedReason           reason;
 49:   TaoLineSearchConvergedReason ls_reason;
 50:   PetscErrorCode               ierr;

 53:   /* Assume that Setup has been called!
 54:      Set the structure for the Jacobian and create a linear solver. */
 55:   delta = ssls->delta;
 56:   rho = ssls->rho;

 58:   TaoComputeVariableBounds(tao);
 59:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
 60:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);
 61:   TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);

 63:   /* Calculate the function value and fischer function value at the
 64:      current iterate */
 65:   TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi);
 66:   VecNorm(ssls->dpsi,NORM_2,&ndpsi);

 68:   while (1) {
 69:     ierr=PetscInfo3(tao, "iter: %D, merit: %g, ndpsi: %g\n",iter, (double)ssls->merit, (double)ndpsi);
 70:     /* Check the termination criteria */
 71:     TaoMonitor(tao,iter++,ssls->merit,ndpsi,0.0,t,&reason);
 72:     if (reason!=TAO_CONTINUE_ITERATING) break;

 74:     /* Calculate direction.  (Really negative of newton direction.  Therefore,
 75:        rest of the code uses -d.) */
 76:     KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre);
 77:     KSPSolve(tao->ksp,ssls->ff,tao->stepdirection);
 78:     KSPGetIterationNumber(tao->ksp,&kspits);
 79:     tao->ksp_its+=kspits;
 80:     VecNorm(tao->stepdirection,NORM_2,&normd);
 81:     VecDot(tao->stepdirection,ssls->dpsi,&innerd);

 83:     /* Make sure that we have a descent direction */
 84:     if (innerd <= delta*pow(normd, rho)) {
 85:       PetscInfo(tao, "newton direction not descent\n");
 86:       VecCopy(ssls->dpsi,tao->stepdirection);
 87:       VecDot(tao->stepdirection,ssls->dpsi,&innerd);
 88:     }

 90:     VecScale(tao->stepdirection, -1.0);
 91:     innerd = -innerd;

 93:     TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
 94:     TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason);
 95:     VecNorm(ssls->dpsi,NORM_2,&ndpsi);
 96:   }
 97:   return(0);
 98: }

100: /* ---------------------------------------------------------- */
101: /*MC
102:    TAOSSILS - semi-smooth infeasible linesearch algorithm for solving
103:        complementarity constraints

105:    Options Database Keys:
106: + -tao_ssls_delta - descent test fraction
107: - -tao_ssls_rho - descent test power

109:    Level: beginner
110: M*/
111: EXTERN_C_BEGIN
114: PetscErrorCode TaoCreate_SSILS(Tao tao)
115: {
116:   TAO_SSLS       *ssls;
118:   const char     *armijo_type = TAOLINESEARCHARMIJO;

121:   PetscNewLog(tao,&ssls);
122:   tao->data = (void*)ssls;
123:   tao->ops->solve=TaoSolve_SSILS;
124:   tao->ops->setup=TaoSetUp_SSILS;
125:   tao->ops->view=TaoView_SSLS;
126:   tao->ops->setfromoptions=TaoSetFromOptions_SSLS;
127:   tao->ops->destroy = TaoDestroy_SSILS;

129:   ssls->delta = 1e-10;
130:   ssls->rho = 2.1;

132:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
133:   TaoLineSearchSetType(tao->linesearch,armijo_type);
134:   TaoLineSearchSetFromOptions(tao->linesearch);
135:   /* Note: linesearch objective and objectivegradient routines are set in solve routine */
136:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);

138:   tao->max_it = 2000;
139:   tao->max_funcs = 4000;
140:   tao->fatol = 0;
141:   tao->frtol = 0;
142:   tao->gttol=0;
143:   tao->grtol=0;
144: #if defined(PETSC_USE_REAL_SINGLE)
145:   tao->gatol = 1.0e-6;
146:   tao->fmin = 1.0e-4;
147: #else
148:   tao->gatol = 1.0e-16;
149:   tao->fmin = 1.0e-8;
150: #endif
151:   return(0);
152: }
153: EXTERN_C_END