/* $Id: ClpGubDynamicMatrix.cpp 1525 2010-02-26 17:27:59Z mjs $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.


#include <cstdio>

#include "CoinPragma.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"

#include "ClpSimplex.hpp"
#include "ClpFactorization.hpp"
#include "ClpQuadraticObjective.hpp"
#include "ClpNonLinearCost.hpp"
// at end to get min/max!
#include "ClpGubDynamicMatrix.hpp"
#include "ClpMessage.hpp"
//#define CLP_DEBUG
//#define CLP_DEBUG_PRINT
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpGubDynamicMatrix::ClpGubDynamicMatrix ()
     : ClpGubMatrix(),
       objectiveOffset_(0.0),
       startColumn_(NULL),
       row_(NULL),
       element_(NULL),
       cost_(NULL),
       fullStart_(NULL),
       id_(NULL),
       dynamicStatus_(NULL),
       lowerColumn_(NULL),
       upperColumn_(NULL),
       lowerSet_(NULL),
       upperSet_(NULL),
       numberGubColumns_(0),
       firstAvailable_(0),
       savedFirstAvailable_(0),
       firstDynamic_(0),
       lastDynamic_(0),
       numberElements_(0)
{
     setType(13);
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpGubDynamicMatrix::ClpGubDynamicMatrix (const ClpGubDynamicMatrix & rhs)
     : ClpGubMatrix(rhs)
{
     objectiveOffset_ = rhs.objectiveOffset_;
     numberGubColumns_ = rhs.numberGubColumns_;
     firstAvailable_ = rhs.firstAvailable_;
     savedFirstAvailable_ = rhs.savedFirstAvailable_;
     firstDynamic_ = rhs.firstDynamic_;
     lastDynamic_ = rhs.lastDynamic_;
     numberElements_ = rhs.numberElements_;
     startColumn_ = ClpCopyOfArray(rhs.startColumn_, numberGubColumns_ + 1);
     CoinBigIndex numberElements = startColumn_[numberGubColumns_];
     row_ = ClpCopyOfArray(rhs.row_, numberElements);;
     element_ = ClpCopyOfArray(rhs.element_, numberElements);;
     cost_ = ClpCopyOfArray(rhs.cost_, numberGubColumns_);
     fullStart_ = ClpCopyOfArray(rhs.fullStart_, numberSets_ + 1);
     id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
     lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_, numberGubColumns_);
     upperColumn_ = ClpCopyOfArray(rhs.upperColumn_, numberGubColumns_);
     dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, numberGubColumns_);
     lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
     upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
}

/* This is the real constructor*/
ClpGubDynamicMatrix::ClpGubDynamicMatrix(ClpSimplex * model, int numberSets,
          int numberGubColumns, const int * starts,
          const double * lower, const double * upper,
          const CoinBigIndex * startColumn, const int * row,
          const double * element, const double * cost,
          const double * lowerColumn, const double * upperColumn,
          const unsigned char * status)
     : ClpGubMatrix()
{
     objectiveOffset_ = model->objectiveOffset();
     model_ = model;
     numberSets_ = numberSets;
     numberGubColumns_ = numberGubColumns;
     fullStart_ = ClpCopyOfArray(starts, numberSets_ + 1);
     lower_ = ClpCopyOfArray(lower, numberSets_);
     upper_ = ClpCopyOfArray(upper, numberSets_);
     int numberColumns = model->numberColumns();
     int numberRows = model->numberRows();
     // Number of columns needed
     int numberGubInSmall = numberSets_ + numberRows + 2 * model->factorizationFrequency() + 2;
     // for small problems this could be too big
     //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
     int numberNeeded = numberGubInSmall + numberColumns;
     firstAvailable_ = numberColumns;
     savedFirstAvailable_ = numberColumns;
     firstDynamic_ = numberColumns;
     lastDynamic_ = numberNeeded;
     startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1);
     CoinBigIndex numberElements = startColumn_[numberGubColumns_];
     row_ = ClpCopyOfArray(row, numberElements);
     element_ = new double[numberElements];
     CoinBigIndex i;
     for (i = 0; i < numberElements; i++)
          element_[i] = element[i];
     cost_ = new double[numberGubColumns_];
     for (i = 0; i < numberGubColumns_; i++) {
          cost_[i] = cost[i];
          // need sorted
          CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i+1], element_ + startColumn_[i]);
     }
     if (lowerColumn) {
          lowerColumn_ = new double[numberGubColumns_];
          for (i = 0; i < numberGubColumns_; i++)
               lowerColumn_[i] = lowerColumn[i];
     } else {
          lowerColumn_ = NULL;
     }
     if (upperColumn) {
          upperColumn_ = new double[numberGubColumns_];
          for (i = 0; i < numberGubColumns_; i++)
               upperColumn_[i] = upperColumn[i];
     } else {
          upperColumn_ = NULL;
     }
     if (upperColumn || lowerColumn) {
          lowerSet_ = new double[numberSets_];
          for (i = 0; i < numberSets_; i++) {
               if (lower[i] > -1.0e20)
                    lowerSet_[i] = lower[i];
               else
                    lowerSet_[i] = -1.0e30;
          }
          upperSet_ = new double[numberSets_];
          for (i = 0; i < numberSets_; i++) {
               if (upper[i] < 1.0e20)
                    upperSet_[i] = upper[i];
               else
                    upperSet_[i] = 1.0e30;
          }
     } else {
          lowerSet_ = NULL;
          upperSet_ = NULL;
     }
     start_ = NULL;
     end_ = NULL;
     dynamicStatus_ = NULL;
     id_ = new int[numberGubInSmall];
     for (i = 0; i < numberGubInSmall; i++)
          id_[i] = -1;
     ClpPackedMatrix* originalMatrixA =
          dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
     assert (originalMatrixA);
     CoinPackedMatrix * originalMatrix = originalMatrixA->getPackedMatrix();
     originalMatrixA->setMatrixNull(); // so can be deleted safely
     // guess how much space needed
     double guess = originalMatrix->getNumElements() + 10;
     guess /= static_cast<double> (numberColumns);
     guess *= 2 * numberGubColumns_;
     numberElements_ = static_cast<int> (CoinMin(guess, 10000000.0));
     numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements();
     matrix_ = originalMatrix;
     flags_ &= ~1;
     // resize model (matrix stays same)
     model->resize(numberRows, numberNeeded);
     if (upperColumn_) {
          // set all upper bounds so we have enough space
          double * columnUpper = model->columnUpper();
          for(i = firstDynamic_; i < lastDynamic_; i++)
               columnUpper[i] = 1.0e10;
     }
     // resize matrix
     // extra 1 is so can keep number of elements handy
     originalMatrix->reserve(numberNeeded, numberElements_, true);
     originalMatrix->reserve(numberNeeded + 1, numberElements_, false);
     originalMatrix->getMutableVectorStarts()[numberColumns] = originalMatrix->getNumElements();
     // redo number of columns
     numberColumns = matrix_->getNumCols();
     backward_ = new int[numberNeeded];
     backToPivotRow_ = new int[numberNeeded];
     // We know a bit better
     delete [] changeCost_;
     changeCost_ = new double [numberRows+numberSets_];
     keyVariable_ = new int[numberSets_];
     // signal to need new ordering
     next_ = NULL;
     for (int iColumn = 0; iColumn < numberNeeded; iColumn++)
          backward_[iColumn] = -1;

     firstGub_ = firstDynamic_;
     lastGub_ = lastDynamic_;
     if (!lowerColumn_ && !upperColumn_)
          gubType_ = 8;
     if (status) {
          status_ = ClpCopyOfArray(status, numberSets_);
     } else {
          status_ = new unsigned char [numberSets_];
          memset(status_, 0, numberSets_);
          int i;
          for (i = 0; i < numberSets_; i++) {
               // make slack key
               setStatus(i, ClpSimplex::basic);
          }
     }
     saveStatus_ = new unsigned char [numberSets_];
     memset(saveStatus_, 0, numberSets_);
     savedKeyVariable_ = new int [numberSets_];
     memset(savedKeyVariable_, 0, numberSets_ * sizeof(int));
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpGubDynamicMatrix::~ClpGubDynamicMatrix ()
{
     delete [] startColumn_;
     delete [] row_;
     delete [] element_;
     delete [] cost_;
     delete [] fullStart_;
     delete [] id_;
     delete [] dynamicStatus_;
     delete [] lowerColumn_;
     delete [] upperColumn_;
     delete [] lowerSet_;
     delete [] upperSet_;
}

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpGubDynamicMatrix &
ClpGubDynamicMatrix::operator=(const ClpGubDynamicMatrix& rhs)
{
     if (this != &rhs) {
          ClpGubMatrix::operator=(rhs);
          delete [] startColumn_;
          delete [] row_;
          delete [] element_;
          delete [] cost_;
          delete [] fullStart_;
          delete [] id_;
          delete [] dynamicStatus_;
          delete [] lowerColumn_;
          delete [] upperColumn_;
          delete [] lowerSet_;
          delete [] upperSet_;
          objectiveOffset_ = rhs.objectiveOffset_;
          numberGubColumns_ = rhs.numberGubColumns_;
          firstAvailable_ = rhs.firstAvailable_;
          savedFirstAvailable_ = rhs.savedFirstAvailable_;
          firstDynamic_ = rhs.firstDynamic_;
          lastDynamic_ = rhs.lastDynamic_;
          numberElements_ = rhs.numberElements_;
          startColumn_ = ClpCopyOfArray(rhs.startColumn_, numberGubColumns_ + 1);
          int numberElements = startColumn_[numberGubColumns_];
          row_ = ClpCopyOfArray(rhs.row_, numberElements);;
          element_ = ClpCopyOfArray(rhs.element_, numberElements);;
          cost_ = ClpCopyOfArray(rhs.cost_, numberGubColumns_);
          fullStart_ = ClpCopyOfArray(rhs.fullStart_, numberSets_ + 1);
          id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
          lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_, numberGubColumns_);
          upperColumn_ = ClpCopyOfArray(rhs.upperColumn_, numberGubColumns_);
          dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, numberGubColumns_);
          lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
          upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
     }
     return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpMatrixBase * ClpGubDynamicMatrix::clone() const
{
     return new ClpGubDynamicMatrix(*this);
}
// Partial pricing
void
ClpGubDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
                                    int & bestSequence, int & numberWanted)
{
     assert(!model->rowScale());
     numberWanted = currentWanted_;
     if (!numberSets_) {
          // no gub
          ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
          return;
     } else {
          // and do some proportion of full set
          int startG2 = static_cast<int> (startFraction * numberSets_);
          int endG2 = static_cast<int> (endFraction * numberSets_ + 0.1);
          endG2 = CoinMin(endG2, numberSets_);
          //printf("gub price - set start %d end %d\n",
          //   startG2,endG2);
          double tolerance = model->currentDualTolerance();
          double * reducedCost = model->djRegion();
          const double * duals = model->dualRowSolution();
          double * cost = model->costRegion();
          double bestDj;
          int numberRows = model->numberRows();
          int numberColumns = lastDynamic_;
          // If nothing found yet can go all the way to end
          int endAll = endG2;
          if (bestSequence < 0 && !startG2)
               endAll = numberSets_;
          if (bestSequence >= 0)
               bestDj = fabs(reducedCost[bestSequence]);
          else
               bestDj = tolerance;
          int saveSequence = bestSequence;
          double djMod = 0.0;
          double infeasibilityCost = model->infeasibilityCost();
          double bestDjMod = 0.0;
          //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
          //     startG2,endG2,numberWanted);
          int bestType = -1;
          int bestSet = -1;
          const double * element = matrix_->getElements();
          const int * row = matrix_->getIndices();
          const CoinBigIndex * startColumn = matrix_->getVectorStarts();
          int * length = matrix_->getMutableVectorLengths();
#if 0
          // make sure first available is clean (in case last iteration rejected)
          cost[firstAvailable_] = 0.0;
          length[firstAvailable_] = 0;
          model->nonLinearCost()->setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0);
          model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
          {
               for (int i = firstAvailable_; i < lastDynamic_; i++)
                    assert(!cost[i]);
          }
#endif
#ifdef CLP_DEBUG
          {
               for (int i = firstDynamic_; i < firstAvailable_; i++) {
                    assert (getDynamicStatus(id_[i-firstDynamic_]) == inSmall);
               }
          }
#endif
          int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
          int minNeg = minimumGoodReducedCosts_ < 0 ? 5 : minimumGoodReducedCosts_;
          for (int iSet = startG2; iSet < endAll; iSet++) {
               if (numberWanted + minNeg < originalWanted_ && iSet > startG2 + minSet) {
                    // give up
                    numberWanted = 0;
                    break;
               } else if (iSet == endG2 && bestSequence >= 0) {
                    break;
               }
               CoinBigIndex j;
               int iBasic = keyVariable_[iSet];
               if (iBasic >= numberColumns) {
                    djMod = - weight(iSet) * infeasibilityCost;
               } else {
                    // get dj without
                    assert (model->getStatus(iBasic) == ClpSimplex::basic);
                    djMod = 0.0;

                    for (j = startColumn[iBasic];
                              j < startColumn[iBasic] + length[iBasic]; j++) {
                         int jRow = row[j];
                         djMod -= duals[jRow] * element[j];
                    }
                    djMod += cost[iBasic];
                    // See if gub slack possible - dj is djMod
                    if (getStatus(iSet) == ClpSimplex::atLowerBound) {
                         double value = -djMod;
                         if (value > tolerance) {
                              numberWanted--;
                              if (value > bestDj) {
                                   // check flagged variable and correct dj
                                   if (!flagged(iSet)) {
                                        bestDj = value;
                                        bestSequence = numberRows + numberColumns + iSet;
                                        bestDjMod = djMod;
                                        bestType = 0;
                                        bestSet = iSet;
                                   } else {
                                        // just to make sure we don't exit before got something
                                        numberWanted++;
                                        abort();
                                   }
                              }
                         }
                    } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
                         double value = djMod;
                         if (value > tolerance) {
                              numberWanted--;
                              if (value > bestDj) {
                                   // check flagged variable and correct dj
                                   if (!flagged(iSet)) {
                                        bestDj = value;
                                        bestSequence = numberRows + numberColumns + iSet;
                                        bestDjMod = djMod;
                                        bestType = 0;
                                        bestSet = iSet;
                                   } else {
                                        // just to make sure we don't exit before got something
                                        numberWanted++;
                                        abort();
                                   }
                              }
                         }
                    }
               }
               for (int iSequence = fullStart_[iSet]; iSequence < fullStart_[iSet+1]; iSequence++) {
                    DynamicStatus status = getDynamicStatus(iSequence);
                    if (status != inSmall) {
                         double value = cost_[iSequence] - djMod;
                         for (j = startColumn_[iSequence];
                                   j < startColumn_[iSequence+1]; j++) {
                              int jRow = row_[j];
                              value -= duals[jRow] * element_[j];
                         }
                         // change sign if at lower bound
                         if (status == atLowerBound)
                              value = -value;
                         if (value > tolerance) {
                              numberWanted--;
                              if (value > bestDj) {
                                   // check flagged variable and correct dj
                                   if (!flagged(iSequence)) {
                                        bestDj = value;
                                        bestSequence = iSequence;
                                        bestDjMod = djMod;
                                        bestType = 1;
                                        bestSet = iSet;
                                   } else {
                                        // just to make sure we don't exit before got something
                                        numberWanted++;
                                   }
                              }
                         }
                    }
               }
               if (numberWanted <= 0) {
                    numberWanted = 0;
                    break;
               }
          }
          // Do packed part before gub and small gub - but lightly
          int saveMinNeg = minimumGoodReducedCosts_;
          int saveSequence2 = bestSequence;
          if (bestSequence >= 0)
               minimumGoodReducedCosts_ = -2;
          int saveLast = lastGub_;
          lastGub_ = firstAvailable_;
          currentWanted_ = numberWanted;
          ClpGubMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
          minimumGoodReducedCosts_ = saveMinNeg;
          lastGub_ = saveLast;
          if (bestSequence != saveSequence2) {
               bestType = -1; // in normal or small gub part
               saveSequence = bestSequence;
          }
          if (bestSequence != saveSequence || bestType >= 0) {
               double * lowerColumn = model->lowerRegion();
               double * upperColumn = model->upperRegion();
               double * solution = model->solutionRegion();
               if (bestType > 0) {
                    // recompute dj and create
                    double value = cost_[bestSequence] - bestDjMod;
                    for (CoinBigIndex jBigIndex = startColumn_[bestSequence];
                              jBigIndex < startColumn_[bestSequence+1]; jBigIndex++) {
                         int jRow = row_[jBigIndex];
                         value -= duals[jRow] * element_[jBigIndex];
                    }
                    double * element =  matrix_->getMutableElements();
                    int * row = matrix_->getMutableIndices();
                    CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
                    int * length = matrix_->getMutableVectorLengths();
                    CoinBigIndex numberElements = startColumn[firstAvailable_];
                    int numberThis = startColumn_[bestSequence+1] - startColumn_[bestSequence];
                    if (numberElements + numberThis > numberElements_) {
                         // need to redo
                         numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
                         matrix_->reserve(numberColumns, numberElements_);
                         element =  matrix_->getMutableElements();
                         row = matrix_->getMutableIndices();
                         // these probably okay but be safe
                         startColumn = matrix_->getMutableVectorStarts();
                         length = matrix_->getMutableVectorLengths();
                    }
                    // already set startColumn[firstAvailable_]=numberElements;
                    length[firstAvailable_] = numberThis;
                    model->costRegion()[firstAvailable_] = cost_[bestSequence];
                    CoinBigIndex base = startColumn_[bestSequence];
                    for (int j = 0; j < numberThis; j++) {
                         row[numberElements] = row_[base+j];
                         element[numberElements++] = element_[base+j];
                    }
                    id_[firstAvailable_-firstDynamic_] = bestSequence;
                    //printf("best %d\n",bestSequence);
                    backward_[firstAvailable_] = bestSet;
                    model->solutionRegion()[firstAvailable_] = 0.0;
                    if (!lowerColumn_ && !upperColumn_) {
                         model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
                         lowerColumn[firstAvailable_] = 0.0;
                         upperColumn[firstAvailable_] = COIN_DBL_MAX;
                    }  else {
                         DynamicStatus status = getDynamicStatus(bestSequence);
                         if (lowerColumn_)
                              lowerColumn[firstAvailable_] = lowerColumn_[bestSequence];
                         else
                              lowerColumn[firstAvailable_] = 0.0;
                         if (upperColumn_)
                              upperColumn[firstAvailable_] = upperColumn_[bestSequence];
                         else
                              upperColumn[firstAvailable_] = COIN_DBL_MAX;
                         if (status == atLowerBound) {
                              solution[firstAvailable_] = lowerColumn[firstAvailable_];
                              model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
                         } else {
                              solution[firstAvailable_] = upperColumn[firstAvailable_];
                              model->setStatus(firstAvailable_, ClpSimplex::atUpperBound);
                         }
                    }
                    model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
                                                   lowerColumn[firstAvailable_],
                                                   upperColumn[firstAvailable_], cost_[bestSequence]);
                    bestSequence = firstAvailable_;
                    // firstAvailable_ only updated if good pivot (in updatePivot)
                    startColumn[firstAvailable_+1] = numberElements;
                    //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
                    reducedCost[bestSequence] = value;
                    gubSlackIn_ = -1;
               } else {
                    // slack - make last column
                    gubSlackIn_ = bestSequence - numberRows - numberColumns;
                    bestSequence = numberColumns + 2 * numberRows;
                    reducedCost[bestSequence] = bestDjMod;
                    //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
                    model->setStatus(bestSequence, getStatus(gubSlackIn_));
                    if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
                         solution[bestSequence] = upper_[gubSlackIn_];
                    else
                         solution[bestSequence] = lower_[gubSlackIn_];
                    lowerColumn[bestSequence] = lower_[gubSlackIn_];
                    upperColumn[bestSequence] = upper_[gubSlackIn_];
                    model->costRegion()[bestSequence] = 0.0;
                    model->nonLinearCost()->setOne(bestSequence, solution[bestSequence], lowerColumn[bestSequence],
                                                   upperColumn[bestSequence], 0.0);
               }
               savedBestSequence_ = bestSequence;
               savedBestDj_ = reducedCost[savedBestSequence_];
          }
          // See if may be finished
          if (!startG2 && bestSequence < 0)
               infeasibilityWeight_ = model_->infeasibilityCost();
          else if (bestSequence >= 0)
               infeasibilityWeight_ = -1.0;
     }
     currentWanted_ = numberWanted;
}
// This is local to Gub to allow synchronization when status is good
int
ClpGubDynamicMatrix::synchronize(ClpSimplex * model, int mode)
{
     int returnNumber = 0;
     switch (mode) {
     case 0: {
#ifdef CLP_DEBUG
          {
               for (int i = 0; i < numberSets_; i++)
                    assert(toIndex_[i] == -1);
          }
#endif
          // lookup array
          int * lookup = new int[lastDynamic_];
          int iColumn;
          int numberColumns = model->numberColumns();
          double * element =  matrix_->getMutableElements();
          int * row = matrix_->getMutableIndices();
          CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
          int * length = matrix_->getMutableVectorLengths();
          double * cost = model->costRegion();
          double * lowerColumn = model->lowerRegion();
          double * upperColumn = model->upperRegion();
          int * pivotVariable = model->pivotVariable();
          CoinBigIndex numberElements = startColumn[firstDynamic_];
          // first just do lookup and basic stuff
          int currentNumber = firstAvailable_;
          firstAvailable_ = firstDynamic_;
          int numberToDo = 0;
          double objectiveChange = 0.0;
          double * solution = model->solutionRegion();
          for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
               int iSet = backward_[iColumn];
               if (toIndex_[iSet] < 0) {
                    toIndex_[iSet] = 0;
                    fromIndex_[numberToDo++] = iSet;
               }
               if (model->getStatus(iColumn) == ClpSimplex::basic || iColumn == keyVariable_[iSet]) {
                    lookup[iColumn] = firstAvailable_;
                    if (iColumn != keyVariable_[iSet]) {
                         int iPivot = backToPivotRow_[iColumn];
                         backToPivotRow_[firstAvailable_] = iPivot;
                         pivotVariable[iPivot] = firstAvailable_;
                    }
                    firstAvailable_++;
               } else {
                    int jColumn = id_[iColumn-firstDynamic_];
                    setDynamicStatus(jColumn, atLowerBound);
                    if (lowerColumn_ || upperColumn_) {
                         if (model->getStatus(iColumn) == ClpSimplex::atUpperBound)
                              setDynamicStatus(jColumn, atUpperBound);
                         // treat solution as if exactly at a bound
                         double value = solution[iColumn];
                         if (fabs(value - lowerColumn[iColumn]) < fabs(value - upperColumn[iColumn]))
                              value = lowerColumn[iColumn];
                         else
                              value = upperColumn[iColumn];
                         objectiveChange += cost[iColumn] * value;
                         // redo lower and upper on sets
                         double shift = value;
                         if (lowerSet_[iSet] > -1.0e20)
                              lower_[iSet] = lowerSet_[iSet] - shift;
                         if (upperSet_[iSet] < 1.0e20)
                              upper_[iSet] = upperSet_[iSet] - shift;
                    }
                    lookup[iColumn] = -1;
               }
          }
          model->setObjectiveOffset(model->objectiveOffset() + objectiveChange);
          firstAvailable_ = firstDynamic_;
          for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
               if (lookup[iColumn] >= 0) {
                    // move
                    int jColumn = id_[iColumn-firstDynamic_];
                    id_[firstAvailable_-firstDynamic_] = jColumn;
                    int numberThis = startColumn_[jColumn+1] - startColumn_[jColumn];
                    length[firstAvailable_] = numberThis;
                    cost[firstAvailable_] = cost[iColumn];
                    lowerColumn[firstAvailable_] = lowerColumn[iColumn];
                    upperColumn[firstAvailable_] = upperColumn[iColumn];
                    double originalLower = lowerColumn_ ? lowerColumn_[jColumn] : 0.0;
                    double originalUpper = upperColumn_ ? upperColumn_[jColumn] : COIN_DBL_MAX;
                    if (originalUpper > 1.0e30)
                         originalUpper = COIN_DBL_MAX;
                    model->nonLinearCost()->setOne(firstAvailable_, solution[iColumn],
                                                   originalLower, originalUpper,
                                                   cost_[jColumn]);
                    CoinBigIndex base = startColumn_[jColumn];
                    for (int j = 0; j < numberThis; j++) {
                         row[numberElements] = row_[base+j];
                         element[numberElements++] = element_[base+j];
                    }
                    model->setStatus(firstAvailable_, model->getStatus(iColumn));
                    backward_[firstAvailable_] = backward_[iColumn];
                    solution[firstAvailable_] = solution[iColumn];
                    firstAvailable_++;
                    startColumn[firstAvailable_] = numberElements;
               }
          }
          // clean up next_
          int * temp = new int [firstAvailable_];
          for (int jSet = 0; jSet < numberToDo; jSet++) {
               int iSet = fromIndex_[jSet];
               toIndex_[iSet] = -1;
               int last = keyVariable_[iSet];
               int j = next_[last];
               bool setTemp = true;
               if (last < lastDynamic_) {
                    last = lookup[last];
                    assert (last >= 0);
                    keyVariable_[iSet] = last;
               } else if (j >= 0) {
                    int newJ = lookup[j];
                    assert (newJ >= 0);
                    j = next_[j];
                    next_[last] = newJ;
                    last = newJ;
               } else {
                    next_[last] = -(iSet + numberColumns + 1);
                    setTemp = false;
               }
               while (j >= 0) {
                    int newJ = lookup[j];
                    assert (newJ >= 0);
                    temp[last] = newJ;
                    last = newJ;
                    j = next_[j];
               }
               if (setTemp)
                    temp[last] = -(keyVariable_[iSet] + 1);
               if (lowerSet_) {
                    // we only need to get lower_ and upper_ correct
                    double shift = 0.0;
                    for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
                         if (getDynamicStatus(j) == atUpperBound)
                              shift += upperColumn_[j];
                         else if (getDynamicStatus(j) == atLowerBound && lowerColumn_)
                              shift += lowerColumn_[j];
                    if (lowerSet_[iSet] > -1.0e20)
                         lower_[iSet] = lowerSet_[iSet] - shift;
                    if (upperSet_[iSet] < 1.0e20)
                         upper_[iSet] = upperSet_[iSet] - shift;
               }
          }
          // move to next_
          CoinMemcpyN(temp + firstDynamic_, (firstAvailable_ - firstDynamic_), next_ + firstDynamic_);
          // if odd iterations may be one out so adjust currentNumber
          currentNumber = CoinMin(currentNumber + 1, lastDynamic_);
          // zero solution
          CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
          // zero cost
          CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
          // zero lengths
          CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
          for ( iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
               model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
               model->setStatus(iColumn, ClpSimplex::atLowerBound);
               backward_[iColumn] = -1;
          }
          delete [] lookup;
          delete [] temp;
          // make sure fromIndex clean
          fromIndex_[0] = -1;
          //#define CLP_DEBUG
#ifdef CLP_DEBUG
          // debug
          {
               int i;
               int numberRows = model->numberRows();
               char * xxxx = new char[numberColumns];
               memset(xxxx, 0, numberColumns);
               for (i = 0; i < numberRows; i++) {
                    int iPivot = pivotVariable[i];
                    assert (model->getStatus(iPivot) == ClpSimplex::basic);
                    if (iPivot < numberColumns && backward_[iPivot] >= 0)
                         xxxx[iPivot] = 1;
               }
               for (i = 0; i < numberSets_; i++) {
                    int key = keyVariable_[i];
                    int iColumn = next_[key];
                    int k = 0;
                    while(iColumn >= 0) {
                         k++;
                         assert (k < 100);
                         assert (backward_[iColumn] == i);
                         iColumn = next_[iColumn];
                    }
                    int stop = -(key + 1);
                    while (iColumn != stop) {
                         assert (iColumn < 0);
                         iColumn = -iColumn - 1;
                         k++;
                         assert (k < 100);
                         assert (backward_[iColumn] == i);
                         iColumn = next_[iColumn];
                    }
                    iColumn = next_[key];
                    while (iColumn >= 0) {
                         assert (xxxx[iColumn]);
                         xxxx[iColumn] = 0;
                         iColumn = next_[iColumn];
                    }
               }
               for (i = 0; i < numberColumns; i++) {
                    if (i < numberColumns && backward_[i] >= 0) {
                         assert (!xxxx[i] || i == keyVariable_[backward_[i]]);
                    }
               }
               delete [] xxxx;
          }
          {
               for (int i = 0; i < numberSets_; i++)
                    assert(toIndex_[i] == -1);
          }
#endif
          savedFirstAvailable_ = firstAvailable_;
     }
     break;
     // flag a variable
     case 1: {
          // id will be sitting at firstAvailable
          int sequence = id_[firstAvailable_-firstDynamic_];
          assert (!flagged(sequence));
          setFlagged(sequence);
          model->clearFlagged(firstAvailable_);
     }
     break;
     // unflag all variables
     case 2: {
          for (int i = 0; i < numberGubColumns_; i++) {
               if (flagged(i)) {
                    unsetFlagged(i);
                    returnNumber++;
               }
          }
     }
     break;
     //  just reset costs and bounds (primal)
     case 3: {
          double * cost = model->costRegion();
          double * solution = model->solutionRegion();
          double * lowerColumn = model->columnLower();
          double * upperColumn = model->columnUpper();
          for (int i = firstDynamic_; i < firstAvailable_; i++) {
               int jColumn = id_[i-firstDynamic_];
               cost[i] = cost_[jColumn];
               if (!lowerColumn_ && !upperColumn_) {
                    lowerColumn[i] = 0.0;
                    upperColumn[i] = COIN_DBL_MAX;
               }  else {
                    if (lowerColumn_)
                         lowerColumn[i] = lowerColumn_[jColumn];
                    else
                         lowerColumn[i] = 0.0;
                    if (upperColumn_)
                         upperColumn[i] = upperColumn_[jColumn];
                    else
                         upperColumn[i] = COIN_DBL_MAX;
               }
               if (model->nonLinearCost())
                    model->nonLinearCost()->setOne(i, solution[i],
                                                   lowerColumn[i],
                                                   upperColumn[i], cost_[jColumn]);
          }
          if (!model->numberIterations() && rhsOffset_) {
               lastRefresh_ = - refreshFrequency_; // force refresh
          }
     }
     break;
     // and get statistics for column generation
     case 4: {
          // In theory we should subtract out ones we have done but ....
          // If key slack then dual 0.0
          // If not then slack could be dual infeasible
          // dj for key is zero so that defines dual on set
          int i;
          int numberColumns = model->numberColumns();
          double * dual = model->dualRowSolution();
          double infeasibilityCost = model->infeasibilityCost();
          double dualTolerance = model->dualTolerance();
          double relaxedTolerance = dualTolerance;
          // we can't really trust infeasibilities if there is dual error
          double error = CoinMin(1.0e-2, model->largestDualError());
          // allow tolerance at least slightly bigger than standard
          relaxedTolerance = relaxedTolerance +  error;
          // but we will be using difference
          relaxedTolerance -= dualTolerance;
          double objectiveOffset = 0.0;
          for (i = 0; i < numberSets_; i++) {
               int kColumn = keyVariable_[i];
               double value = 0.0;
               if (kColumn < numberColumns) {
                    kColumn = id_[kColumn-firstDynamic_];
                    // dj without set
                    value = cost_[kColumn];
                    for (CoinBigIndex j = startColumn_[kColumn];
                              j < startColumn_[kColumn+1]; j++) {
                         int iRow = row_[j];
                         value -= dual[iRow] * element_[j];
                    }
                    double infeasibility = 0.0;
                    if (getStatus(i) == ClpSimplex::atLowerBound) {
                         if (-value > dualTolerance)
                              infeasibility = -value - dualTolerance;
                    } else if (getStatus(i) == ClpSimplex::atUpperBound) {
                         if (value > dualTolerance)
                              infeasibility = value - dualTolerance;
                    }
                    if (infeasibility > 0.0) {
                         sumDualInfeasibilities_ += infeasibility;
                         if (infeasibility > relaxedTolerance)
                              sumOfRelaxedDualInfeasibilities_ += infeasibility;
                         numberDualInfeasibilities_ ++;
                    }
               } else {
                    // slack key - may not be feasible
                    assert (getStatus(i) == ClpSimplex::basic);
                    // negative as -1.0 for slack
                    value = -weight(i) * infeasibilityCost;
               }
               // Now subtract out from all
               for (CoinBigIndex k = fullStart_[i]; k < fullStart_[i+1]; k++) {
                    if (getDynamicStatus(k) != inSmall) {
                         double djValue = cost_[k] - value;
                         for (CoinBigIndex j = startColumn_[k];
                                   j < startColumn_[k+1]; j++) {
                              int iRow = row_[j];
                              djValue -= dual[iRow] * element_[j];
                         }
                         double infeasibility = 0.0;
                         double shift = 0.0;
                         if (getDynamicStatus(k) == atLowerBound) {
                              if (lowerColumn_)
                                   shift = lowerColumn_[k];
                              if (djValue < -dualTolerance)
                                   infeasibility = -djValue - dualTolerance;
                         } else {
                              // at upper bound
                              shift = upperColumn_[k];
                              if (djValue > dualTolerance)
                                   infeasibility = djValue - dualTolerance;
                         }
                         objectiveOffset += shift * cost_[k];
                         if (infeasibility > 0.0) {
                              sumDualInfeasibilities_ += infeasibility;
                              if (infeasibility > relaxedTolerance)
                                   sumOfRelaxedDualInfeasibilities_ += infeasibility;
                              numberDualInfeasibilities_ ++;
                         }
                    }
               }
          }
          model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
     }
     break;
     // see if time to re-factorize
     case 5: {
          if (firstAvailable_ > numberSets_ + model->numberRows() + model->factorizationFrequency())
               returnNumber = 4;
     }
     break;
     // return 1 if there may be changing bounds on variable (column generation)
     case 6: {
          returnNumber = (lowerColumn_ != NULL || upperColumn_ != NULL) ? 1 : 0;
#if 0
          if (!returnNumber) {
               // may be gub slacks
               for (int i = 0; i < numberSets_; i++) {
                    if (upper_[i] > lower_[i]) {
                         returnNumber = 1;
                         break;
                    }
               }
          }
#endif
     }
     break;
     // restore firstAvailable_
     case 7: {
          int iColumn;
          int * length = matrix_->getMutableVectorLengths();
          double * cost = model->costRegion();
          double * solution = model->solutionRegion();
          int currentNumber = firstAvailable_;
          firstAvailable_ = savedFirstAvailable_;
          // zero solution
          CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
          // zero cost
          CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
          // zero lengths
          CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
          for ( iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
               model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
               model->setStatus(iColumn, ClpSimplex::atLowerBound);
               backward_[iColumn] = -1;
          }
     }
     break;
     // make sure set is clean
     case 8: {
          int sequenceIn = model->sequenceIn();
          if (sequenceIn < model->numberColumns()) {
               int iSet = backward_[sequenceIn];
               if (iSet >= 0 && lowerSet_) {
                    // we only need to get lower_ and upper_ correct
                    double shift = 0.0;
                    for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
                         if (getDynamicStatus(j) == atUpperBound)
                              shift += upperColumn_[j];
                         else if (getDynamicStatus(j) == atLowerBound && lowerColumn_)
                              shift += lowerColumn_[j];
                    if (lowerSet_[iSet] > -1.0e20)
                         lower_[iSet] = lowerSet_[iSet] - shift;
                    if (upperSet_[iSet] < 1.0e20)
                         upper_[iSet] = upperSet_[iSet] - shift;
               }
               if (sequenceIn == firstAvailable_) {
                    // not really in small problem
                    int iBig = id_[sequenceIn-firstDynamic_];
                    if (model->getStatus(sequenceIn) == ClpSimplex::atLowerBound)
                         setDynamicStatus(iBig, atLowerBound);
                    else
                         setDynamicStatus(iBig, atUpperBound);
               }
          }
     }
     break;
     // adjust lower,upper
     case 9: {
          int sequenceIn = model->sequenceIn();
          if (sequenceIn >= firstDynamic_ && sequenceIn < lastDynamic_ && lowerSet_) {
               int iSet = backward_[sequenceIn];
               assert (iSet >= 0);
               int inBig = id_[sequenceIn-firstDynamic_];
               const double * solution = model->solutionRegion();
               setDynamicStatus(inBig, inSmall);
               if (lowerSet_[iSet] > -1.0e20)
                    lower_[iSet] += solution[sequenceIn];
               if (upperSet_[iSet] < 1.0e20)
                    upper_[iSet] += solution[sequenceIn];
               model->setObjectiveOffset(model->objectiveOffset() -
                                         solution[sequenceIn]*cost_[inBig]);
          }
     }
     }
     return returnNumber;
}
// Add a new variable to a set
void
ClpGubDynamicMatrix::insertNonBasic(int sequence, int iSet)
{
     int last = keyVariable_[iSet];
     int j = next_[last];
     while (j >= 0) {
          last = j;
          j = next_[j];
     }
     next_[last] = -(sequence + 1);
     next_[sequence] = j;
}
// Sets up an effective RHS and does gub crash if needed
void
ClpGubDynamicMatrix::useEffectiveRhs(ClpSimplex * model, bool cheapest)
{
     // Do basis - cheapest or slack if feasible (unless cheapest set)
     int longestSet = 0;
     int iSet;
     for (iSet = 0; iSet < numberSets_; iSet++)
          longestSet = CoinMax(longestSet, fullStart_[iSet+1] - fullStart_[iSet]);

     double * upper = new double[longestSet+1];
     double * cost = new double[longestSet+1];
     double * lower = new double[longestSet+1];
     double * solution = new double[longestSet+1];
     assert (!next_);
     delete [] next_;
     int numberColumns = model->numberColumns();
     next_ = new int[numberColumns+numberSets_+CoinMax(2*longestSet, lastDynamic_-firstDynamic_)];
     char * mark = new char[numberColumns];
     memset(mark, 0, numberColumns);
     for (int iColumn = 0; iColumn < numberColumns; iColumn++)
          next_[iColumn] = COIN_INT_MAX;
     int i;
     int * keys = new int[numberSets_];
     int * back = new int[numberGubColumns_];
     CoinFillN(back, numberGubColumns_, -1);
     for (i = 0; i < numberSets_; i++)
          keys[i] = COIN_INT_MAX;
     delete [] dynamicStatus_;
     dynamicStatus_ = new unsigned char [numberGubColumns_];
     memset(dynamicStatus_, 0, numberGubColumns_); // for clarity
     for (i = 0; i < numberGubColumns_; i++)
          setDynamicStatus(i, atLowerBound);
     // set up chains
     for (i = firstDynamic_; i < lastDynamic_; i++) {
          if (id_[i-firstDynamic_] >= 0) {
               if (model->getStatus(i) == ClpSimplex::basic)
                    mark[i] = 1;
               int iSet = backward_[i];
               assert (iSet >= 0);
               int iNext = keys[iSet];
               next_[i] = iNext;
               keys[iSet] = i;
               back[id_[i-firstDynamic_]] = i;
          } else {
               model->setStatus(i, ClpSimplex::atLowerBound);
               backward_[i] = -1;
          }
     }
     double * columnSolution = model->solutionRegion();
     int numberRows = getNumRows();
     toIndex_ = new int[numberSets_];
     for (iSet = 0; iSet < numberSets_; iSet++)
          toIndex_[iSet] = -1;
     fromIndex_ = new int [numberRows+numberSets_];
     double tolerance = model->primalTolerance();
     double * element =  matrix_->getMutableElements();
     int * row = matrix_->getMutableIndices();
     CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
     int * length = matrix_->getMutableVectorLengths();
     double objectiveOffset = 0.0;
     for (iSet = 0; iSet < numberSets_; iSet++) {
          int j;
          int numberBasic = 0;
          int iBasic = -1;
          int iStart = fullStart_[iSet];
          int iEnd = fullStart_[iSet+1];
          // find one with smallest length
          int smallest = numberRows + 1;
          double value = 0.0;
          j = keys[iSet];
          while (j != COIN_INT_MAX) {
               if (model->getStatus(j) == ClpSimplex::basic) {
                    if (length[j] < smallest) {
                         smallest = length[j];
                         iBasic = j;
                    }
                    numberBasic++;
               }
               value += columnSolution[j];
               j = next_[j];
          }
          bool done = false;
          if (numberBasic > 1 || (numberBasic == 1 && getStatus(iSet) == ClpSimplex::basic)) {
               if (getStatus(iSet) == ClpSimplex::basic)
                    iBasic = iSet + numberColumns; // slack key - use
               done = true;
          } else if (numberBasic == 1) {
               // see if can be key
               double thisSolution = columnSolution[iBasic];
               if (thisSolution < 0.0) {
                    value -= thisSolution;
                    thisSolution = 0.0;
                    columnSolution[iBasic] = thisSolution;
               }
               // try setting slack to a bound
               assert (upper_[iSet] < 1.0e20 || lower_[iSet] > -1.0e20);
               double cost1 = COIN_DBL_MAX;
               int whichBound = -1;
               if (upper_[iSet] < 1.0e20) {
                    // try slack at ub
                    double newBasic = thisSolution + upper_[iSet] - value;
                    if (newBasic >= -tolerance) {
                         // can go
                         whichBound = 1;
                         cost1 = newBasic * cost_[iBasic];
                         // But if exact then may be good solution
                         if (fabs(upper_[iSet] - value) < tolerance)
                              cost1 = -COIN_DBL_MAX;
                    }
               }
               if (lower_[iSet] > -1.0e20) {
                    // try slack at lb
                    double newBasic = thisSolution + lower_[iSet] - value;
                    if (newBasic >= -tolerance) {
                         // can go but is it cheaper
                         double cost2 = newBasic * cost_[iBasic];
                         // But if exact then may be good solution
                         if (fabs(lower_[iSet] - value) < tolerance)
                              cost2 = -COIN_DBL_MAX;
                         if (cost2 < cost1)
                              whichBound = 0;
                    }
               }
               if (whichBound != -1) {
                    // key
                    done = true;
                    if (whichBound) {
                         // slack to upper
                         columnSolution[iBasic] = thisSolution + upper_[iSet] - value;
                         setStatus(iSet, ClpSimplex::atUpperBound);
                    } else {
                         // slack to lower
                         columnSolution[iBasic] = thisSolution + lower_[iSet] - value;
                         setStatus(iSet, ClpSimplex::atLowerBound);
                    }
               }
          }
          if (!done) {
               if (!cheapest) {
                    // see if slack can be key
                    if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
                         done = true;
                         setStatus(iSet, ClpSimplex::basic);
                         iBasic = iSet + numberColumns;
                    }
               }
               if (!done) {
                    // set non basic if there was one
                    if (iBasic >= 0)
                         model->setStatus(iBasic, ClpSimplex::atLowerBound);
                    // find cheapest
                    int numberInSet = iEnd - iStart;
                    if (!lowerColumn_) {
                         CoinZeroN(lower, numberInSet);
                    } else {
                         for (int j = 0; j < numberInSet; j++)
                              lower[j] = lowerColumn_[j+iStart];
                    }
                    if (!upperColumn_) {
                         CoinFillN(upper, numberInSet, COIN_DBL_MAX);
                    } else {
                         for (int j = 0; j < numberInSet; j++)
                              upper[j] = upperColumn_[j+iStart];
                    }
                    CoinFillN(solution, numberInSet, 0.0);
                    // and slack
                    iBasic = numberInSet;
                    solution[iBasic] = -value;
                    lower[iBasic] = -upper_[iSet];
                    upper[iBasic] = -lower_[iSet];
                    int kphase;
                    if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
                         // feasible
                         kphase = 1;
                         cost[iBasic] = 0.0;
                         for (int j = 0; j < numberInSet; j++)
                              cost[j] = cost_[j+iStart];
                    } else {
                         // infeasible
                         kphase = 0;
                         // remember bounds are flipped so opposite to natural
                         if (value < lower_[iSet] - tolerance)
                              cost[iBasic] = 1.0;
                         else
                              cost[iBasic] = -1.0;
                         CoinZeroN(cost, numberInSet);
                    }
                    double dualTolerance = model->dualTolerance();
                    for (int iphase = kphase; iphase < 2; iphase++) {
                         if (iphase) {
                              cost[numberInSet] = 0.0;
                              for (int j = 0; j < numberInSet; j++)
                                   cost[j] = cost_[j+iStart];
                         }
                         // now do one row lp
                         bool improve = true;
                         while (improve) {
                              improve = false;
                              double dual = cost[iBasic];
                              int chosen = -1;
                              double best = dualTolerance;
                              int way = 0;
                              for (int i = 0; i <= numberInSet; i++) {
                                   double dj = cost[i] - dual;
                                   double improvement = 0.0;
                                   double distance = 0.0;
                                   if (iphase || i < numberInSet)
                                        assert (solution[i] >= lower[i] && solution[i] <= upper[i]);
                                   if (dj > dualTolerance)
                                        improvement = dj * (solution[i] - lower[i]);
                                   else if (dj < -dualTolerance)
                                        improvement = dj * (solution[i] - upper[i]);
                                   if (improvement > best) {
                                        best = improvement;
                                        chosen = i;
                                        if (dj < 0.0) {
                                             way = 1;
                                             distance = upper[i] - solution[i];
                                        } else {
                                             way = -1;
                                             distance = solution[i] - lower[i];
                                        }
                                   }
                              }
                              if (chosen >= 0) {
                                   improve = true;
                                   // now see how far
                                   if (way > 0) {
                                        // incoming increasing so basic decreasing
                                        // if phase 0 then go to nearest bound
                                        double distance = upper[chosen] - solution[chosen];
                                        double basicDistance;
                                        if (!iphase) {
                                             assert (iBasic == numberInSet);
                                             assert (solution[iBasic] > upper[iBasic]);
                                             basicDistance = solution[iBasic] - upper[iBasic];
                                        } else {
                                             basicDistance = solution[iBasic] - lower[iBasic];
                                        }
                                        // need extra coding for unbounded
                                        assert (CoinMin(distance, basicDistance) < 1.0e20);
                                        if (distance > basicDistance) {
                                             // incoming becomes basic
                                             solution[chosen] += basicDistance;
                                             if (!iphase)
                                                  solution[iBasic] = upper[iBasic];
                                             else
                                                  solution[iBasic] = lower[iBasic];
                                             iBasic = chosen;
                                        } else {
                                             // flip
                                             solution[chosen] = upper[chosen];
                                             solution[iBasic] -= distance;
                                        }
                                   } else {
                                        // incoming decreasing so basic increasing
                                        // if phase 0 then go to nearest bound
                                        double distance = solution[chosen] - lower[chosen];
                                        double basicDistance;
                                        if (!iphase) {
                                             assert (iBasic == numberInSet);
                                             assert (solution[iBasic] < lower[iBasic]);
                                             basicDistance = lower[iBasic] - solution[iBasic];
                                        } else {
                                             basicDistance = upper[iBasic] - solution[iBasic];
                                        }
                                        // need extra coding for unbounded - for now just exit
                                        if (CoinMin(distance, basicDistance) > 1.0e20) {
                                             printf("unbounded on set %d\n", iSet);
                                             iphase = 1;
                                             iBasic = numberInSet;
                                             break;
                                        }
                                        if (distance > basicDistance) {
                                             // incoming becomes basic
                                             solution[chosen] -= basicDistance;
                                             if (!iphase)
                                                  solution[iBasic] = lower[iBasic];
                                             else
                                                  solution[iBasic] = upper[iBasic];
                                             iBasic = chosen;
                                        } else {
                                             // flip
                                             solution[chosen] = lower[chosen];
                                             solution[iBasic] += distance;
                                        }
                                   }
                                   if (!iphase) {
                                        if(iBasic < numberInSet)
                                             break; // feasible
                                        else if (solution[iBasic] >= lower[iBasic] &&
                                                  solution[iBasic] <= upper[iBasic])
                                             break; // feasible (on flip)
                                   }
                              }
                         }
                    }
                    // do solution i.e. bounds
                    if (lowerColumn_ || upperColumn_) {
                         for (int j = 0; j < numberInSet; j++) {
                              if (j != iBasic) {
                                   objectiveOffset += solution[j] * cost[j];
                                   if (lowerColumn_ && upperColumn_) {
                                        if (fabs(solution[j] - lowerColumn_[j+iStart]) >
                                                  fabs(solution[j] - upperColumn_[j+iStart]))
                                             setDynamicStatus(j + iStart, atUpperBound);
                                   } else if (upperColumn_ && solution[j] > 0.0) {
                                        setDynamicStatus(j + iStart, atUpperBound);
                                   } else {
                                        setDynamicStatus(j + iStart, atLowerBound);
                                   }
                              }
                         }
                    }
                    // convert iBasic back and do bounds
                    if (iBasic == numberInSet) {
                         // slack basic
                         setStatus(iSet, ClpSimplex::basic);
                         iBasic = iSet + numberColumns;
                    } else {
                         iBasic += fullStart_[iSet];
                         if (back[iBasic] >= 0) {
                              // exists
                              iBasic = back[iBasic];
                         } else {
                              // create
                              CoinBigIndex numberElements = startColumn[firstAvailable_];
                              int numberThis = startColumn_[iBasic+1] - startColumn_[iBasic];
                              if (numberElements + numberThis > numberElements_) {
                                   // need to redo
                                   numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
                                   matrix_->reserve(numberColumns, numberElements_);
                                   element =  matrix_->getMutableElements();
                                   row = matrix_->getMutableIndices();
                                   // these probably okay but be safe
                                   startColumn = matrix_->getMutableVectorStarts();
                                   length = matrix_->getMutableVectorLengths();
                              }
                              length[firstAvailable_] = numberThis;
                              model->costRegion()[firstAvailable_] = cost_[iBasic];
                              if (lowerColumn_)
                                   model->lowerRegion()[firstAvailable_] = lowerColumn_[iBasic];
                              else
                                   model->lowerRegion()[firstAvailable_] = 0.0;
                              if (upperColumn_)
                                   model->upperRegion()[firstAvailable_] = upperColumn_[iBasic];
                              else
                                   model->upperRegion()[firstAvailable_] = COIN_DBL_MAX;
                              columnSolution[firstAvailable_] = solution[iBasic-fullStart_[iSet]];
                              CoinBigIndex base = startColumn_[iBasic];
                              for (int j = 0; j < numberThis; j++) {
                                   row[numberElements] = row_[base+j];
                                   element[numberElements++] = element_[base+j];
                              }
                              // already set startColumn[firstAvailable_]=numberElements;
                              id_[firstAvailable_-firstDynamic_] = iBasic;
                              setDynamicStatus(iBasic, inSmall);
                              backward_[firstAvailable_] = iSet;
                              iBasic = firstAvailable_;
                              firstAvailable_++;
                              startColumn[firstAvailable_] = numberElements;
                         }
                         model->setStatus(iBasic, ClpSimplex::basic);
                         // remember bounds flipped
                         if (upper[numberInSet] == lower[numberInSet])
                              setStatus(iSet, ClpSimplex::isFixed);
                         else if (solution[numberInSet] == upper[numberInSet])
                              setStatus(iSet, ClpSimplex::atLowerBound);
                         else if (solution[numberInSet] == lower[numberInSet])
                              setStatus(iSet, ClpSimplex::atUpperBound);
                         else
                              abort();
                    }
                    for (j = iStart; j < iEnd; j++) {
                         int iBack = back[j];
                         if (iBack >= 0) {
                              if (model->getStatus(iBack) != ClpSimplex::basic) {
                                   int inSet = j - iStart;
                                   columnSolution[iBack] = solution[inSet];
                                   if (upper[inSet] == lower[inSet])
                                        model->setStatus(iBack, ClpSimplex::isFixed);
                                   else if (solution[inSet] == upper[inSet])
                                        model->setStatus(iBack, ClpSimplex::atUpperBound);
                                   else if (solution[inSet] == lower[inSet])
                                        model->setStatus(iBack, ClpSimplex::atLowerBound);
                              }
                         }
                    }
               }
          }
          keyVariable_[iSet] = iBasic;
     }
     model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
     delete [] lower;
     delete [] solution;
     delete [] upper;
     delete [] cost;
     // make sure matrix is in good shape
     matrix_->orderMatrix();
     // create effective rhs
     delete [] rhsOffset_;
     rhsOffset_ = new double[numberRows];
     // and redo chains
     memset(mark, 0, numberColumns);
     for (int iColumnX = 0; iColumnX < firstAvailable_; iColumnX++)
          next_[iColumnX] = COIN_INT_MAX;
     for (i = 0; i < numberSets_; i++) {
          keys[i] = COIN_INT_MAX;
          int iKey = keyVariable_[i];
          if (iKey < numberColumns)
               model->setStatus(iKey, ClpSimplex::basic);
     }
     // set up chains
     for (i = 0; i < firstAvailable_; i++) {
          if (model->getStatus(i) == ClpSimplex::basic)
               mark[i] = 1;
          int iSet = backward_[i];
          if (iSet >= 0) {
               int iNext = keys[iSet];
               next_[i] = iNext;
               keys[iSet] = i;
          }
     }
     for (i = 0; i < numberSets_; i++) {
          if (keys[i] != COIN_INT_MAX) {
               // something in set
               int j;
               if (getStatus(i) != ClpSimplex::basic) {
                    // make sure fixed if it is
                    if (upper_[i] == lower_[i])
                         setStatus(i, ClpSimplex::isFixed);
                    // slack not key - choose one with smallest length
                    int smallest = numberRows + 1;
                    int key = -1;
                    j = keys[i];
                    while (1) {
                         if (mark[j] && length[j] < smallest) {
                              key = j;
                              smallest = length[j];
                         }
                         if (next_[j] != COIN_INT_MAX) {
                              j = next_[j];
                         } else {
                              // correct end
                              next_[j] = -(keys[i] + 1);
                              break;
                         }
                    }
                    if (key >= 0) {
                         keyVariable_[i] = key;
                    } else {
                         // nothing basic - make slack key
                         //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
                         // fudge to avoid const problem
                         status_[i] = 1;
                    }
               } else {
                    // slack key
                    keyVariable_[i] = numberColumns + i;
                    int j;
                    double sol = 0.0;
                    j = keys[i];
                    while (1) {
                         sol += columnSolution[j];
                         if (next_[j] != COIN_INT_MAX) {
                              j = next_[j];
                         } else {
                              // correct end
                              next_[j] = -(keys[i] + 1);
                              break;
                         }
                    }
                    if (sol > upper_[i] + tolerance) {
                         setAbove(i);
                    } else if (sol < lower_[i] - tolerance) {
                         setBelow(i);
                    } else {
                         setFeasible(i);
                    }
               }
               // Create next_
               int key = keyVariable_[i];
               redoSet(model, key, keys[i], i);
          } else {
               // nothing in set!
               next_[i+numberColumns] = -(i + numberColumns + 1);
               keyVariable_[i] = numberColumns + i;
               double sol = 0.0;
               if (sol > upper_[i] + tolerance) {
                    setAbove(i);
               } else if (sol < lower_[i] - tolerance) {
                    setBelow(i);
               } else {
                    setFeasible(i);
               }
          }
     }
     delete [] keys;
     delete [] mark;
     delete [] back;
     rhsOffset(model, true);
}
/* Returns effective RHS if it is being used.  This is used for long problems
   or big gub or anywhere where going through full columns is
   expensive.  This may re-compute */
