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

  5: PetscErrorCode TaoSetUp_SSFLS(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->w);
 14:   VecDuplicate(tao->solution,&ssls->ff);
 15:   VecDuplicate(tao->solution,&ssls->dpsi);
 16:   VecDuplicate(tao->solution,&ssls->da);
 17:   VecDuplicate(tao->solution,&ssls->db);
 18:   VecDuplicate(tao->solution,&ssls->t1);
 19:   VecDuplicate(tao->solution,&ssls->t2);
 20:   if (!tao->XL) {
 21:     VecDuplicate(tao->solution,&tao->XL);
 22:   }
 23:   if (!tao->XU) {
 24:     VecDuplicate(tao->solution,&tao->XU);
 25:   }
 26:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
 27:   return(0);
 28: }

 32: static PetscErrorCode TaoSolve_SSFLS(Tao tao)
 33: {
 34:   TAO_SSLS                     *ssls = (TAO_SSLS *)tao->data;
 35:   PetscReal                    psi, ndpsi, normd, innerd, t=0;
 36:   PetscReal                    delta, rho;
 37:   PetscInt                     iter=0,kspits;
 38:   TaoConvergedReason           reason;
 39:   TaoLineSearchConvergedReason ls_reason;
 40:   PetscErrorCode               ierr;

 43:   /* Assume that Setup has been called!
 44:      Set the structure for the Jacobian and create a linear solver. */
 45:   delta = ssls->delta;
 46:   rho = ssls->rho;

 48:   TaoComputeVariableBounds(tao);
 49:   /* Project solution inside bounds */
 50:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
 51:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);
 52:   TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);

 54:   /* Calculate the function value and fischer function value at the
 55:      current iterate */
 56:   TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi);
 57:   VecNorm(ssls->dpsi,NORM_2,&ndpsi);

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

 65:     /* Calculate direction.  (Really negative of newton direction.  Therefore,
 66:        rest of the code uses -d.) */
 67:     KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre);
 68:     KSPSolve(tao->ksp,ssls->ff,tao->stepdirection);
 69:     KSPGetIterationNumber(tao->ksp,&kspits);
 70:     tao->ksp_its+=kspits;

 72:     VecCopy(tao->stepdirection,ssls->w);
 73:     VecScale(ssls->w,-1.0);
 74:     VecBoundGradientProjection(ssls->w,tao->solution,tao->XL,tao->XU,ssls->w);

 76:     VecNorm(ssls->w,NORM_2,&normd);
 77:     VecDot(ssls->w,ssls->dpsi,&innerd);

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

 86:     VecScale(tao->stepdirection, -1.0);
 87:     innerd = -innerd;

 89:     TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
 90:     TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason);
 91:     VecNorm(ssls->dpsi,NORM_2,&ndpsi);
 92:   }
 93:   return(0);
 94: }

 98: PetscErrorCode TaoDestroy_SSFLS(Tao tao)
 99: {
100:   TAO_SSLS       *ssls = (TAO_SSLS *)tao->data;

104:   VecDestroy(&ssls->ff);
105:   VecDestroy(&ssls->w);
106:   VecDestroy(&ssls->dpsi);
107:   VecDestroy(&ssls->da);
108:   VecDestroy(&ssls->db);
109:   VecDestroy(&ssls->t1);
110:   VecDestroy(&ssls->t2);
111:   PetscFree(tao->data);
112:   return(0);
113: }

115: /* ---------------------------------------------------------- */
116: /*MC
117:    TAOSSFLS - Semi-smooth feasible linesearch algorithm for solving
118:        complementarity constraints

120:    Options Database Keys:
121: + -tao_ssls_delta - descent test fraction
122: - -tao_ssls_rho - descent test power

124:    Level: beginner
125: M*/

127: EXTERN_C_BEGIN
130: PetscErrorCode TaoCreate_SSFLS(Tao tao)
131: {
132:   TAO_SSLS       *ssls;
134:   const char     *armijo_type = TAOLINESEARCHARMIJO;

137:   PetscNewLog(tao,&ssls);
138:   tao->data = (void*)ssls;
139:   tao->ops->solve=TaoSolve_SSFLS;
140:   tao->ops->setup=TaoSetUp_SSFLS;
141:   tao->ops->view=TaoView_SSLS;
142:   tao->ops->setfromoptions=TaoSetFromOptions_SSLS;
143:   tao->ops->destroy = TaoDestroy_SSFLS;

145:   ssls->delta = 1e-10;
146:   ssls->rho = 2.1;

148:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
149:   TaoLineSearchSetType(tao->linesearch,armijo_type);
150:   TaoLineSearchSetFromOptions(tao->linesearch);
151:   /* Linesearch objective and objectivegradient routines are  set in solve routine */
152:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);

154:   tao->max_it = 2000;
155:   tao->max_funcs = 4000;
156:   tao->fatol = 0;
157:   tao->frtol = 0;
158:   tao->grtol=0;
159:   tao->grtol=0;
160: #if defined(PETSC_USE_REAL_SINGLE)
161:   tao->gatol = 1.0e-6;
162:   tao->fmin = 1.0e-4;
163: #else
164:   tao->gatol = 1.0e-16;
165:   tao->fmin = 1.0e-8;
166: #endif
167:   return(0);
168: }
169: EXTERN_C_END