Actual source code: ssfls.c

      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:   PetscCall(TaoParametersInitialize(tao));
    147:   PetscObjectParameterSetDefault(tao, max_it, 2000);
    148:   PetscObjectParameterSetDefault(tao, max_funcs, 4000);
    149:   PetscObjectParameterSetDefault(tao, gttol, 0);
    150:   PetscObjectParameterSetDefault(tao, grtol, 0);
    151:   PetscObjectParameterSetDefault(tao, gatol, PetscDefined(USE_REAL_SINGLE) ? 1.0e-6 : 1.0e-16);
    152:   PetscObjectParameterSetDefault(tao, fmin, PetscDefined(USE_REAL_SINGLE) ? 1.0e-4 : 1.0e-8);
    153:   PetscFunctionReturn(PETSC_SUCCESS);
    154: }