double *
ClpGubDynamicMatrix::rhsOffset(ClpSimplex * model, bool forceRefresh,
                               bool
#ifdef CLP_DEBUG
                               check
#endif
                              )
{
     //forceRefresh=true;
     //check=false;
#ifdef CLP_DEBUG
     double * saveE = NULL;
     if (rhsOffset_ && check) {
          int numberRows = model->numberRows();
          saveE = new double[numberRows];
     }
#endif
     if (rhsOffset_) {
#ifdef CLP_DEBUG
          if (check) {
               // no need - but check anyway
               int numberRows = model->numberRows();
               double * rhs = new double[numberRows];
               int numberColumns = model->numberColumns();
               int iRow;
               CoinZeroN(rhs, numberRows);
               // do ones at bounds before gub
               const double * smallSolution = model->solutionRegion();
               const double * element = matrix_->getElements();
               const int * row = matrix_->getIndices();
               const CoinBigIndex * startColumn = matrix_->getVectorStarts();
               const int * length = matrix_->getVectorLengths();
               int iColumn;
               for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
                    if (model->getStatus(iColumn) != ClpSimplex::basic) {
                         double value = smallSolution[iColumn];
                         for (CoinBigIndex j = startColumn[iColumn];
                                   j < startColumn[iColumn] + length[iColumn]; j++) {
                              int jRow = row[j];
                              rhs[jRow] -= value * element[j];
                         }
                    }
               }
               if (lowerColumn_ || upperColumn_) {
                    double * solution = new double [numberGubColumns_];
                    for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
                         double value = 0.0;
                         if(getDynamicStatus(iColumn) == atUpperBound)
                              value = upperColumn_[iColumn];
                         else if (lowerColumn_)
                              value = lowerColumn_[iColumn];
                         solution[iColumn] = value;
                    }
                    // ones at bounds in small and gub
                    for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
                         int jFull = id_[iColumn-firstDynamic_];
                         solution[jFull] = smallSolution[iColumn];
                    }
                    // zero all basic in small model
                    int * pivotVariable = model->pivotVariable();
                    for (iRow = 0; iRow < numberRows; iRow++) {
                         int iColumn = pivotVariable[iRow];
                         if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
                              int iSequence = id_[iColumn-firstDynamic_];
                              solution[iSequence] = 0.0;
                         }
                    }
                    // and now compute value to use for key
                    ClpSimplex::Status iStatus;
                    for (int iSet = 0; iSet < numberSets_; iSet++) {
                         iColumn = keyVariable_[iSet];
                         if (iColumn < numberColumns) {
                              int iSequence = id_[iColumn-firstDynamic_];
                              solution[iSequence] = 0.0;
                              double b = 0.0;
                              // key is structural - where is slack
                              iStatus = getStatus(iSet);
                              assert (iStatus != ClpSimplex::basic);
                              if (iStatus == ClpSimplex::atLowerBound)
                                   b = lowerSet_[iSet];
                              else
                                   b = upperSet_[iSet];
                              // subtract out others at bounds
                              for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
                                   b -= solution[j];
                              solution[iSequence] = b;
                         }
                    }
                    for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
                         double value = solution[iColumn];
                         if (value) {
                              for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
                                   int iRow = row_[j];
                                   rhs[iRow] -= element_[j] * value;
                              }
                         }
                    }
                    // now do lower and upper bounds on sets
                    for (int iSet = 0; iSet < numberSets_; iSet++) {
                         iColumn = keyVariable_[iSet];
                         double shift = 0.0;
                         for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++) {
                              if (getDynamicStatus(j) != inSmall && j != iColumn) {
                                   if (getDynamicStatus(j) == atLowerBound) {
                                        if (lowerColumn_)
                                             shift += lowerColumn_[j];
                                   } else {
                                        shift += upperColumn_[j];
                                   }
                              }
                         }
                         if (lowerSet_[iSet] > -1.0e20)
                              assert(fabs(lower_[iSet] - (lowerSet_[iSet] - shift)) < 1.0e-3);
                         if (upperSet_[iSet] < 1.0e20)
                              assert(fabs(upper_[iSet] - ( upperSet_[iSet] - shift)) < 1.0e-3);
                    }
                    delete [] solution;
               } else {
                    // no bounds
                    ClpSimplex::Status iStatus;
                    for (int iSet = 0; iSet < numberSets_; iSet++) {
                         int iColumn = keyVariable_[iSet];
                         if (iColumn < numberColumns) {
                              int iSequence = id_[iColumn-firstDynamic_];
                              double b = 0.0;
                              // key is structural - where is slack
                              iStatus = getStatus(iSet);
                              assert (iStatus != ClpSimplex::basic);
                              if (iStatus == ClpSimplex::atLowerBound)
                                   b = lower_[iSet];
                              else
                                   b = upper_[iSet];
                              if (b) {
                                   for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence+1]; j++) {
                                        int iRow = row_[j];
                                        rhs[iRow] -= element_[j] * b;
                                   }
                              }
                         }
                    }
               }
               for (iRow = 0; iRow < numberRows; iRow++) {
                    if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
                         printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
               }
               CoinMemcpyN(rhs, numberRows, saveE);
               delete [] rhs;
          }
