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:   PetscFree(tao->data);
 29:   return 0;
 30: }

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

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

 44:   TaoComputeVariableBounds(tao);
 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:   tao->reason = TAO_CONTINUE_ITERATING;
 55:   while (PETSC_TRUE) {
 56:     PetscInfo(tao, "iter: %D, merit: %g, ndpsi: %g\n",tao->niter, (double)ssls->merit, (double)ndpsi);
 57:     /* Check the termination criteria */
 58:     TaoLogConvergenceHistory(tao,ssls->merit,ndpsi,0.0,tao->ksp_its);
 59:     TaoMonitor(tao,tao->niter,ssls->merit,ndpsi,0.0,t);
 60:     (*tao->ops->convergencetest)(tao,tao->cnvP);
 61:     if (tao->reason!=TAO_CONTINUE_ITERATING) break;

 63:     /* Call general purpose update function */
 64:     if (tao->ops->update) {
 65:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
 66:     }
 67:     tao->niter++;

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

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

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

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

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

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

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

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

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

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

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