Actual source code: asfls.c

      1: #include "src/complementarity/impls/ssls/ssls.h"
      2: // Context for ASXLS
      3: //   -- active-set        - reduced matrices formed
      4: //                          - inherit properties of original system
      5: //   -- semismooth (S)  - function not differentiable
      6: //                      - merit function continuously differentiable
      7: //                      - Fischer-Burmeister reformulation of complementarity
      8: //                        - Billups composition for two finite bounds
      9: //   -- infeasible (I)  - iterates not guaranteed to remain within bounds
     10: //   -- feasible (F)    - iterates guaranteed to remain within bounds
     11: //   -- linesearch (LS) - Armijo rule on direction
     12: //
     13: // Many other reformulations are possible and combinations of
     14: // feasible/infeasible and linesearch/trust region are possible.
     15: //
     16: // Basic theory
     17: //   Fischer-Burmeister reformulation is semismooth with a continuously
     18: //   differentiable merit function and strongly semismooth if the F has
     19: //   lipschitz continuous derivatives.
     20: //
     21: //   Every accumulation point generated by the algorithm is a stationary
     22: //   point for the merit function.  Stationary points of the merit function
     23: //   are solutions of the complementarity problem if
     24: //     a.  the stationary point has a BD-regular subdifferential, or
     25: //     b.  the Schur complement F'/F'_ff is a P_0-matrix where ff is the
     26: //         index set corresponding to the free variables.
     27: //
     28: //   If one of the accumulation points has a BD-regular subdifferential then
     29: //     a.  the entire sequence converges to this accumulation point at
     30: //         a local q-superlinear rate
     31: //     b.  if in addition the reformulation is strongly semismooth near
     32: //         this accumulation point, then the algorithm converges at a
     33: //         local q-quadratic rate.
     34: //
     35: // The theory for the feasible version follows from the feasible descent
     36: // algorithm framework.
     37: //
     38: // References:
     39: //   Billups, "Algorithms for Complementarity Problems and Generalized
     40: //     Equations," Ph.D thesis, University of Wisconsin - Madison, 1995.
     41: //   De Luca, Facchinei, Kanzow, "A Semismooth Equation Approach to the
     42: //     Solution of Nonlinear Complementarity Problems," Mathematical
     43: //     Programming, 75, pages 407-439, 1996.
     44: //   Ferris, Kanzow, Munson, "Feasible Descent Algorithms for Mixed
     45: //     Complementarity Problems," Mathematical Programming, 86,
     46: //     pages 475-497, 1999.
     47: //   Fischer, "A Special Newton-type Optimization Method," Optimization,
     48: //     24, pages 269-284, 1992
     49: //   Munson, Facchinei, Ferris, Fischer, Kanzow, "The Semismooth Algorithm
     50: //     for Large Scale Complementarity Problems," Technical Report 99-06,
     51: //     University of Wisconsin - Madison, 1999.

     54: /*------------------------------------------------------------*/
     57: static int TaoSolve_ASFLS(TAO_SOLVER tao, void *solver)
     58: {
     59:   TAO_SSLS *asls = (TAO_SSLS *)solver;
     60:   TaoVec *x, *l, *u, *ff, *dpsi, *d, *w, *t1, *t2, *da, *db;
     61:   TaoMat *J,*Jsub;
     62:   TaoVec *dxfree,*r1, *r2, *r3;
     63:   //  TaoLinearSolver *lsolver;
     64:   TaoIndexSet *FreeVariableSet,*FixedVariableSet;
     65:   double psi, ndpsi, normd, innerd, t=0;
     66:   double delta, rho;
     67:   int iter=0, info, nn, nf;
     68:   TaoTerminateReason reason;
     69:   TaoTruth flag;

     71:   TaoFunctionBegin;

     73:   // Assume that Setup has been called!
     74:   // Set the structure for the Jacobian and create a linear solver.

     76:   delta = asls->delta;
     77:   rho = asls->rho;

     79:   info = TaoGetSolution(tao, &x); CHKERRQ(info);
     80:   info = TaoGetVariableBounds(tao, &l, &u); CHKERRQ(info);
     81:   info = TaoEvaluateVariableBounds(tao,l,u); CHKERRQ(info);
     82:   info = TaoGetJacobian(tao, &J); CHKERRQ(info);

     84:   info = x->GetDimension(&nn); CHKERRQ(info);

     86:   info = x->PointwiseMaximum(x, l); CHKERRQ(info);
     87:   info = x->PointwiseMinimum(x, u); CHKERRQ(info);

     89:   info = x->CreateIndexSet(&FreeVariableSet); CHKERRQ(info);
     90:   info = x->CreateIndexSet(&FixedVariableSet); CHKERRQ(info);
     91:   info = J->CreateReducedMatrix(FreeVariableSet,FreeVariableSet,&Jsub); CHKERRQ(info);
     92:   info = x->Clone(&dxfree); CHKERRQ(info);
     93:   info = x->Clone(&r1); CHKERRQ(info);
     94:   info = x->Clone(&r2); CHKERRQ(info);
     95:   info = x->Clone(&r3); CHKERRQ(info);

     97:   //  f = asls->f;
     98:   ff = asls->ff;
     99:   dpsi = asls->dpsi;
    100:   d = asls->d;
    101:   w = asls->w;
    102:   t1 = asls->t1;
    103:   t2 = asls->t2;
    104:   da = asls->da;
    105:   db = asls->db;

    107:   info = TaoSetMeritFunction(tao, Tao_SSLS_Function, Tao_ASLS_FunctionGradient,
    108:                              TAO_NULL, TAO_NULL, asls); CHKERRQ(info);

    110:   // Calculate the function value and fischer function value at the
    111:   // current iterate

    113:   info = TaoComputeMeritFunctionGradient(tao, x, &psi, dpsi); CHKERRQ(info);
    114:   info = dpsi->Norm2(&ndpsi); CHKERRQ(info);

    116:   while (1) {

    118:     // Check the termination criteria

    120:     info = TaoMonitor(tao, iter++, asls->merit, ndpsi, 0.0, t, &reason); CHKERRQ(info);
    121:     if (TAO_CONTINUE_ITERATING != reason) break;

    123:     // We are going to solve a linear system of equations.  We need to
    124:     // set the tolerances for the solve so that we maintain an asymptotic
    125:     // rate of convergence that is superlinear.
    126:     // Note: these tolerances are for the reduced system.  We really need
    127:     // to make sure that the full system satisfies the full-space conditions.

    129:     // This rule gives superlinear asymptotic convergence
    130:     // asls->atol = TaoMin(0.5, asls->merit*sqrt(asls->merit));
    131:     // asls->rtol = 0.0;

    133:     // This rule gives quadratic asymptotic convergence
    134:     // asls->atol = TaoMin(0.5, asls->merit*asls->merit);
    135:     // asls->rtol = 0.0;

    137:     // Calculate a free and fixed set of variables.  The fixed set of
    138:     // variables are those for the d_b is approximately equal to zero.
    139:     // The definition of approximately changes as we approact the solution
    140:     // to the problem.

    142:     // No one rule is guaranteed to work in all cases.  The following
    143:     // definition is based on the norm of the Jacobian matrix.  If the
    144:     // norm is large, the tolerance becomes smaller.

    146:     info = J->Norm1(&(asls->identifier)); CHKERRQ(info);
    147:     // asls->identifier = asls->atol / (1 + asls->identifier);
    148:     asls->identifier = TaoMin(asls->merit, 1e-2) / (1 + asls->identifier);
    149:     // asls->identifier = 1e-4 / (1 + asls->identifier);
    150:     // asls->identifier = 1e-4;
    151:     // asls->identifier = 1e-8;

    153:     info = t1->SetToConstant(-asls->identifier); CHKERRQ(info);
    154:     info = t2->SetToConstant( asls->identifier); CHKERRQ(info);

    156:     info = FixedVariableSet->WhichBetweenOrEqual(t1, db, t2); CHKERRQ(info);
    157:     info = FreeVariableSet->ComplementOf(FixedVariableSet); CHKERRQ(info);

    159:     info = FixedVariableSet->GetSize(&nf);
    160:     printf("Fixed size: %d\n", nf);

    162:     // We now have our partition.  Now calculate the direction in the
    163:     // fixed variable space.  This direction is the gradient direction.

    165:     info = r1->SetReducedVec(ff, FixedVariableSet); CHKERRQ(info);
    166:     info = r2->SetReducedVec(da, FixedVariableSet); CHKERRQ(info);
    167:     info = r1->PointwiseDivide(r1, r2); CHKERRQ(info);

    169:     info = d->SetToZero(); CHKERRQ(info);
    170:     info = d->ReducedXPY(r1, FixedVariableSet);CHKERRQ(info);

    172:     // Our direction in the FixedVariableSet is fixed.  Calculate the
    173:     // information needed for the step in the FreeVariableSet.  To
    174:     // do this, we need to know the diagonal perturbation and the
    175:     // right hand side.

    177:     info = r1->SetReducedVec(da, FreeVariableSet); CHKERRQ(info);
    178:     info = r2->SetReducedVec(ff, FreeVariableSet); CHKERRQ(info);
    179:     info = r3->SetReducedVec(db, FreeVariableSet); CHKERRQ(info);

    181:     info = r1->PointwiseDivide(r1, r3); CHKERRQ(info);
    182:     info = r2->PointwiseDivide(r2, r3); CHKERRQ(info);

    184:     // r1 is the diagonal perturbation
    185:     // r2 is the right hand side
    186:     // r3 is no longer needed

    188:     // Now need to modify r2 for our direction choice in the fixed
    189:     // variable set:  calculate t1 = J*d, take the reduced vector
    190:     // of t1 and modify r2.

    192:     info = J->Multiply(d, t1); CHKERRQ(info);
    193:     info = r3->SetReducedVec(t1, FreeVariableSet); CHKERRQ(info);
    194:     info = r2->Axpy(-1.0, r3);

    196:     // Calculate the reduced problem matrix and the direction

    198:     info = Jsub->SetReducedMatrix(J, FreeVariableSet, FreeVariableSet); CHKERRQ(info);
    199:     info = Jsub->AddDiagonal(r1); CHKERRQ(info);
    200:     info = dxfree->SetReducedVec(d, FreeVariableSet);CHKERRQ(info);
    201:     info = dxfree->SetToZero(); CHKERRQ(info);

    203:     // Set the convergence for the reduced problem solver.
    204:     // info = TaoGetLinearSolver(tao, &lsolver); CHKERRQ(info);
    205:     // info = lsolver->SetTolerances(0.1*asls->rtol, 0.1*asls->atol, 1e30, 10*nn); CHKERRQ(info);

    207:     // Calculate the reduced direction.  (Really negative of Newton
    208:     // direction.  Therefore, rest of the code uses -d.)

    210:     info = TaoPreLinearSolve(tao, Jsub); CHKERRQ(info);
    211:     info = TaoLinearSolve(tao, Jsub, r2, dxfree, &flag); CHKERRQ(info);

    213:     // Add the direction in the free variables back into the real direction.

    215:     info = d->ReducedXPY(dxfree, FreeVariableSet);CHKERRQ(info);

    217:    // Calculate the norm of the relative residual in the real space

    219:     // info = t1->PointwiseMultiply(da, d); CHKERRQ(info);
    220:     // info = J->Multiply(d, t2); CHKERRQ(info);
    221:     // info = t2->PointwiseMultiply(db, t2); CHKERRQ(info);
    222:     // info = t1->Axpy( 1.0, t2); CHKERRQ(info);
    223:     // info = t1->Axpy(-1.0, ff); CHKERRQ(info);
    224:     // info = t1->Norm2(&two_norm); CHKERRQ(info);

    226:     // if (two_norm >= asls->atol) {
    227:     //   printf("Bad absolute residual: actual: %5.4e needed: %5.4e\n",
    228:     //          two_norm, asls->atol);
    229:     // }

    231:     // Check the projected real direction for descent and if not, use the
    232:     // negative gradient direction.

    234:     info = w->CopyFrom(d); CHKERRQ(info);
    235:     info = w->Negate(); CHKERRQ(info);
    236:     info = w->BoundGradientProjection(w, l, x, u); CHKERRQ(info);

    238:     info = w->Norm2(&normd); CHKERRQ(info);
    239:     info = w->Dot(dpsi, &innerd); CHKERRQ(info);

    241:     if (innerd >= -delta*pow(normd, rho)) {
    242:       info = PetscLogInfo((tao, "TaoSolve_SSLS: %d: newton direction not descent\n", iter)); CHKERRQ(info);
    243:       info = d->CopyFrom(dpsi); CHKERRQ(info);
    244:       info = d->Dot(dpsi, &innerd); CHKERRQ(info);
    245:     }
    246:     info = d->Negate(); CHKERRQ(info);
    247:     innerd = -innerd;

    249:     // We now have a correct descent direction.  Apply a linesearch to
    250:     // find the new iterate.

    252:     t = 1;
    253:     info = TaoLineSearchApply(tao, x, dpsi, d, w,
    254:                               &psi, &t, &tao->lsflag); CHKERRQ(info);
    255:     info = dpsi->Norm2(&ndpsi); CHKERRQ(info);
    256:   }

    258:   info = TaoMatDestroy(Jsub);CHKERRQ(info);
    259:   info = TaoVecDestroy(dxfree);CHKERRQ(info);
    260:   info = TaoVecDestroy(r1);CHKERRQ(info);
    261:   info = TaoVecDestroy(r2);CHKERRQ(info);
    262:   info = TaoVecDestroy(r3);CHKERRQ(info);
    263:   info = TaoIndexSetDestroy(FreeVariableSet);CHKERRQ(info);
    264:   info = TaoIndexSetDestroy(FixedVariableSet);CHKERRQ(info);

    266:   TaoFunctionReturn(0);
    267: }

    269: /* ---------------------------------------------------------- */
    273: int TaoCreate_ASFLS(TAO_SOLVER tao)
    274: {
    275:   TAO_SSLS *asls;
    276:   int        info;

    278:   TaoFunctionBegin;

    280:   info = TaoNew(TAO_SSLS, &asls); CHKERRQ(info);
    281:   info = PetscLogObjectMemory(tao, sizeof(TAO_SSLS)); CHKERRQ(info);

    283:   asls->delta = 1e-10;
    284:   asls->rho = 2.1;

    286:   asls->identifier = 1e-5;

    288:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_ASFLS,(void*)asls); CHKERRQ(info);
    289:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_SSLS,TaoSetDown_SSLS); CHKERRQ(info);
    290:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_SSLS); CHKERRQ(info);
    291:   info=TaoSetTaoViewRoutine(tao,TaoView_SSLS); CHKERRQ(info);

    293:   info = TaoCreateProjectedArmijoLineSearch(tao); CHKERRQ(info);

    295:   info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
    296:   info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);

    298:   info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
    299:   info = TaoSetGradientTolerances(tao,1.0e-16,0.0,0.0); CHKERRQ(info);
    300:   info = TaoSetFunctionLowerBound(tao,1.0e-8); CHKERRQ(info);

    302:   TaoFunctionReturn(0);
    303: }