#endif
          if (forceRefresh || (refreshFrequency_ && model->numberIterations() >=
                               lastRefresh_ + refreshFrequency_)) {
               int numberRows = model->numberRows();
               int numberColumns = model->numberColumns();
               int iRow;
               CoinZeroN(rhsOffset_, numberRows);
               // do ones at bounds before gub
               const double * smallSolution = model->solutionRegion();
               const double * element = matrix_->getElements();
               const int * row = matrix_->getIndices();
               const CoinBigIndex * startColumn = matrix_->getVectorStarts();
               const int * length = matrix_->getVectorLengths();
               int iColumn;
               for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
                    if (model->getStatus(iColumn) != ClpSimplex::basic) {
                         double value = smallSolution[iColumn];
                         for (CoinBigIndex j = startColumn[iColumn];
                                   j < startColumn[iColumn] + length[iColumn]; j++) {
                              int jRow = row[j];
                              rhsOffset_[jRow] -= value * element[j];
                         }
                    }
               }
               if (lowerColumn_ || upperColumn_) {
                    double * solution = new double [numberGubColumns_];
                    for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
                         double value = 0.0;
                         if(getDynamicStatus(iColumn) == atUpperBound)
                              value = upperColumn_[iColumn];
                         else if (lowerColumn_)
                              value = lowerColumn_[iColumn];
                         solution[iColumn] = value;
                    }
                    // ones in gub and in small problem
                    for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
                         int jFull = id_[iColumn-firstDynamic_];
                         solution[jFull] = smallSolution[iColumn];
                    }
                    // zero all basic in small model
                    int * pivotVariable = model->pivotVariable();
                    for (iRow = 0; iRow < numberRows; iRow++) {
                         int iColumn = pivotVariable[iRow];
                         if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
                              int iSequence = id_[iColumn-firstDynamic_];
                              solution[iSequence] = 0.0;
                         }
                    }
                    // and now compute value to use for key
                    ClpSimplex::Status iStatus;
                    int iSet;
                    for ( iSet = 0; iSet < numberSets_; iSet++) {
                         iColumn = keyVariable_[iSet];
                         if (iColumn < numberColumns) {
                              int iSequence = id_[iColumn-firstDynamic_];
                              solution[iSequence] = 0.0;
                              double b = 0.0;
                              // key is structural - where is slack
                              iStatus = getStatus(iSet);
                              assert (iStatus != ClpSimplex::basic);
                              if (iStatus == ClpSimplex::atLowerBound)
                                   b = lowerSet_[iSet];
                              else
                                   b = upperSet_[iSet];
                              // subtract out others at bounds
                              for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
                                   b -= solution[j];
                              solution[iSequence] = b;
                         }
                    }
                    for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
                         double value = solution[iColumn];
                         if (value) {
                              for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
                                   int iRow = row_[j];
                                   rhsOffset_[iRow] -= element_[j] * value;
                              }
                         }
                    }
                    // now do lower and upper bounds on sets
                    // and offset
                    double objectiveOffset = 0.0;
                    for ( iSet = 0; iSet < numberSets_; iSet++) {
                         iColumn = keyVariable_[iSet];
                         double shift = 0.0;
                         for (CoinBigIndex j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++) {
                              if (getDynamicStatus(j) != inSmall) {
                                   double value = 0.0;
                                   if (getDynamicStatus(j) == atLowerBound) {
                                        if (lowerColumn_)
                                             value = lowerColumn_[j];
                                   } else {
                                        value = upperColumn_[j];
                                   }
                                   if (j != iColumn)
                                        shift += value;
                                   objectiveOffset += value * cost_[j];
                              }
                         }
                         if (lowerSet_[iSet] > -1.0e20)
                              lower_[iSet] = lowerSet_[iSet] - shift;
                         if (upperSet_[iSet] < 1.0e20)
                              upper_[iSet] = upperSet_[iSet] - shift;
                    }
                    delete [] solution;
                    model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
               } else {
                    // no bounds
                    ClpSimplex::Status iStatus;
                    for (int iSet = 0; iSet < numberSets_; iSet++) {
                         int iColumn = keyVariable_[iSet];
                         if (iColumn < numberColumns) {
                              int iSequence = id_[iColumn-firstDynamic_];
                              double b = 0.0;
                              // key is structural - where is slack
                              iStatus = getStatus(iSet);
                              assert (iStatus != ClpSimplex::basic);
                              if (iStatus == ClpSimplex::atLowerBound)
                                   b = lower_[iSet];
                              else
                                   b = upper_[iSet];
                              if (b) {
                                   for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence+1]; j++) {
                                        int iRow = row_[j];
                                        rhsOffset_[iRow] -= element_[j] * b;
                                   }
                              }
                         }
                    }
               }
