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

  3: PetscErrorCode TaoSetUp_SSFLS(Tao tao)
  4: {
  5:   TAO_SSLS       *ssls = (TAO_SSLS *)tao->data;

  9:   VecDuplicate(tao->solution,&tao->gradient);
 10:   VecDuplicate(tao->solution,&tao->stepdirection);
 11:   VecDuplicate(tao->solution,&ssls->w);
 12:   VecDuplicate(tao->solution,&ssls->ff);
 13:   VecDuplicate(tao->solution,&ssls->dpsi);
 14:   VecDuplicate(tao->solution,&ssls->da);
 15:   VecDuplicate(tao->solution,&ssls->db);
 16:   VecDuplicate(tao->solution,&ssls->t1);
 17:   VecDuplicate(tao->solution,&ssls->t2);
 18:   if (!tao->XL) {
 19:     VecDuplicate(tao->solution,&tao->XL);
 20:   }
 21:   if (!tao->XU) {
 22:     VecDuplicate(tao->solution,&tao->XU);
 23:   }
 24:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
 25:   return(0);
 26: }

 28: static PetscErrorCode TaoSolve_SSFLS(Tao tao)
 29: {
 30:   TAO_SSLS                     *ssls = (TAO_SSLS *)tao->data;
 31:   PetscReal                    psi, ndpsi, normd, innerd, t=0;
 32:   PetscReal                    delta, rho;
 33:   TaoConvergedReason           reason;
 34:   TaoLineSearchConvergedReason ls_reason;
 35:   PetscErrorCode               ierr;

 38:   /* Assume that Setup has been called!
 39:      Set the structure for the Jacobian and create a linear solver. */
 40:   delta = ssls->delta;
 41:   rho = ssls->rho;

 43:   TaoComputeVariableBounds(tao);
 44:   /* Project solution inside bounds */
 45:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
 46:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);
 47:   TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);

 49:   /* Calculate the function value and fischer function value at the
 50:      current iterate */
 51:   TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi);
 52:   VecNorm(ssls->dpsi,NORM_2,&ndpsi);

 54:   while (PETSC_TRUE) {
 55:     ierr=PetscInfo3(tao, "iter: %D, merit: %g, ndpsi: %g\n",tao->niter, (double)ssls->merit, (double)ndpsi);
 56:     /* Check the termination criteria */
 57:     TaoMonitor(tao,tao->niter,ssls->merit,ndpsi,0.0,t,&reason);
 58:     if (reason!=TAO_CONTINUE_ITERATING) break;
 59:     tao->niter++;

 61:     /* Calculate direction.  (Really negative of newton direction.  Therefore,
 62:        rest of the code uses -d.) */
 63:     KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre);
 64:     KSPSolve(tao->ksp,ssls->ff,tao->stepdirection);
 65:     KSPGetIterationNumber(tao->ksp,&tao->ksp_its);
 66:     tao->ksp_tot_its+=tao->ksp_its;

 68:     VecCopy(tao->stepdirection,ssls->w);
 69:     VecScale(ssls->w,-1.0);
 70:     VecBoundGradientProjection(ssls->w,tao->solution,tao->XL,tao->XU,ssls->w);

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

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

 82:     VecScale(tao->stepdirection, -1.0);
 83:     innerd = -innerd;

 85:     TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
 86:     TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason);
 87:     VecNorm(ssls->dpsi,NORM_2,&ndpsi);
 88:   }
 89:   return(0);
 90: }

 92: PetscErrorCode TaoDestroy_SSFLS(Tao tao)
 93: {
 94:   TAO_SSLS       *ssls = (TAO_SSLS *)tao->data;

 98:   VecDestroy(&ssls->ff);
 99:   VecDestroy(&ssls->w);
100:   VecDestroy(&ssls->dpsi);
101:   VecDestroy(&ssls->da);
102:   VecDestroy(&ssls->db);
103:   VecDestroy(&ssls->t1);
104:   VecDestroy(&ssls->t2);
105:   PetscFree(tao->data);
106:   return(0);
107: }

109: /* ---------------------------------------------------------- */
110: /*MC
111:    TAOSSFLS - Semi-smooth feasible linesearch algorithm for solving
112:        complementarity constraints

114:    Options Database Keys:
115: + -tao_ssls_delta - descent test fraction
116: - -tao_ssls_rho - descent test power

118:    Level: beginner
119: M*/

121: PETSC_EXTERN PetscErrorCode TaoCreate_SSFLS(Tao tao)
122: {
123:   TAO_SSLS       *ssls;
125:   const char     *armijo_type = TAOLINESEARCHARMIJO;

128:   PetscNewLog(tao,&ssls);
129:   tao->data = (void*)ssls;
130:   tao->ops->solve=TaoSolve_SSFLS;
131:   tao->ops->setup=TaoSetUp_SSFLS;
132:   tao->ops->view=TaoView_SSLS;
133:   tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
134:   tao->ops->destroy = TaoDestroy_SSFLS;

136:   ssls->delta = 1e-10;
137:   ssls->rho = 2.1;

139:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
140:   TaoLineSearchSetType(tao->linesearch,armijo_type);
141:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
142:   TaoLineSearchSetFromOptions(tao->linesearch);
143:   /* Linesearch objective and objectivegradient routines are  set in solve routine */
144:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
145:   KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);

147:   /* Override default settings (unless already changed) */
148:   if (!tao->max_it_changed) tao->max_it = 2000;
149:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
150:   if (!tao->gttol_changed) tao->gttol = 0;
151:   if (!tao->grtol_changed) tao->grtol = 0;
152: #if defined(PETSC_USE_REAL_SINGLE)
153:   if (!tao->gatol_changed) tao->gatol = 1.0e-6;
154:   if (!tao->fmin_changed)  tao->fmin = 1.0e-4;
155: #else
156:   if (!tao->gatol_changed) tao->gatol = 1.0e-16;
157:   if (!tao->fmin_changed)  tao->fmin = 1.0e-8;
158: #endif
159:   return(0);
160: }