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

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

  7:   PetscFunctionBegin;
  8:   PetscCall(VecDuplicate(tao->solution, &tao->gradient));
  9:   PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
 10:   PetscCall(VecDuplicate(tao->solution, &ssls->w));
 11:   PetscCall(VecDuplicate(tao->solution, &ssls->ff));
 12:   PetscCall(VecDuplicate(tao->solution, &ssls->dpsi));
 13:   PetscCall(VecDuplicate(tao->solution, &ssls->da));
 14:   PetscCall(VecDuplicate(tao->solution, &ssls->db));
 15:   PetscCall(VecDuplicate(tao->solution, &ssls->t1));
 16:   PetscCall(VecDuplicate(tao->solution, &ssls->t2));
 17:   PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));
 18:   PetscFunctionReturn(PETSC_SUCCESS);
 19: }

 21: static PetscErrorCode TaoSolve_SSFLS(Tao tao)
 22: {
 23:   TAO_SSLS                    *ssls = (TAO_SSLS *)tao->data;
 24:   PetscReal                    psi, ndpsi, normd, innerd, t = 0;
 25:   PetscReal                    delta, rho;
 26:   TaoLineSearchConvergedReason ls_reason;

 28:   PetscFunctionBegin;
 29:   /* Assume that Setup has been called!
 30:      Set the structure for the Jacobian and create a linear solver. */
 31:   delta = ssls->delta;
 32:   rho   = ssls->rho;

 34:   PetscCall(TaoComputeVariableBounds(tao));
 35:   /* Project solution inside bounds */
 36:   PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
 37:   PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, Tao_SSLS_FunctionGradient, tao));
 38:   PetscCall(TaoLineSearchSetObjectiveRoutine(tao->linesearch, Tao_SSLS_Function, tao));

 40:   /* Calculate the function value and fischer function value at the
 41:      current iterate */
 42:   PetscCall(TaoLineSearchComputeObjectiveAndGradient(tao->linesearch, tao->solution, &psi, ssls->dpsi));
 43:   PetscCall(VecNorm(ssls->dpsi, NORM_2, &ndpsi));

 45:   tao->reason = TAO_CONTINUE_ITERATING;
 46:   while (PETSC_TRUE) {
 47:     PetscCall(PetscInfo(tao, "iter: %" PetscInt_FMT ", merit: %g, ndpsi: %g\n", tao->niter, (double)ssls->merit, (double)ndpsi));
 48:     /* Check the termination criteria */
 49:     PetscCall(TaoLogConvergenceHistory(tao, ssls->merit, ndpsi, 0.0, tao->ksp_its));
 50:     PetscCall(TaoMonitor(tao, tao->niter, ssls->merit, ndpsi, 0.0, t));
 51:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
 52:     if (tao->reason != TAO_CONTINUE_ITERATING) break;

 54:     /* Call general purpose update function */
 55:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
 56:     tao->niter++;

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

 65:     PetscCall(VecCopy(tao->stepdirection, ssls->w));
 66:     PetscCall(VecScale(ssls->w, -1.0));
 67:     PetscCall(VecBoundGradientProjection(ssls->w, tao->solution, tao->XL, tao->XU, ssls->w));

 69:     PetscCall(VecNorm(ssls->w, NORM_2, &normd));
 70:     PetscCall(VecDot(ssls->w, ssls->dpsi, &innerd));

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

 79:     PetscCall(VecScale(tao->stepdirection, -1.0));
 80:     innerd = -innerd;

 82:     PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
 83:     PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &psi, ssls->dpsi, tao->stepdirection, &t, &ls_reason));
 84:     PetscCall(VecNorm(ssls->dpsi, NORM_2, &ndpsi));
 85:   }
 86:   PetscFunctionReturn(PETSC_SUCCESS);
 87: }

 89: static PetscErrorCode TaoDestroy_SSFLS(Tao tao)
 90: {
 91:   TAO_SSLS *ssls = (TAO_SSLS *)tao->data;

 93:   PetscFunctionBegin;
 94:   PetscCall(VecDestroy(&ssls->ff));
 95:   PetscCall(VecDestroy(&ssls->w));
 96:   PetscCall(VecDestroy(&ssls->dpsi));
 97:   PetscCall(VecDestroy(&ssls->da));
 98:   PetscCall(VecDestroy(&ssls->db));
 99:   PetscCall(VecDestroy(&ssls->t1));
100:   PetscCall(VecDestroy(&ssls->t2));
101:   PetscCall(KSPDestroy(&tao->ksp));
102:   PetscCall(PetscFree(tao->data));
103:   PetscFunctionReturn(PETSC_SUCCESS);
104: }

106: /* ---------------------------------------------------------- */
107: /*MC
108:    TAOSSFLS - Semi-smooth feasible linesearch algorithm for solving
109:        complementarity constraints

111:    Options Database Keys:
112: + -tao_ssls_delta - descent test fraction
113: - -tao_ssls_rho - descent test power

115:    Level: beginner
116: M*/

118: PETSC_EXTERN PetscErrorCode TaoCreate_SSFLS(Tao tao)
119: {
120:   TAO_SSLS   *ssls;
121:   const char *armijo_type = TAOLINESEARCHARMIJO;

123:   PetscFunctionBegin;
124:   PetscCall(PetscNew(&ssls));
125:   tao->data                = (void *)ssls;
126:   tao->ops->solve          = TaoSolve_SSFLS;
127:   tao->ops->setup          = TaoSetUp_SSFLS;
128:   tao->ops->view           = TaoView_SSLS;
129:   tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
130:   tao->ops->destroy        = TaoDestroy_SSFLS;

132:   ssls->delta = 1e-10;
133:   ssls->rho   = 2.1;

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

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