#ifdef CLP_DEBUG
               if (saveE) {
                    for (iRow = 0; iRow < numberRows; iRow++) {
                         if (fabs(saveE[iRow] - rhsOffset_[iRow]) > 1.0e-3)
                              printf("** %d - old eff %g new %g\n", iRow, saveE[iRow], rhsOffset_[iRow]);
                    }
                    delete [] saveE;
               }
#endif
               lastRefresh_ = model->numberIterations();
          }
     }
     return rhsOffset_;
}
/*
  update information for a pivot (and effective rhs)
*/
int
ClpGubDynamicMatrix::updatePivot(ClpSimplex * model, double oldInValue, double oldOutValue)
{

     // now update working model
     int sequenceIn = model->sequenceIn();
     int sequenceOut = model->sequenceOut();
     bool doPrinting = (model->messageHandler()->logLevel() == 63);
     bool print = false;
     int iSet;
     int trueIn = -1;
     int trueOut = -1;
     int numberRows = model->numberRows();
     int numberColumns = model->numberColumns();
     if (sequenceIn == firstAvailable_) {
          if (doPrinting)
               printf("New variable ");
          if (sequenceIn != sequenceOut) {
               insertNonBasic(firstAvailable_, backward_[firstAvailable_]);
               setDynamicStatus(id_[sequenceIn-firstDynamic_], inSmall);
               firstAvailable_++;
          } else {
               int bigSequence = id_[sequenceIn-firstDynamic_];
               if (model->getStatus(sequenceIn) == ClpSimplex::atUpperBound)
                    setDynamicStatus(bigSequence, atUpperBound);
               else
                    setDynamicStatus(bigSequence, atLowerBound);
          }
          synchronize(model, 8);
     }
     if (sequenceIn < lastDynamic_) {
          iSet = backward_[sequenceIn];
          if (iSet >= 0) {
               int bigSequence = id_[sequenceIn-firstDynamic_];
               trueIn = bigSequence + numberRows + numberColumns + numberSets_;
               if (doPrinting)
                    printf(" incoming set %d big seq %d", iSet, bigSequence);
               print = true;
          }
     } else if (sequenceIn >= numberRows + numberColumns) {
          trueIn = numberRows + numberColumns + gubSlackIn_;
     }
     if (sequenceOut < lastDynamic_) {
          iSet = backward_[sequenceOut];
          if (iSet >= 0) {
               int bigSequence = id_[sequenceOut-firstDynamic_];
               trueOut = bigSequence + firstDynamic_;
               if (getDynamicStatus(bigSequence) != inSmall) {
                    if (model->getStatus(sequenceOut) == ClpSimplex::atUpperBound)
                         setDynamicStatus(bigSequence, atUpperBound);
                    else
                         setDynamicStatus(bigSequence, atLowerBound);
               }
               if (doPrinting)
                    printf(" ,outgoing set %d big seq %d,", iSet, bigSequence);
               print = true;
               model->setSequenceIn(sequenceOut);
               synchronize(model, 8);
               model->setSequenceIn(sequenceIn);
          }
     }
     if (print && doPrinting)
          printf("\n");
     ClpGubMatrix::updatePivot(model, oldInValue, oldOutValue);
     // Redo true in and out
     if (trueIn >= 0)
          trueSequenceIn_ = trueIn;
     if (trueOut >= 0)
          trueSequenceOut_ = trueOut;
     if (doPrinting && 0) {
          for (int i = 0; i < numberSets_; i++) {
               printf("set %d key %d lower %g upper %g\n", i, keyVariable_[i], lower_[i], upper_[i]);
               for (int j = fullStart_[i]; j < fullStart_[i+1]; j++)
                    if (getDynamicStatus(j) == atUpperBound) {
                         bool print = true;
                         for (int k = firstDynamic_; k < firstAvailable_; k++) {
                              if (id_[k-firstDynamic_] == j)
                                   print = false;
                              if (id_[k-firstDynamic_] == j)
                                   assert(getDynamicStatus(j) == inSmall);
                         }
                         if (print)
                              printf("variable %d at ub\n", j);
                    }
          }
     }
#ifdef CLP_DEBUG
     char * inSmall = new char [numberGubColumns_];
     memset(inSmall, 0, numberGubColumns_);
     for (int i = 0; i < numberGubColumns_; i++)
          if (getDynamicStatus(i) == ClpGubDynamicMatrix::inSmall)
               inSmall[i] = 1;
     for (int i = firstDynamic_; i < firstAvailable_; i++) {
          int k = id_[i-firstDynamic_];
          inSmall[k] = 0;
     }
     for (int i = 0; i < numberGubColumns_; i++)
          assert (!inSmall[i]);
     delete [] inSmall;
#endif
     return 0;
}
void
ClpGubDynamicMatrix::times(double scalar,
                           const double * x, double * y) const
{
     if (model_->specialOptions() != 16) {
          ClpPackedMatrix::times(scalar, x, y);
     } else {
          int iRow;
          int numberColumns = model_->numberColumns();
          int numberRows = model_->numberRows();
          const double * element =  matrix_->getElements();
          const int * row = matrix_->getIndices();
          const CoinBigIndex * startColumn = matrix_->getVectorStarts();
          const int * length = matrix_->getVectorLengths();
          int * pivotVariable = model_->pivotVariable();
          int numberToDo = 0;
          for (iRow = 0; iRow < numberRows; iRow++) {
               y[iRow] -= scalar * rhsOffset_[iRow];
               int iColumn = pivotVariable[iRow];
               if (iColumn < numberColumns) {
                    int iSet = backward_[iColumn];
                    if (iSet >= 0 && toIndex_[iSet] < 0) {
                         toIndex_[iSet] = 0;
                         fromIndex_[numberToDo++] = iSet;
                    }
                    CoinBigIndex j;
                    double value = scalar * x[iColumn];
                    if (value) {
                         for (j = startColumn[iColumn];
                                   j < startColumn[iColumn] + length[iColumn]; j++) {
                              int jRow = row[j];
                              y[jRow] += value * element[j];
                         }
                    }
               }
          }
          // and gubs which are interacting
          for (int jSet = 0; jSet < numberToDo; jSet++) {
               int iSet = fromIndex_[jSet];
               toIndex_[iSet] = -1;
               int iKey = keyVariable_[iSet];
               if (iKey < numberColumns) {
                    double valueKey;
                    if (getStatus(iSet) == ClpSimplex::atLowerBound)
                         valueKey = lower_[iSet];
                    else
                         valueKey = upper_[iSet];
                    double value = scalar * (x[iKey] - valueKey);
                    if (value) {
                         for (CoinBigIndex j = startColumn[iKey];
                                   j < startColumn[iKey] + length[iKey]; j++) {
                              int jRow = row[j];
                              y[jRow] += value * element[j];
                         }
                    }
               }
          }
     }
}
/* Just for debug - may be extended to other matrix types later.
   Returns number and sum of primal infeasibilities.
*/
int
ClpGubDynamicMatrix::checkFeasible(ClpSimplex * /*model*/, double & sum) const
{
     int numberRows = model_->numberRows();
     double * rhs = new double[numberRows];
     int numberColumns = model_->numberColumns();
     int iRow;
     CoinZeroN(rhs, numberRows);
     // do ones at bounds before gub
     const double * smallSolution = model_->solutionRegion();
     const double * element = matrix_->getElements();
     const int * row = matrix_->getIndices();
     const CoinBigIndex * startColumn = matrix_->getVectorStarts();
     const int * length = matrix_->getVectorLengths();
     int iColumn;
     int numberInfeasible = 0;
     const double * rowLower = model_->rowLower();
     const double * rowUpper = model_->rowUpper();
     sum = 0.0;
     for (iRow = 0; iRow < numberRows; iRow++) {
          double value = smallSolution[numberColumns+iRow];
          if (value < rowLower[iRow] - 1.0e-5 ||
                    value > rowUpper[iRow] + 1.0e-5) {
               //printf("row %d %g %g %g\n",
               //     iRow,rowLower[iRow],value,rowUpper[iRow]);
               numberInfeasible++;
               sum += CoinMax(rowLower[iRow] - value, value - rowUpper[iRow]);
          }
          rhs[iRow] = value;
     }
     const double * columnLower = model_->columnLower();
     const double * columnUpper = model_->columnUpper();
     for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
          double value = smallSolution[iColumn];
          if (value < columnLower[iColumn] - 1.0e-5 ||
                    value > columnUpper[iColumn] + 1.0e-5) {
               //printf("column %d %g %g %g\n",
               //     iColumn,columnLower[iColumn],value,columnUpper[iColumn]);
               numberInfeasible++;
               sum += CoinMax(columnLower[iColumn] - value, value - columnUpper[iColumn]);
          }
          for (CoinBigIndex j = startColumn[iColumn];
                    j < startColumn[iColumn] + length[iColumn]; j++) {
               int jRow = row[j];
               rhs[jRow] -= value * element[j];
          }
     }
     double * solution = new double [numberGubColumns_];
     for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
          double value = 0.0;
          if(getDynamicStatus(iColumn) == atUpperBound)
               value = upperColumn_[iColumn];
          else if (lowerColumn_)
               value = lowerColumn_[iColumn];
          solution[iColumn] = value;
     }
     // ones in small and gub
     for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
          int jFull = id_[iColumn-firstDynamic_];
          solution[jFull] = smallSolution[iColumn];
     }
     // fill in all basic in small model
     int * pivotVariable = model_->pivotVariable();
     for (iRow = 0; iRow < numberRows; iRow++) {
          int iColumn = pivotVariable[iRow];
          if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
               int iSequence = id_[iColumn-firstDynamic_];
               solution[iSequence] = smallSolution[iColumn];
          }
     }
     // and now compute value to use for key
     ClpSimplex::Status iStatus;
     for (int iSet = 0; iSet < numberSets_; iSet++) {
          iColumn = keyVariable_[iSet];
          if (iColumn < numberColumns) {
               int iSequence = id_[iColumn-firstDynamic_];
               solution[iSequence] = 0.0;
               double b = 0.0;
               // key is structural - where is slack
               iStatus = getStatus(iSet);
               assert (iStatus != ClpSimplex::basic);
               if (iStatus == ClpSimplex::atLowerBound)
                    b = lower_[iSet];
               else
                    b = upper_[iSet];
               // subtract out others at bounds
               for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
                    b -= solution[j];
               solution[iSequence] = b;
          }
     }
     for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
          double value = solution[iColumn];
          if ((lowerColumn_ && value < lowerColumn_[iColumn] - 1.0e-5) ||
                    (!lowerColumn_ && value < -1.0e-5) ||
                    (upperColumn_ && value > upperColumn_[iColumn] + 1.0e-5)) {
               //printf("column %d %g %g %g\n",
               //     iColumn,lowerColumn_[iColumn],value,upperColumn_[iColumn]);
               numberInfeasible++;
          }
          if (value) {
               for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
                    int iRow = row_[j];
                    rhs[iRow] -= element_[j] * value;
               }
          }
     }
     for (iRow = 0; iRow < numberRows; iRow++) {
          if (fabs(rhs[iRow]) > 1.0e-5)
               printf("rhs mismatch %d %g\n", iRow, rhs[iRow]);
     }
     delete [] solution;
     delete [] rhs;
     return numberInfeasible;
}
// Cleans data after setWarmStart
void
ClpGubDynamicMatrix::cleanData(ClpSimplex * model)
{
     // and redo chains
     int numberColumns = model->numberColumns();
     int iColumn;
     // do backward
     int * mark = new int [numberGubColumns_];
     for (iColumn = 0; iColumn < numberGubColumns_; iColumn++)
          mark[iColumn] = -1;
     int i;
     for (i = 0; i < firstDynamic_; i++) {
          assert (backward_[i] == -1);
          next_[i] = -1;
     }
     for (i = firstDynamic_; i < firstAvailable_; i++) {
          iColumn = id_[i-firstDynamic_];
          mark[iColumn] = i;
     }
     for (i = 0; i < numberSets_; i++) {
          int iKey = keyVariable_[i];
          int lastNext = -1;
          int firstNext = -1;
          for (CoinBigIndex k = fullStart_[i]; k < fullStart_[i+1]; k++) {
               iColumn = mark[k];
               if (iColumn >= 0) {
                    if (iColumn != iKey) {
                         if (lastNext >= 0)
                              next_[lastNext] = iColumn;
                         else
                              firstNext = iColumn;
                         lastNext = iColumn;
                    }
                    backward_[iColumn] = i;
               }
          }
          setFeasible(i);
          if (firstNext >= 0) {
               // others
               next_[iKey] = firstNext;
               next_[lastNext] = -(iKey + 1);
          } else if (iKey < numberColumns) {
               next_[iKey] = -(iKey + 1);
          }
     }
     delete [] mark;
     // fill matrix
     double * element =  matrix_->getMutableElements();
     int * row = matrix_->getMutableIndices();
     CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
     int * length = matrix_->getMutableVectorLengths();
     CoinBigIndex numberElements = startColumn[firstDynamic_];
     for (i = firstDynamic_; i < firstAvailable_; i++) {
          int iColumn = id_[i-firstDynamic_];
          int numberThis = startColumn_[iColumn+1] - startColumn_[iColumn];
          length[i] = numberThis;
          for (CoinBigIndex jBigIndex = startColumn_[iColumn];
                    jBigIndex < startColumn_[iColumn+1]; jBigIndex++) {
               row[numberElements] = row_[jBigIndex];
               element[numberElements++] = element_[jBigIndex];
          }
          startColumn[i+1] = numberElements;
     }
}
