Actual source code: ssils.c

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

  3: PetscErrorCode TaoSetUp_SSILS(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->ff);
 10:   VecDuplicate(tao->solution, &ssls->dpsi);
 11:   VecDuplicate(tao->solution, &ssls->da);
 12:   VecDuplicate(tao->solution, &ssls->db);
 13:   VecDuplicate(tao->solution, &ssls->t1);
 14:   VecDuplicate(tao->solution, &ssls->t2);
 15:   return 0;
 16: }

 18: PetscErrorCode TaoDestroy_SSILS(Tao tao)
 19: {
 20:   TAO_SSLS *ssls = (TAO_SSLS *)tao->data;

 22:   VecDestroy(&ssls->ff);
 23:   VecDestroy(&ssls->dpsi);
 24:   VecDestroy(&ssls->da);
 25:   VecDestroy(&ssls->db);
 26:   VecDestroy(&ssls->t1);
 27:   VecDestroy(&ssls->t2);
 28:   KSPDestroy(&tao->ksp);
 29:   PetscFree(tao->data);
 30:   return 0;
 31: }

 33: static PetscErrorCode TaoSolve_SSILS(Tao tao)
 34: {
 35:   TAO_SSLS                    *ssls = (TAO_SSLS *)tao->data;
 36:   PetscReal                    psi, ndpsi, normd, innerd, t = 0;
 37:   PetscReal                    delta, rho;
 38:   TaoLineSearchConvergedReason ls_reason;

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

 45:   TaoComputeVariableBounds(tao);
 46:   VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
 47:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, Tao_SSLS_FunctionGradient, tao);
 48:   TaoLineSearchSetObjectiveRoutine(tao->linesearch, Tao_SSLS_Function, tao);

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

 55:   tao->reason = TAO_CONTINUE_ITERATING;
 56:   while (PETSC_TRUE) {
 57:     PetscInfo(tao, "iter: %" PetscInt_FMT ", merit: %g, ndpsi: %g\n", tao->niter, (double)ssls->merit, (double)ndpsi);
 58:     /* Check the termination criteria */
 59:     TaoLogConvergenceHistory(tao, ssls->merit, ndpsi, 0.0, tao->ksp_its);
 60:     TaoMonitor(tao, tao->niter, ssls->merit, ndpsi, 0.0, t);
 61:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
 62:     if (tao->reason != TAO_CONTINUE_ITERATING) break;

 64:     /* Call general purpose update function */
 65:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
 66:     tao->niter++;

 68:     /* Calculate direction.  (Really negative of newton direction.  Therefore,
 69:        rest of the code uses -d.) */
 70:     KSPSetOperators(tao->ksp, tao->jacobian, tao->jacobian_pre);
 71:     KSPSolve(tao->ksp, ssls->ff, tao->stepdirection);
 72:     KSPGetIterationNumber(tao->ksp, &tao->ksp_its);
 73:     tao->ksp_tot_its += tao->ksp_its;
 74:     VecNorm(tao->stepdirection, NORM_2, &normd);
 75:     VecDot(tao->stepdirection, ssls->dpsi, &innerd);

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

 84:     VecScale(tao->stepdirection, -1.0);
 85:     innerd = -innerd;

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

 94: /* ---------------------------------------------------------- */
 95: /*MC
 96:    TAOSSILS - semi-smooth infeasible linesearch algorithm for solving
 97:        complementarity constraints

 99:    Options Database Keys:
100: + -tao_ssls_delta - descent test fraction
101: - -tao_ssls_rho - descent test power

103:    Level: beginner
104: M*/
105: PETSC_EXTERN PetscErrorCode TaoCreate_SSILS(Tao tao)
106: {
107:   TAO_SSLS   *ssls;
108:   const char *armijo_type = TAOLINESEARCHARMIJO;

110:   PetscNew(&ssls);
111:   tao->data                = (void *)ssls;
112:   tao->ops->solve          = TaoSolve_SSILS;
113:   tao->ops->setup          = TaoSetUp_SSILS;
114:   tao->ops->view           = TaoView_SSLS;
115:   tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
116:   tao->ops->destroy        = TaoDestroy_SSILS;

118:   ssls->delta = 1e-10;
119:   ssls->rho   = 2.1;

121:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
122:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
123:   TaoLineSearchSetType(tao->linesearch, armijo_type);
124:   TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);
125:   TaoLineSearchSetFromOptions(tao->linesearch);
126:   /* Note: linesearch objective and objectivegradient routines are set in solve routine */
127:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
128:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
129:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);

131:   /* Override default settings (unless already changed) */
132:   if (!tao->max_it_changed) tao->max_it = 2000;
133:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
134:   if (!tao->gttol_changed) tao->gttol = 0;
135:   if (!tao->grtol_changed) tao->grtol = 0;
136: #if defined(PETSC_USE_REAL_SINGLE)
137:   if (!tao->gatol_changed) tao->gatol = 1.0e-6;
138:   if (!tao->fmin_changed) tao->fmin = 1.0e-4;
139: #else
140:   if (!tao->gatol_changed) tao->gatol = 1.0e-16;
141:   if (!tao->fmin_changed) tao->fmin = 1.0e-8;
142: #endif
143:   return 0;
144: }