Chapter 7.  Advanced Solver Uses

CBC uses a generic OsiSolverInterface and its resolve capability. This does not give much flexibility so advanced users can inherit from their interface of choice. This section illustrates how to implement a specialized solver for a long thin problem, e.g., fast0507 again. As with the other examples in the Guide, the sample code is not guaranteed to be the fastest way to solve the problem. The main purpose of the example is to illustrate techniques. The full source is in CbcSolver2.hpp and CbcSolver2.cpp located in the CBC Samples directory, see Chapter 8, More Samples .

The method initialSolve is called a few times in CBC, and provides a convenient starting point. The modelPtr_ derives from OsiClpSolverInterface.

The resolve() method is more complicated than initialSolve(). The main pieces of data are a counter count_ which is incremented each solve and an integer array node_ which stores the last time a variable was active in a solution. For the first few solves, the normal Dual Simplex is called and node_ array is updated.

After the first few solves, only those variables which took part in a solution in the last so many solves are used. As fast0507 is a set covering problem, any rows which are already covered can be taken out.

Example 7.3. Create Small Sub-Problem

    int * whichRow = new int[numberRows]; // Array to say which rows used
    int * whichColumn = new int [numberColumns]; // Array to say which columns used
    int i;
    const double * lower = modelPtr_->columnLower();
    const double * upper = modelPtr_->columnUpper();
    setBasis(basis_,modelPtr_); // Set basis
    int nNewCol=0; // Number of columns in small model
    // Column copy of matrix
    const double * element = modelPtr_->matrix()->getElements();
    const int * row = modelPtr_->matrix()->getIndices();
    const CoinBigIndex * columnStart = modelPtr_->matrix()->getVectorStarts();
    const int * columnLength = modelPtr_->matrix()->getVectorLengths();
    int * rowActivity = new int[numberRows]; // Number of columns with entries in each row
    int * rowActivity2 = new int[numberRows]; // Lower bound on row activity for each row
    char * mark = (char *) modelPtr_->dualColumnSolution(); // Get some space to mark columns
    for (i=0;i<numberColumns;i++) {
      bool choose = (node_[i]>count_-memory_&&node_[i]>0); // Choose if used recently
      // Take if used recently or active in some sense
      if ((choose&&upper[i])
	  ||lower[i]>0.0) {
        mark[i]=1; // mark as used
	whichColumn[nNewCol++]=i; // add to list
        CoinBigIndex j;
        double value = upper[i];
        if (value) {
          for (j=columnStart[i];
               j<columnStart[i]+columnLength[i];j++) {
            int iRow=row[j];
            assert (element[j]==1.0);
            rowActivity[iRow] ++; // This variable can cover this row
          if (lower[i]>0.0) {
            for (j=columnStart[i];
                 j<columnStart[i]+columnLength[i];j++) {
              int iRow=row[j];
              rowActivity2[iRow] ++; // This row redundant
    int nOK=0; // Use to count rows which can be covered
    int nNewRow=0; // Use to make list of rows needed
    for (i=0;i<numberRows;i++) {
      if (rowActivity[i])
      if (!rowActivity2[i])
        whichRow[nNewRow++]=i; // not satisfied
        modelPtr_->setRowStatus(i,ClpSimplex::basic); // make slack basic
    if (nOK<numberRows) {
      // The variables we have do not cover rows - see if we can find any that do
      for (i=0;i<numberColumns;i++) {
        if (!mark[i]&&upper[i]) {
          CoinBigIndex j;
          int good=0;
          for (j=columnStart[i];
               j<columnStart[i]+columnLength[i];j++) {
            int iRow=row[j];
            if (!rowActivity[iRow]) {
              rowActivity[iRow] ++;
          if (good) {
            nOK+=good; // This covers - put in list
    delete [] rowActivity;
    delete [] rowActivity2;
    if (nOK<numberRows) {
      // By inspection the problem is infeasible - no need to solve
      delete [] whichRow;
      delete [] whichColumn;
      printf("infeasible by inspection\n");
    // Now make up a small model with the right rows and columns
    ClpSimplex *  temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn);

If the variables cover the rows, then the problem is feasible (no cuts are being used). (If the rows were equality constraints, then this might not be the case. More work would be needed.) After the solution to the subproblem, the reduced costs of the full problem are checked. If the reduced cost of any variable not in the subproblem is negative, the code goes back to the full problem and cleans up with Primal Simplex.

Example 7.4. Check Optimal Solution

    temp->setDualObjectiveLimit(1.0e50); // Switch off dual cutoff as problem is restricted
    temp->dual(); // solve
    double * solution = modelPtr_->primalColumnSolution(); // put back solution
    const double * solution2 = temp->primalColumnSolution();
    for (i=0;i<nNewCol;i++) {
      int iColumn = whichColumn[i];
    double * rowSolution = modelPtr_->primalRowSolution();
    const double * rowSolution2 = temp->primalRowSolution();
    double * dual = modelPtr_->dualRowSolution();
    const double * dual2 = temp->dualRowSolution();
    for (i=0;i<nNewRow;i++) {
      int iRow=whichRow[i];
    // See if optimal
    double * dj = modelPtr_->dualColumnSolution();
    // get reduced cost for large problem
    // this assumes minimization
    int nBad=0;
    for (i=0;i<numberColumns;i++) {
      if (modelPtr_->getStatus(i)==ClpSimplex::atLowerBound
    // If necessary clean up with primal (and save some statistics)
    if (nBad) {
      iterationsBad_ += modelPtr_->numberIterations();

The array node_ is updated, as for the first few solves. To give some idea of the effect of this tactic, the problem fast0507 has 63,009 variables but the small problem never has more than 4,000 variables. In only about ten percent of solves was it necessary to resolve, and then the average number of iterations on full problem was less than 20.