Actual source code: ssils.c

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

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

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

     23:   PetscFunctionBegin;
     24:   PetscCall(VecDestroy(&ssls->ff));
     25:   PetscCall(VecDestroy(&ssls->dpsi));
     26:   PetscCall(VecDestroy(&ssls->da));
     27:   PetscCall(VecDestroy(&ssls->db));
     28:   PetscCall(VecDestroy(&ssls->t1));
     29:   PetscCall(VecDestroy(&ssls->t2));
     30:   PetscCall(KSPDestroy(&tao->ksp));
     31:   PetscCall(PetscFree(tao->data));
     32:   PetscFunctionReturn(PETSC_SUCCESS);
     33: }

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

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

     48:   PetscCall(TaoComputeVariableBounds(tao));
     49:   PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
     50:   PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, Tao_SSLS_FunctionGradient, tao));
     51:   PetscCall(TaoLineSearchSetObjectiveRoutine(tao->linesearch, Tao_SSLS_Function, tao));

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

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

     67:     /* Call general purpose update function */
     68:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
     69:     tao->niter++;

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

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

     87:     PetscCall(VecScale(tao->stepdirection, -1.0));
     88:     innerd = -innerd;

     90:     PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
     91:     PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &psi, ssls->dpsi, tao->stepdirection, &t, &ls_reason));
     92:     PetscCall(VecNorm(ssls->dpsi, NORM_2, &ndpsi));
     93:   }
     94:   PetscFunctionReturn(PETSC_SUCCESS);
     95: }

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

    102:    Options Database Keys:
    103: + -tao_ssls_delta - descent test fraction
    104: - -tao_ssls_rho - descent test power

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

    113:   PetscFunctionBegin;
    114:   PetscCall(PetscNew(&ssls));
    115:   tao->data                = (void *)ssls;
    116:   tao->ops->solve          = TaoSolve_SSILS;
    117:   tao->ops->setup          = TaoSetUp_SSILS;
    118:   tao->ops->view           = TaoView_SSLS;
    119:   tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
    120:   tao->ops->destroy        = TaoDestroy_SSILS;

    122:   ssls->delta = 1e-10;
    123:   ssls->rho   = 2.1;

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

    135:   /* Override default settings (unless already changed) */
    136:   PetscCall(TaoParametersInitialize(tao));
    137:   PetscObjectParameterSetDefault(tao, max_it, 2000);
    138:   PetscObjectParameterSetDefault(tao, max_funcs, 4000);
    139:   PetscObjectParameterSetDefault(tao, gttol, 0);
    140:   PetscObjectParameterSetDefault(tao, grtol, 0);
    141:   PetscObjectParameterSetDefault(tao, gatol, PetscDefined(USE_REAL_SINGLE) ? 1.0e-6 : 1.0e-16);
    142:   PetscObjectParameterSetDefault(tao, fmin, PetscDefined(USE_REAL_SINGLE) ? 1.0e-4 : 1.0e-8);
    143:   PetscFunctionReturn(PETSC_SUCCESS);
    144: }