/* $Id$ */
// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#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;
              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;
                } else {
                  way = -1;
                }
              }
            }
            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;
  }
}

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/