Actual source code: ssfls.c

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

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

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

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

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

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

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

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

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

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

 63:     VecCopy(tao->stepdirection, ssls->w);
 64:     VecScale(ssls->w, -1.0);
 65:     VecBoundGradientProjection(ssls->w, tao->solution, tao->XL, tao->XU, ssls->w);

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

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

 77:     VecScale(tao->stepdirection, -1.0);
 78:     innerd = -innerd;

 80:     TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
 81:     TaoLineSearchApply(tao->linesearch, tao->solution, &psi, ssls->dpsi, tao->stepdirection, &t, &ls_reason);
 82:     VecNorm(ssls->dpsi, NORM_2, &ndpsi);
 83:   }
 84:   return 0;
 85: }

 87: PetscErrorCode TaoDestroy_SSFLS(Tao tao)
 88: {
 89:   TAO_SSLS *ssls = (TAO_SSLS *)tao->data;

 91:   VecDestroy(&ssls->ff);
 92:   VecDestroy(&ssls->w);
 93:   VecDestroy(&ssls->dpsi);
 94:   VecDestroy(&ssls->da);
 95:   VecDestroy(&ssls->db);
 96:   VecDestroy(&ssls->t1);
 97:   VecDestroy(&ssls->t2);
 98:   KSPDestroy(&tao->ksp);
 99:   PetscFree(tao->data);
100:   return 0;
101: }

103: /* ---------------------------------------------------------- */
104: /*MC
105:    TAOSSFLS - Semi-smooth feasible linesearch algorithm for solving
106:        complementarity constraints

108:    Options Database Keys:
109: + -tao_ssls_delta - descent test fraction
110: - -tao_ssls_rho - descent test power

112:    Level: beginner
113: M*/

115: PETSC_EXTERN PetscErrorCode TaoCreate_SSFLS(Tao tao)
116: {
117:   TAO_SSLS   *ssls;
118:   const char *armijo_type = TAOLINESEARCHARMIJO;

120:   PetscNew(&ssls);
121:   tao->data                = (void *)ssls;
122:   tao->ops->solve          = TaoSolve_SSFLS;
123:   tao->ops->setup          = TaoSetUp_SSFLS;
124:   tao->ops->view           = TaoView_SSLS;
125:   tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
126:   tao->ops->destroy        = TaoDestroy_SSFLS;

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

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

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