/* $Id: ClpDynamicMatrix.cpp 1941 2013-04-10 16:52:27Z stefan $ */
// Copyright (C) 2004, 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 "ClpDynamicMatrix.hpp"
#include "ClpMessage.hpp"
//#define CLP_DEBUG
//#define CLP_DEBUG_PRINT
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpDynamicMatrix::ClpDynamicMatrix ()
     : ClpPackedMatrix(),
       sumDualInfeasibilities_(0.0),
       sumPrimalInfeasibilities_(0.0),
       sumOfRelaxedDualInfeasibilities_(0.0),
       sumOfRelaxedPrimalInfeasibilities_(0.0),
       savedBestGubDual_(0.0),
       savedBestSet_(0),
       backToPivotRow_(NULL),
       keyVariable_(NULL),
       toIndex_(NULL),
       fromIndex_(NULL),
       numberSets_(0),
       numberActiveSets_(0),
       objectiveOffset_(0.0),
       lowerSet_(NULL),
       upperSet_(NULL),
       status_(NULL),
       model_(NULL),
       firstAvailable_(0),
       firstAvailableBefore_(0),
       firstDynamic_(0),
       lastDynamic_(0),
       numberStaticRows_(0),
       numberElements_(0),
       numberDualInfeasibilities_(0),
       numberPrimalInfeasibilities_(0),
       noCheck_(-1),
       infeasibilityWeight_(0.0),
       numberGubColumns_(0),
       maximumGubColumns_(0),
       maximumElements_(0),
       startSet_(NULL),
       next_(NULL),
       startColumn_(NULL),
       row_(NULL),
       element_(NULL),
       cost_(NULL),
       id_(NULL),
       dynamicStatus_(NULL),
       columnLower_(NULL),
       columnUpper_(NULL)
{
     setType(15);
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpDynamicMatrix::ClpDynamicMatrix (const ClpDynamicMatrix & rhs)
     : ClpPackedMatrix(rhs)
{
     objectiveOffset_ = rhs.objectiveOffset_;
     numberSets_ = rhs.numberSets_;
     numberActiveSets_ = rhs.numberActiveSets_;
     firstAvailable_ = rhs.firstAvailable_;
     firstAvailableBefore_ = rhs.firstAvailableBefore_;
     firstDynamic_ = rhs.firstDynamic_;
     lastDynamic_ = rhs.lastDynamic_;
     numberStaticRows_ = rhs.numberStaticRows_;
     numberElements_ = rhs.numberElements_;
     backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, lastDynamic_);
     keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
     toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
     fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + 1 - numberStaticRows_);
     lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
     upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
     status_ = ClpCopyOfArray(rhs.status_, static_cast<int>(2*numberSets_+4*sizeof(int)));
     model_ = rhs.model_;
     sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
     sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
     sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
     sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
     numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
     numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
     savedBestGubDual_ = rhs.savedBestGubDual_;
     savedBestSet_ = rhs.savedBestSet_;
     noCheck_ = rhs.noCheck_;
     infeasibilityWeight_ = rhs.infeasibilityWeight_;
     // Now secondary data
     numberGubColumns_ = rhs.numberGubColumns_;
     maximumGubColumns_ = rhs.maximumGubColumns_;
     maximumElements_ = rhs.maximumElements_;
     startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1);
     next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_);
     startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1);
     row_ = ClpCopyOfArray(rhs.row_, maximumElements_);
     element_ = ClpCopyOfArray(rhs.element_, maximumElements_);
     cost_ = ClpCopyOfArray(rhs.cost_, maximumGubColumns_);
     id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
     columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_);
     columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_);
     dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_);
}

/* This is the real constructor*/
ClpDynamicMatrix::ClpDynamicMatrix(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 * columnLower, const double * columnUpper,
                                   const unsigned char * status,
                                   const unsigned char * dynamicStatus)
     : ClpPackedMatrix()
{
     setType(15);
     objectiveOffset_ = model->objectiveOffset();
     model_ = model;
     numberSets_ = numberSets;
     numberGubColumns_ = numberGubColumns;
     maximumGubColumns_ = numberGubColumns_;
     if (numberGubColumns_)
          maximumElements_ = startColumn[numberGubColumns_];
     else
          maximumElements_ = 0;
     startSet_ = new int [numberSets_+1];
     next_ = new int [maximumGubColumns_];
     // fill in startSet and next
     int iSet;
     if (numberGubColumns_) {
          for (iSet = 0; iSet < numberSets_; iSet++) {
               int first = starts[iSet];
               int last = starts[iSet+1] - 1;
               startSet_[iSet] = first;
               for (int i = first; i < last; i++)
                    next_[i] = i + 1;
               next_[last] = -iSet - 1;
          }
	  startSet_[numberSets_] = starts[numberSets_];
     }
     int numberColumns = model->numberColumns();
     int numberRows = model->numberRows();
     numberStaticRows_ = numberRows;
     savedBestGubDual_ = 0.0;
     savedBestSet_ = 0;
     // Number of columns needed
     int frequency = model->factorizationFrequency();
     int numberGubInSmall = numberRows + frequency + CoinMin(frequency, numberSets_) + 4;
     // But we may have two per row + one for incoming (make it two) 
     numberGubInSmall = CoinMax(2*numberRows+2,numberGubInSmall);
     // for small problems this could be too big
     //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
     int numberNeeded = numberGubInSmall + numberColumns;
     firstAvailable_ = numberColumns;
     firstAvailableBefore_ = firstAvailable_;
     firstDynamic_ = numberColumns;
     lastDynamic_ = numberNeeded;
     startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1);
     if (!numberGubColumns_) {
       if (!startColumn_)
	 startColumn_ = new CoinBigIndex [1];
       startColumn_[0] = 0;
     }
     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];
          // I don't think I need sorted but ...
          CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i+1], element_ + startColumn_[i]);
     }
     if (columnLower) {
          columnLower_ = new double[numberGubColumns_];
          for (i = 0; i < numberGubColumns_; i++)
               columnLower_[i] = columnLower[i];
     } else {
          columnLower_ = NULL;
     }
     if (columnUpper) {
          columnUpper_ = new double[numberGubColumns_];
          for (i = 0; i < numberGubColumns_; i++)
               columnUpper_[i] = columnUpper[i];
     } else {
          columnUpper_ = NULL;
     }
     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;
     }
     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 = numberElements;
     guess /= static_cast<double> (numberColumns);
     guess *= 2 * numberGubInSmall;
     numberElements_ = static_cast<int> (guess);
     numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements();
     matrix_ = originalMatrix;
     //delete originalMatrixA;
     flags_ &= ~1;
     // resize model (matrix stays same)
     // modify frequency
     if (frequency>=50)
       frequency = 50+(frequency-50)/2;
     int newRowSize = numberRows + CoinMin(numberSets_, frequency+numberRows) + 1;
     model->resize(newRowSize, numberNeeded);
     for (i = numberRows; i < newRowSize; i++)
          model->setRowStatus(i, ClpSimplex::basic);
     if (columnUpper_) {
          // 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();
     originalMatrix->setDimensions(newRowSize, -1);
     numberActiveColumns_ = firstDynamic_;
     // redo number of columns
     numberColumns = matrix_->getNumCols();
     backToPivotRow_ = new int[numberNeeded];
     keyVariable_ = new int[numberSets_];
     if (status) {
          status_ = ClpCopyOfArray(status, static_cast<int>(2*numberSets_+4*sizeof(int)));
          assert (dynamicStatus);
          dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2*numberGubColumns_);
     } else {
          status_ = new unsigned char [2*numberSets_+4*sizeof(int)];
          memset(status_, 0, numberSets_);
          int i;
          for (i = 0; i < numberSets_; i++) {
               // make slack key
               setStatus(i, ClpSimplex::basic);
          }
          dynamicStatus_ = new unsigned char [2*numberGubColumns_];
          memset(dynamicStatus_, 0, numberGubColumns_); // for clarity
          for (i = 0; i < numberGubColumns_; i++)
               setDynamicStatus(i, atLowerBound);
     }
     toIndex_ = new int[numberSets_];
     for (iSet = 0; iSet < numberSets_; iSet++)
          toIndex_[iSet] = -1;
     fromIndex_ = new int [newRowSize-numberStaticRows_+1];
     numberActiveSets_ = 0;
     rhsOffset_ = NULL;
     if (numberGubColumns_) {
          if (!status) {
               gubCrash();
          } else {
               initialProblem();
          }
     }
     noCheck_ = -1;
     infeasibilityWeight_ = 0.0;
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpDynamicMatrix::~ClpDynamicMatrix ()
{
     delete [] backToPivotRow_;
     delete [] keyVariable_;
     delete [] toIndex_;
     delete [] fromIndex_;
     delete [] lowerSet_;
     delete [] upperSet_;
     delete [] status_;
     delete [] startSet_;
     delete [] next_;
     delete [] startColumn_;
     delete [] row_;
     delete [] element_;
     delete [] cost_;
     delete [] id_;
     delete [] dynamicStatus_;
     delete [] columnLower_;
     delete [] columnUpper_;
}

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpDynamicMatrix &
ClpDynamicMatrix::operator=(const ClpDynamicMatrix& rhs)
{
     if (this != &rhs) {
          ClpPackedMatrix::operator=(rhs);
          delete [] backToPivotRow_;
          delete [] keyVariable_;
          delete [] toIndex_;
          delete [] fromIndex_;
          delete [] lowerSet_;
          delete [] upperSet_;
          delete [] status_;
          delete [] startSet_;
          delete [] next_;
          delete [] startColumn_;
          delete [] row_;
          delete [] element_;
          delete [] cost_;
          delete [] id_;
          delete [] dynamicStatus_;
          delete [] columnLower_;
          delete [] columnUpper_;
          objectiveOffset_ = rhs.objectiveOffset_;
          numberSets_ = rhs.numberSets_;
          numberActiveSets_ = rhs.numberActiveSets_;
          firstAvailable_ = rhs.firstAvailable_;
          firstAvailableBefore_ = rhs.firstAvailableBefore_;
          firstDynamic_ = rhs.firstDynamic_;
          lastDynamic_ = rhs.lastDynamic_;
          numberStaticRows_ = rhs.numberStaticRows_;
          numberElements_ = rhs.numberElements_;
          backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, lastDynamic_);
          keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
          toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
          fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + 1 - numberStaticRows_);
          lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
          upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
          status_ = ClpCopyOfArray(rhs.status_, static_cast<int>(2*numberSets_+4*sizeof(int)));
          model_ = rhs.model_;
          sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
          sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
          sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
          sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
          numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
          numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
          savedBestGubDual_ = rhs.savedBestGubDual_;
          savedBestSet_ = rhs.savedBestSet_;
          noCheck_ = rhs.noCheck_;
          infeasibilityWeight_ = rhs.infeasibilityWeight_;
          // Now secondary data
          numberGubColumns_ = rhs.numberGubColumns_;
          maximumGubColumns_ = rhs.maximumGubColumns_;
          maximumElements_ = rhs.maximumElements_;
          startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1);
          next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_);
          startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1);
          row_ = ClpCopyOfArray(rhs.row_, maximumElements_);
          element_ = ClpCopyOfArray(rhs.element_, maximumElements_);
          cost_ = ClpCopyOfArray(rhs.cost_, maximumGubColumns_);
          id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
          columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_);
          columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_);
          dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_);
     }
     return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpMatrixBase * ClpDynamicMatrix::clone() const
{
     return new ClpDynamicMatrix(*this);
}
// Partial pricing
void
ClpDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
                                 int & bestSequence, int & numberWanted)
{
     numberWanted = currentWanted_;
     assert(!model->rowScale());
     if (numberSets_) {
          // Do packed part before gub
          // always???
          //printf("normal packed price - start %d end %d (passed end %d, first %d)\n",
          ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
     } else {
          // no gub
          ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
          return;
     }
     if (numberWanted > 0) {
          // 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 bestDj;
          int numberRows = model->numberRows();
          int slackOffset = lastDynamic_ + numberRows;
          int structuralOffset = slackOffset + numberSets_;
          // If nothing found yet can go all the way to end
          int endAll = endG2;
          if (bestSequence < 0 && !startG2)
               endAll = numberSets_;
          if (bestSequence >= 0) {
               if (bestSequence != savedBestSequence_)
                    bestDj = fabs(reducedCost[bestSequence]); // dj from slacks or permanent
               else
                    bestDj = savedBestDj_;
          } else {
               bestDj = tolerance;
          }
          int saveSequence = bestSequence;
          double djMod = 0.0;
          double bestDjMod = 0.0;
          //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
          //     startG2,endG2,numberWanted);
          int bestSet = -1;
#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
          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;
               }
               int gubRow = toIndex_[iSet];
               if (gubRow >= 0) {
                    djMod = duals[gubRow+numberStaticRows_]; // have I got sign right?
               } else {
                    int iBasic = keyVariable_[iSet];
                    if (iBasic >= maximumGubColumns_) {
                         djMod = 0.0; // set not in
                    } else {
                         // get dj without
                         djMod = 0.0;
                         for (CoinBigIndex j = startColumn_[iBasic];
                                   j < startColumn_[iBasic+1]; 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 = slackOffset + iSet;
                                             bestDjMod = djMod;
                                             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 = slackOffset + iSet;
                                             bestDjMod = djMod;
                                             bestSet = iSet;
                                        } else {
                                             // just to make sure we don't exit before got something
                                             numberWanted++;
                                             abort();
                                        }
                                   }
                              }
                         }
                    }
               }
               int iSequence = startSet_[iSet];
               while (iSequence >= 0) {
                    DynamicStatus status = getDynamicStatus(iSequence);
                    if (status == atLowerBound || status == atUpperBound) {
                         double value = cost_[iSequence] - djMod;
                         for (CoinBigIndex 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)) {
				     if (false/*status == atLowerBound
						&&keyValue(iSet)<1.0e-7*/) {
				       // can't come in because
				       // of ones at ub
				       numberWanted++;
				     } else {
				     
                                        bestDj = value;
                                        bestSequence = structuralOffset + iSequence;
                                        bestDjMod = djMod;
                                        bestSet = iSet;
				     }
                                   } else {
                                        // just to make sure we don't exit before got something
                                        numberWanted++;
                                   }
                              }
                         }
                    }
                    iSequence = next_[iSequence]; //onto next in set
               }
               if (numberWanted <= 0) {
                    numberWanted = 0;
                    break;
               }
          }
          if (bestSequence != saveSequence) {
               savedBestGubDual_ = bestDjMod;
               savedBestDj_ = bestDj;
               savedBestSequence_ = bestSequence;
               savedBestSet_ = bestSet;
          }
          // See if may be finished
          if (!startG2 && bestSequence < 0)
               infeasibilityWeight_ = model_->infeasibilityCost();
          else if (bestSequence >= 0)
               infeasibilityWeight_ = -1.0;
     }
     currentWanted_ = numberWanted;
}
/* 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 *
ClpDynamicMatrix::rhsOffset(ClpSimplex * model, bool forceRefresh,
                            bool
#ifdef CLP_DEBUG
                            check
#endif
                           )
{
     // forceRefresh=true;printf("take out forceRefresh\n");
     if (!model_->numberIterations())
          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 iRow;
               int iSet;
               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;
               double objectiveOffset = 0.0;
               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 (columnLower_ || columnUpper_) {
                    double * solution = new double [numberGubColumns_];
                    for (iSet = 0; iSet < numberSets_; iSet++) {
                         int j = startSet_[iSet];
                         while (j >= 0) {
                              double value = 0.0;
                              if (getDynamicStatus(j) != inSmall) {
                                   if (getDynamicStatus(j) == atLowerBound) {
                                        if (columnLower_)
                                             value = columnLower_[j];
                                   } else if (getDynamicStatus(j) == atUpperBound) {
                                        value = columnUpper_[j];
                                   } else if (getDynamicStatus(j) == soloKey) {
                                        value = keyValue(iSet);
                                   }
                                   objectiveOffset += value * cost_[j];
                              }
                              solution[j] = value;
                              j = next_[j]; //onto next in set
                         }
                    }
                    // ones in gub and in small problem
                    for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
                         if (model_->getStatus(iColumn) != ClpSimplex::basic) {
                              int jFull = id_[iColumn-firstDynamic_];
                              solution[jFull] = smallSolution[iColumn];
                         }
                    }
                    for (iSet = 0; iSet < numberSets_; iSet++) {
                         int kRow = toIndex_[iSet];
                         if (kRow >= 0)
                              kRow += numberStaticRows_;
                         int j = startSet_[iSet];
                         while (j >= 0) {
                              double value = solution[j];
                              if (value) {
                                   for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) {
                                        int iRow = row_[k];
                                        rhs[iRow] -= element_[k] * value;
                                   }
                                   if (kRow >= 0)
                                        rhs[kRow] -= value;
                              }
                              j = next_[j]; //onto next in set
                         }
                    }
                    delete [] solution;
               } else {
                    // bounds
                    ClpSimplex::Status iStatus;
                    for (int iSet = 0; iSet < numberSets_; iSet++) {
                         int kRow = toIndex_[iSet];
                         if (kRow < 0) {
                              int iColumn = keyVariable_[iSet];
                              if (iColumn < maximumGubColumns_) {
                                   // key is not treated as basic
                                   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];
                                   if (b) {
                                        objectiveOffset += b * cost_[iColumn];
                                        for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
                                             int iRow = row_[j];
                                             rhs[iRow] -= element_[j] * b;
                                        }
                                   }
                              }
                         }
                    }
               }
               if (fabs(model->objectiveOffset() - objectiveOffset_ - objectiveOffset) > 1.0e-1)
                    printf("old offset %g, true %g\n", model->objectiveOffset(),
                           objectiveOffset_ - objectiveOffset);
               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 iSet;
               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;
               double objectiveOffset = 0.0;
               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 (columnLower_ || columnUpper_) {
                    double * solution = new double [numberGubColumns_];
                    for (iSet = 0; iSet < numberSets_; iSet++) {
                         int j = startSet_[iSet];
                         while (j >= 0) {
                              double value = 0.0;
                              if (getDynamicStatus(j) != inSmall) {
                                   if (getDynamicStatus(j) == atLowerBound) {
                                        if (columnLower_)
                                             value = columnLower_[j];
                                   } else if (getDynamicStatus(j) == atUpperBound) {
                                        value = columnUpper_[j];
					assert (value<1.0e30);
                                   } else if (getDynamicStatus(j) == soloKey) {
                                        value = keyValue(iSet);
				   }
                                   objectiveOffset += value * cost_[j];
                              }
                              solution[j] = value;
                              j = next_[j]; //onto next in set
                         }
                    }
                    // ones in gub and in small problem
                    for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
                         if (model_->getStatus(iColumn) != ClpSimplex::basic) {
                              int jFull = id_[iColumn-firstDynamic_];
                              solution[jFull] = smallSolution[iColumn];
                         }
                    }
                    for (iSet = 0; iSet < numberSets_; iSet++) {
                         int kRow = toIndex_[iSet];
                         if (kRow >= 0) 
			   kRow += numberStaticRows_;
			 int j = startSet_[iSet];
			 while (j >= 0) {
			   //? use DynamicStatus status = getDynamicStatus(j);
			   double value = solution[j];
			   if (value) {
			     for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) {
			       int iRow = row_[k];
			       rhsOffset_[iRow] -= element_[k] * value;
			     }
			     if (kRow >= 0)
			       rhsOffset_[kRow] -= value;
			   }
			   j = next_[j]; //onto next in set
			 }
                    }
                    delete [] solution;
               } else {
                    // bounds
                    ClpSimplex::Status iStatus;
                    for (int iSet = 0; iSet < numberSets_; iSet++) {
                         int kRow = toIndex_[iSet];
                         if (kRow < 0) {
                              int iColumn = keyVariable_[iSet];
                              if (iColumn < maximumGubColumns_) {
                                   // key is not treated as basic
                                   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];
                                   if (b) {
                                        objectiveOffset += b * cost_[iColumn];
                                        for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
                                             int iRow = row_[j];
                                             rhsOffset_[iRow] -= element_[j] * b;
                                        }
                                   }
                              }
                         }
                    }
               }
               model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
#ifdef CLP_DEBUG
               if (saveE) {
                    int iRow;
                    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
ClpDynamicMatrix::updatePivot(ClpSimplex * model, double oldInValue, double oldOutValue)
{

     // now update working model
     int sequenceIn = model->sequenceIn();
     int sequenceOut = model->sequenceOut();
     int numberColumns = model->numberColumns();
     if (sequenceIn != sequenceOut && sequenceIn < numberColumns)
          backToPivotRow_[sequenceIn] = model->pivotRow();
     if (sequenceIn >= firstDynamic_ && sequenceIn < numberColumns) {
          int bigSequence = id_[sequenceIn-firstDynamic_];
          if (getDynamicStatus(bigSequence) != inSmall) {
               firstAvailable_++;
               setDynamicStatus(bigSequence, inSmall);
          }
     }
     // make sure slack is synchronized
     if (sequenceIn >= numberColumns + numberStaticRows_) {
          int iDynamic = sequenceIn - numberColumns - numberStaticRows_;
          int iSet = fromIndex_[iDynamic];
          setStatus(iSet, model->getStatus(sequenceIn));
     }
     if (sequenceOut >= numberColumns + numberStaticRows_) {
          int iDynamic = sequenceOut - numberColumns - numberStaticRows_;
          int iSet = fromIndex_[iDynamic];
          // out may have gone through barrier - so check
          double valueOut = model->lowerRegion()[sequenceOut];
          if (fabs(valueOut - static_cast<double> (lowerSet_[iSet])) <
                    fabs(valueOut - static_cast<double> (upperSet_[iSet])))
               setStatus(iSet, ClpSimplex::atLowerBound);
          else
               setStatus(iSet, ClpSimplex::atUpperBound);
          if (lowerSet_[iSet] == upperSet_[iSet])
               setStatus(iSet, ClpSimplex::isFixed);
#if 0
          if (getStatus(iSet) != model->getStatus(sequenceOut))
               printf("** set %d status %d, var status %d\n", iSet,
                      getStatus(iSet), model->getStatus(sequenceOut));
#endif
     }
     ClpMatrixBase::updatePivot(model, oldInValue, oldOutValue);
#ifdef CLP_DEBUG
     char * inSmall = new char [numberGubColumns_];
     memset(inSmall, 0, numberGubColumns_);
     const double * solution = model->solutionRegion();
     for (int i = 0; i < numberGubColumns_; i++)
          if (getDynamicStatus(i) == ClpDynamicMatrix::inSmall)
               inSmall[i] = 1;
     for (int i = firstDynamic_; i < firstAvailable_; i++) {
          int k = id_[i-firstDynamic_];
          inSmall[k] = 0;
          //if (k>=23289&&k<23357&&solution[i])
          //printf("var %d (in small %d) has value %g\n",k,i,solution[i]);
     }
     for (int i = 0; i < numberGubColumns_; i++)
          assert (!inSmall[i]);
     delete [] inSmall;
#ifndef NDEBUG
     for (int i = 0; i < numberActiveSets_; i++) {
          int iSet = fromIndex_[i];
          assert (toIndex_[iSet] == i);
          //if (getStatus(iSet)!=model->getRowStatus(i+numberStaticRows_))
          //printf("*** set %d status %d, var status %d\n",iSet,
          //     getStatus(iSet),model->getRowStatus(i+numberStaticRows_));
          //assert (model->getRowStatus(i+numberStaticRows_)==getStatus(iSet));
          //if (iSet==1035) {
          //printf("rhs for set %d (%d) is %g %g - cost %g\n",iSet,i,model->lowerRegion(0)[i+numberStaticRows_],
          //     model->upperRegion(0)[i+numberStaticRows_],model->costRegion(0)[i+numberStaticRows_]);
          //}
     }
#endif
#endif
     if (numberStaticRows_+numberActiveSets_<model->numberRows())
       return 0;
     else
       return 1;
}
/*
     utility dual function for dealing with dynamic constraints
     mode=n see ClpGubMatrix.hpp for definition
     Remember to update here when settled down
*/
void
ClpDynamicMatrix::dualExpanded(ClpSimplex * model,
                               CoinIndexedVector * /*array*/,
                               double * /*other*/, int mode)
{
     switch (mode) {
          // modify costs before transposeUpdate
     case 0:
          break;
          // create duals for key variables (without check on dual infeasible)
     case 1:
          break;
          // as 1 but check slacks and compute djs
     case 2: {
          // do pivots
          int * pivotVariable = model->pivotVariable();
          int numberRows = numberStaticRows_ + numberActiveSets_;
          int numberColumns = model->numberColumns();
          for (int iRow = 0; iRow < numberRows; iRow++) {
               int iPivot = pivotVariable[iRow];
               if (iPivot < numberColumns)
                    backToPivotRow_[iPivot] = iRow;
          }
          if (noCheck_ >= 0) {
               if (infeasibilityWeight_ != model_->infeasibilityCost()) {
                    // don't bother checking
                    sumDualInfeasibilities_ = 100.0;
                    numberDualInfeasibilities_ = 1;
                    sumOfRelaxedDualInfeasibilities_ = 100.0;
                    return;
               }
          }
          // 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;
          double * dual = model->dualRowSolution();
          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;
          sumDualInfeasibilities_ = 0.0;
          numberDualInfeasibilities_ = 0;
          sumOfRelaxedDualInfeasibilities_ = 0.0;
          for (i = 0; i < numberSets_; i++) {
               double value = 0.0;
               int gubRow = toIndex_[i];
               if (gubRow < 0) {
                    int kColumn = keyVariable_[i];
                    if (kColumn < maximumGubColumns_) {
                         // 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 {
                    value = dual[gubRow+numberStaticRows_];
               }
               // Now subtract out from all
               int k = startSet_[i];
               while (k >= 0) {
                    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;
                         if (getDynamicStatus(k) == atLowerBound) {
                              if (djValue < -dualTolerance)
                                   infeasibility = -djValue - dualTolerance;
                         } else if (getDynamicStatus(k) == atUpperBound) {
                              // at upper bound
                              if (djValue > dualTolerance)
                                   infeasibility = djValue - dualTolerance;
                         }
                         if (infeasibility > 0.0) {
                              sumDualInfeasibilities_ += infeasibility;
                              if (infeasibility > relaxedTolerance)
                                   sumOfRelaxedDualInfeasibilities_ += infeasibility;
                              numberDualInfeasibilities_ ++;
                         }
                    }
                    k = next_[k]; //onto next in set
               }
          }
     }
     infeasibilityWeight_ = -1.0;
     break;
     // Report on infeasibilities of key variables
     case 3: {
          model->setSumDualInfeasibilities(model->sumDualInfeasibilities() +
                                           sumDualInfeasibilities_);
          model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() +
                                              numberDualInfeasibilities_);
          model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() +
                    sumOfRelaxedDualInfeasibilities_);
     }
     break;
     // modify costs before transposeUpdate for partial pricing
     case 4:
          break;
     }
}
/*
     general utility function for dealing with dynamic constraints
     mode=n see ClpGubMatrix.hpp for definition
     Remember to update here when settled down
*/
int
ClpDynamicMatrix::generalExpanded(ClpSimplex * model, int mode, int &number)
{
     int returnCode = 0;
#if 0 //ndef NDEBUG
     {
       int numberColumns = model->numberColumns();
       int numberRows = model->numberRows();
       int * pivotVariable = model->pivotVariable();
       if (pivotVariable&&model->numberIterations()) {
	 for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
	   assert (pivotVariable[i]==i+numberColumns);
	 }
       }
     }
#endif
     switch (mode) {
          // Fill in pivotVariable
     case 0: {
          // If no effective rhs - form it
          if (!rhsOffset_) {
               rhsOffset_ = new double[model->numberRows()];
               rhsOffset(model, true);
          }
          int i;
          int numberBasic = number;
          int numberColumns = model->numberColumns();
          // Use different array so can build from true pivotVariable_
          //int * pivotVariable = model->pivotVariable();
          int * pivotVariable = model->rowArray(0)->getIndices();
          for (i = 0; i < numberColumns; i++) {
               if (model->getColumnStatus(i) == ClpSimplex::basic)
                    pivotVariable[numberBasic++] = i;
          }
          number = numberBasic;
     }
     break;
     // Do initial extra rows + maximum basic
     case 2: {
          number = model->numberRows();
     }
     break;
     // Before normal replaceColumn
     case 3: {
          if (numberActiveSets_ + numberStaticRows_ == model_->numberRows()) {
               // no space - re-factorize
               returnCode = 4;
               number = -1; // say no need for normal replaceColumn
          }
     }
     break;
     // To see if can dual or primal
     case 4: {
          returnCode = 1;
     }
     break;
     // save status
     case 5: {
       memcpy(status_+numberSets_,status_,numberSets_);
       memcpy(status_+2*numberSets_,&numberActiveSets_,sizeof(int));
       memcpy(dynamicStatus_+maximumGubColumns_,
	      dynamicStatus_,maximumGubColumns_);
     }
     break;
     // restore status
     case 6: {
       memcpy(status_,status_+numberSets_,numberSets_);
       memcpy(&numberActiveSets_,status_+2*numberSets_,sizeof(int));
       memcpy(dynamicStatus_,dynamicStatus_+maximumGubColumns_,
	      maximumGubColumns_);
       initialProblem();
     }
     break;
     // unflag all variables
     case 8: {
          for (int i = 0; i < numberGubColumns_; i++) {
               if (flagged(i)) {
                    unsetFlagged(i);
                    returnCode++;
               }
          }
     }
     break;
     // redo costs in primal
     case 9: {
          double * cost = model->costRegion();
          double * solution = model->solutionRegion();
          double * columnLower = model->lowerRegion();
          double * columnUpper = model->upperRegion();
          int i;
          bool doCosts = (number & 4) != 0;
          bool doBounds = (number & 1) != 0;
          for ( i = firstDynamic_; i < firstAvailable_; i++) {
               int jColumn = id_[i-firstDynamic_];
               if (doBounds) {
                    if (!columnLower_ && !columnUpper_) {
                         columnLower[i] = 0.0;
                         columnUpper[i] = COIN_DBL_MAX;
                    }  else {
                         if (columnLower_)
                              columnLower[i] = columnLower_[jColumn];
                         else
                              columnLower[i] = 0.0;
                         if (columnUpper_)
                              columnUpper[i] = columnUpper_[jColumn];
                         else
                              columnUpper[i] = COIN_DBL_MAX;
                    }
               }
               if (doCosts) {
                    cost[i] = cost_[jColumn];
                    // Original bounds
                    if (model->nonLinearCost())
                         model->nonLinearCost()->setOne(i, solution[i],
                                                        this->columnLower(jColumn),
                                                        this->columnUpper(jColumn), cost_[jColumn]);
               }
          }
          // and active sets
          for (i = 0; i < numberActiveSets_; i++ ) {
               int iSet = fromIndex_[i];
               int iSequence = lastDynamic_ + numberStaticRows_ + i;
               if (doBounds) {
                    if (lowerSet_[iSet] > -1.0e20)
                         columnLower[iSequence] = lowerSet_[iSet];
                    else
                         columnLower[iSequence] = -COIN_DBL_MAX;
                    if (upperSet_[iSet] < 1.0e20)
                         columnUpper[iSequence] = upperSet_[iSet];
                    else
                         columnUpper[iSequence] = COIN_DBL_MAX;
               }
               if (doCosts) {
                    if (model->nonLinearCost()) {
                         double trueLower;
                         if (lowerSet_[iSet] > -1.0e20)
                              trueLower = lowerSet_[iSet];
                         else
                              trueLower = -COIN_DBL_MAX;
                         double trueUpper;
                         if (upperSet_[iSet] < 1.0e20)
                              trueUpper = upperSet_[iSet];
                         else
                              trueUpper = COIN_DBL_MAX;
                         model->nonLinearCost()->setOne(iSequence, solution[iSequence],
                                                        trueLower, trueUpper, 0.0);
                    }
               }
          }
     }
     break;
     // return 1 if there may be changing bounds on variable (column generation)
     case 10: {
          // return 1 as bounds on rhs will change
          returnCode = 1;
     }
     break;
     // make sure set is clean
     case 7: {
          // first flag
          if (number >= firstDynamic_ && number < lastDynamic_) {
	    int sequence = id_[number-firstDynamic_];
	    setFlagged(sequence);
	    //model->clearFlagged(number);
	  } else if (number>=model_->numberColumns()+numberStaticRows_) {
	    // slack
	    int iSet = fromIndex_[number-model_->numberColumns()-
				  numberStaticRows_];
	    setFlaggedSlack(iSet);
	    //model->clearFlagged(number);
	  }
     }
     case 11: {
          //int sequenceIn = model->sequenceIn();
          if (number >= firstDynamic_ && number < lastDynamic_) {
	    //assert (number == model->sequenceIn());
               // take out variable (but leave key)
               double * cost = model->costRegion();
               double * columnLower = model->lowerRegion();
               double * columnUpper = model->upperRegion();
               double * solution = model->solutionRegion();
               int * length = matrix_->getMutableVectorLengths();
#ifndef NDEBUG
               if (length[number]) {
                    int * row = matrix_->getMutableIndices();
                    CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
                    int iRow = row[startColumn[number] + length[number] - 1];
                    int which = iRow - numberStaticRows_;
                    assert (which >= 0);
                    int iSet = fromIndex_[which];
                    assert (toIndex_[iSet] == which);
               }
#endif
               // no need firstAvailable_--;
               solution[firstAvailable_] = 0.0;
               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);
               columnLower[firstAvailable_] = 0.0;
               columnUpper[firstAvailable_] = COIN_DBL_MAX;

               // not really in small problem
               int iBig = id_[number-firstDynamic_];
               if (model->getStatus(number) == ClpSimplex::atLowerBound) {
                    setDynamicStatus(iBig, atLowerBound);
                    if (columnLower_)
                         modifyOffset(number, columnLower_[iBig]);
               } else {
                    setDynamicStatus(iBig, atUpperBound);
                    modifyOffset(number, columnUpper_[iBig]);
               }
	  } else if (number>=model_->numberColumns()+numberStaticRows_) {
	    // slack
	    int iSet = fromIndex_[number-model_->numberColumns()-
				  numberStaticRows_];
	    printf("what now - set %d\n",iSet);
          }
     }
     break;
     default:
          break;
     }
     return returnCode;
}
/* Purely for column generation and similar ideas.  Allows
   matrix and any bounds or costs to be updated (sensibly).
   Returns non-zero if any changes.
*/
int
ClpDynamicMatrix::refresh(ClpSimplex * model)
{
     // If at beginning then create initial problem
     if (firstDynamic_ == firstAvailable_) {
          initialProblem();
          return 1;
     } else if (!model->nonLinearCost()) {
          // will be same as last time
          return 1;
     }
#ifndef NDEBUG
     {
       int numberColumns = model->numberColumns();
       int numberRows = model->numberRows();
       int * pivotVariable = model->pivotVariable();
       for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
	 assert (pivotVariable[i]==i+numberColumns);
       }
     }
#endif
     // lookup array
     int * active = new int [numberActiveSets_];
     CoinZeroN(active, numberActiveSets_);
     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 * columnLower = model->lowerRegion();
     double * columnUpper = model->upperRegion();
     CoinBigIndex numberElements = startColumn[firstDynamic_];
     // first just do lookup and basic stuff
     int currentNumber = firstAvailable_;
     firstAvailable_ = firstDynamic_;
     double objectiveChange = 0.0;
     double * solution = model->solutionRegion();
     int currentNumberActiveSets = numberActiveSets_;
     for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
          int iRow = row[startColumn[iColumn] + length[iColumn] - 1];
          int which = iRow - numberStaticRows_;
          assert (which >= 0);
          int iSet = fromIndex_[which];
          assert (toIndex_[iSet] == which);
          if (model->getStatus(iColumn) == ClpSimplex::basic) {
               active[which]++;
               // may as well make key
               keyVariable_[iSet] = id_[iColumn-firstDynamic_];
          }
     }
     int i;
     numberActiveSets_ = 0;
     int numberDeleted = 0;
     for (i = 0; i < currentNumberActiveSets; i++) {
          int iRow = i + numberStaticRows_;
          int numberActive = active[i];
          int iSet = fromIndex_[i];
          if (model->getRowStatus(iRow) == ClpSimplex::basic) {
               numberActive++;
               // may as well make key
               keyVariable_[iSet] = maximumGubColumns_ + iSet;
          }
          if (numberActive > 1) {
               // keep
               active[i] = numberActiveSets_;
               numberActiveSets_++;
          } else {
               active[i] = -1;
          }
     }

     for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
          int iRow = row[startColumn[iColumn] + length[iColumn] - 1];
          int which = iRow - numberStaticRows_;
          int jColumn = id_[iColumn-firstDynamic_];
          if (model->getStatus(iColumn) == ClpSimplex::atLowerBound ||
                    model->getStatus(iColumn) == ClpSimplex::atUpperBound) {
               double value = solution[iColumn];
               bool toLowerBound = true;
	       assert (jColumn>=0);assert (iColumn>=0);
               if (columnUpper_) {
                    if (!columnLower_) {
                         if (fabs(value - columnUpper_[jColumn]) < fabs(value))
                              toLowerBound = false;
                    } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[jColumn])) {
                         toLowerBound = false;
                    }
               }
               if (toLowerBound) {
                    // throw out to lower bound
                    if (columnLower_) {
                         setDynamicStatus(jColumn, atLowerBound);
                         // treat solution as if exactly at a bound
                         double value = columnLower[iColumn];
                         objectiveChange += cost[iColumn] * value;
                    } else {
                         setDynamicStatus(jColumn, atLowerBound);
                    }
               } else {
                    // throw out to upper bound
                    setDynamicStatus(jColumn, atUpperBound);
                    double value = columnUpper[iColumn];
                    objectiveChange += cost[iColumn] * value;
               }
               numberDeleted++;
          } else {
               assert(model->getStatus(iColumn) != ClpSimplex::superBasic); // deal with later
               int iPut = active[which];
               if (iPut >= 0) {
                    // move
                    id_[firstAvailable_-firstDynamic_] = jColumn;
                    int numberThis = startColumn_[jColumn+1] - startColumn_[jColumn] + 1;
                    length[firstAvailable_] = numberThis;
                    cost[firstAvailable_] = cost[iColumn];
                    columnLower[firstAvailable_] = columnLower[iColumn];
                    columnUpper[firstAvailable_] = columnUpper[iColumn];
                    model->nonLinearCost()->setOne(firstAvailable_, solution[iColumn], 0.0, COIN_DBL_MAX,
                                                   cost_[jColumn]);
                    CoinBigIndex base = startColumn_[jColumn];
                    for (int j = 0; j < numberThis - 1; j++) {
                         row[numberElements] = row_[base+j];
                         element[numberElements++] = element_[base+j];
                    }
                    row[numberElements] = iPut + numberStaticRows_;
                    element[numberElements++] = 1.0;
                    model->setStatus(firstAvailable_, model->getStatus(iColumn));
                    solution[firstAvailable_] = solution[iColumn];
                    firstAvailable_++;
                    startColumn[firstAvailable_] = numberElements;
               }
          }
     }
     // 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);
          columnLower[iColumn] = 0.0;
          columnUpper[iColumn] = COIN_DBL_MAX;
     }
     // move rhs and set rest to infinite
     numberActiveSets_ = 0;
     for (i = 0; i < currentNumberActiveSets; i++) {
          int iSet = fromIndex_[i];
          assert (toIndex_[iSet] == i);
          int iRow = i + numberStaticRows_;
          int iPut = active[i];
          if (iPut >= 0) {
               int in = iRow + numberColumns;
               int out = iPut + numberColumns + numberStaticRows_;
               solution[out] = solution[in];
               columnLower[out] = columnLower[in];
               columnUpper[out] = columnUpper[in];
               cost[out] = cost[in];
               double trueLower;
               if (lowerSet_[iSet] > -1.0e20)
                    trueLower = lowerSet_[iSet];
               else
                    trueLower = -COIN_DBL_MAX;
               double trueUpper;
               if (upperSet_[iSet] < 1.0e20)
                    trueUpper = upperSet_[iSet];
               else
                    trueUpper = COIN_DBL_MAX;
               model->nonLinearCost()->setOne(out, solution[out],
                                              trueLower, trueUpper, 0.0);
               model->setStatus(out, model->getStatus(in));
               toIndex_[iSet] = numberActiveSets_;
               rhsOffset_[numberActiveSets_+numberStaticRows_] = rhsOffset_[i+numberStaticRows_];
               fromIndex_[numberActiveSets_++] = fromIndex_[i];
          } else {
               // adjust offset
               // put out as key
               int jColumn = keyVariable_[iSet];
               toIndex_[iSet] = -1;
               if (jColumn < maximumGubColumns_) {
                    setDynamicStatus(jColumn, soloKey);
                    double value = keyValue(iSet);
                    objectiveChange += cost_[jColumn] * value;
                    modifyOffset(jColumn, -value);
               }
          }
     }
     model->setObjectiveOffset(model->objectiveOffset() - objectiveChange);
     delete [] active;
     for (i = numberActiveSets_; i < currentNumberActiveSets; i++) {
          int iSequence = i + numberStaticRows_ + numberColumns;
          solution[iSequence] = 0.0;
          columnLower[iSequence] = -COIN_DBL_MAX;
          columnUpper[iSequence] = COIN_DBL_MAX;
          cost[iSequence] = 0.0;
          model->nonLinearCost()->setOne(iSequence, solution[iSequence],
                                         columnLower[iSequence],
                                         columnUpper[iSequence], 0.0);
          model->setStatus(iSequence, ClpSimplex::basic);
          rhsOffset_[i+numberStaticRows_] = 0.0;
     }
     if (currentNumberActiveSets != numberActiveSets_ || numberDeleted) {
          // deal with pivotVariable
          int * pivotVariable = model->pivotVariable();
          int sequence = firstDynamic_;
          int iRow = 0;
          int base1 = firstDynamic_;
          int base2 = lastDynamic_;
          int base3 = numberColumns + numberStaticRows_;
          int numberRows = numberStaticRows_ + currentNumberActiveSets;
          while (sequence < firstAvailable_) {
               int iPivot = pivotVariable[iRow];
               while (iPivot < base1 || (iPivot >= base2 && iPivot < base3)) {
                    iPivot = pivotVariable[++iRow];
               }
               pivotVariable[iRow++] = sequence;
               sequence++;
          }
          // move normal basic ones down
          int iPut = iRow;
          for (; iRow < numberRows; iRow++) {
               int iPivot = pivotVariable[iRow];
               if (iPivot < base1 || (iPivot >= base2 && iPivot < base3)) {
                    pivotVariable[iPut++] = iPivot;
               }
          }
          // and basic others
          for (i = 0; i < numberActiveSets_; i++) {
               if (model->getRowStatus(i + numberStaticRows_) == ClpSimplex::basic) {
                    pivotVariable[iPut++] = i + base3;
               }
          }
	  if (iPut<numberStaticRows_+numberActiveSets_) {
	    printf("lost %d sets\n",
		   numberStaticRows_+numberActiveSets_-iPut);
	    iPut = numberStaticRows_+numberActiveSets_;
	  }
          for (i = numberActiveSets_; i < currentNumberActiveSets; i++) {
               pivotVariable[iPut++] = i + base3;
          }
          //assert (iPut == numberRows);
     }
#ifdef CLP_DEBUG
#if 0
     printf("row for set 244 is %d, row status %d value %g ", toIndex_[244], status_[244],
            keyValue(244));
     int jj = startSet_[244];
     while (jj >= 0) {
          printf("var %d status %d ", jj, dynamicStatus_[jj]);
          jj = next_[jj];
     }
     printf("\n");
#endif
#ifdef CLP_DEBUG
     {
          // debug coding to analyze set
          int i;
          int kSet = -6;
          if (kSet >= 0) {
               int * back = new int [numberGubColumns_];
               for (i = 0; i < numberGubColumns_; i++)
                    back[i] = -1;
               for (i = firstDynamic_; i < firstAvailable_; i++)
                    back[id_[i-firstDynamic_]] = i;
               int sequence = startSet_[kSet];
               if (toIndex_[kSet] < 0) {
                    printf("Not in - Set %d - slack status %d, key %d\n", kSet, status_[kSet], keyVariable_[kSet]);
                    while (sequence >= 0) {
                         printf("( %d status %d ) ", sequence, getDynamicStatus(sequence));
                         sequence = next_[sequence];
                    }
               } else {
                    int iRow = numberStaticRows_ + toIndex_[kSet];
                    printf("In - Set %d - slack status %d, key %d offset %g slack %g\n", kSet, status_[kSet], keyVariable_[kSet],
                           rhsOffset_[iRow], model->solutionRegion(0)[iRow]);
                    while (sequence >= 0) {
                         int iBack = back[sequence];
                         printf("( %d status %d value %g) ", sequence, getDynamicStatus(sequence), model->solutionRegion()[iBack]);
                         sequence = next_[sequence];
                    }
               }
               printf("\n");
               delete [] back;
          }
     }
#endif
     int n = numberActiveSets_;
     for (i = 0; i < numberSets_; i++) {
          if (toIndex_[i] < 0) {
               //assert(keyValue(i)>=lowerSet_[i]&&keyValue(i)<=upperSet_[i]);
               n++;
          }
	  int k=0;
	  for (int j=startSet_[i];j<startSet_[i+1];j++) {
	    if (getDynamicStatus(j)==soloKey)
	      k++;
	  }
	  assert (k<2);
     }
     assert (n == numberSets_);
#endif
     return 1;
}
void
ClpDynamicMatrix::times(double scalar,
                        const double * x, double * y) const
{
     if (model_->specialOptions() != 16) {
          ClpPackedMatrix::times(scalar, x, y);
     } else {
          int iRow;
          const double * element =  matrix_->getElements();
          const int * row = matrix_->getIndices();
          const CoinBigIndex * startColumn = matrix_->getVectorStarts();
          const int * length = matrix_->getVectorLengths();
          int * pivotVariable = model_->pivotVariable();
          for (iRow = 0; iRow < numberStaticRows_ + numberActiveSets_; iRow++) {
               y[iRow] -= scalar * rhsOffset_[iRow];
               int iColumn = pivotVariable[iRow];
               if (iColumn < lastDynamic_) {
                    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];
                         }
                    }
               }
          }
     }
}
// Modifies rhs offset
void
ClpDynamicMatrix::modifyOffset(int sequence, double amount)
{
     if (amount) {
          assert (rhsOffset_);
          CoinBigIndex j;
          for (j = startColumn_[sequence]; j < startColumn_[sequence+1]; j++) {
               int iRow = row_[j];
               rhsOffset_[iRow] += amount * element_[j];
          }
     }
}
// Gets key value when none in small
double
ClpDynamicMatrix::keyValue(int iSet) const
{
     double value = 0.0;
     if (toIndex_[iSet] < 0) {
          int key = keyVariable_[iSet];
          if (key < maximumGubColumns_) {
               if (getStatus(iSet) == ClpSimplex::atLowerBound)
                    value = lowerSet_[iSet];
               else
                    value = upperSet_[iSet];
               int numberKey = 0;
               int j = startSet_[iSet];
               while (j >= 0) {
                    DynamicStatus status = getDynamicStatus(j);
                    assert (status != inSmall);
                    if (status == soloKey) {
                         numberKey++;
                    } else if (status == atUpperBound) {
                         value -= columnUpper_[j];
                    } else if (columnLower_) {
                         value -= columnLower_[j];
                    }
                    j = next_[j]; //onto next in set
               }
               assert (numberKey == 1);
          } else {
               int j = startSet_[iSet];
               while (j >= 0) {
                    DynamicStatus status = getDynamicStatus(j);
                    assert (status != inSmall);
                    assert (status != soloKey);
                    if (status == atUpperBound) {
                         value += columnUpper_[j];
                    } else if (columnLower_) {
                         value += columnLower_[j];
                    }
                    j = next_[j]; //onto next in set
               }
#if 0
	       // slack must be feasible
	       double oldValue=value;
	       value = CoinMax(value,lowerSet_[iSet]);
               value = CoinMin(value,upperSet_[iSet]);
	       if (value!=oldValue)
		 printf("using %g (not %g) for slack on set %d (%g,%g)\n",
			value,oldValue,iSet,lowerSet_[iSet],upperSet_[iSet]);
#endif
            }
     }
     return value;
}
// Switches off dj checking each factorization (for BIG models)
void
ClpDynamicMatrix::switchOffCheck()
{
     noCheck_ = 0;
     infeasibilityWeight_ = 0.0;
}
/* Creates a variable.  This is called after partial pricing and may modify matrix.
   May update bestSequence.
*/
void
ClpDynamicMatrix::createVariable(ClpSimplex * model, int & bestSequence)
{
     int numberRows = model->numberRows();
     int slackOffset = lastDynamic_ + numberRows;
     int structuralOffset = slackOffset + numberSets_;
     int bestSequence2 = savedBestSequence_ - structuralOffset;
     if (bestSequence >= slackOffset) {
          double * columnLower = model->lowerRegion();
          double * columnUpper = model->upperRegion();
          double * solution = model->solutionRegion();
          double * reducedCost = model->djRegion();
          const double * duals = model->dualRowSolution();
          if (toIndex_[savedBestSet_] < 0) {
               // need to put key into basis
               int newRow = numberActiveSets_ + numberStaticRows_;
               model->dualRowSolution()[newRow] = savedBestGubDual_;
               double valueOfKey = keyValue(savedBestSet_); // done before toIndex_ set
               toIndex_[savedBestSet_] = numberActiveSets_;
               fromIndex_[numberActiveSets_++] = savedBestSet_;
               int iSequence = lastDynamic_ + newRow;
               // we need to get lower and upper correct
               double shift = 0.0;
               int j = startSet_[savedBestSet_];
               while (j >= 0) {
                    if (getDynamicStatus(j) == atUpperBound)
                         shift += columnUpper_[j];
                    else if (getDynamicStatus(j) == atLowerBound && columnLower_)
                         shift += columnLower_[j];
                    j = next_[j]; //onto next in set
               }
               if (lowerSet_[savedBestSet_] > -1.0e20)
                    columnLower[iSequence] = lowerSet_[savedBestSet_];
               else
                    columnLower[iSequence] = -COIN_DBL_MAX;
               if (upperSet_[savedBestSet_] < 1.0e20)
                    columnUpper[iSequence] = upperSet_[savedBestSet_];
               else
                    columnUpper[iSequence] = COIN_DBL_MAX;
#ifdef CLP_DEBUG
               if (model_->logLevel() == 63) {
                    printf("%d in in set %d, key is %d rhs %g %g - keyvalue %g\n",
                           bestSequence2, savedBestSet_, keyVariable_[savedBestSet_],
                           columnLower[iSequence], columnUpper[iSequence], valueOfKey);
                    int j = startSet_[savedBestSet_];
                    while (j >= 0) {
                         if (getDynamicStatus(j) == atUpperBound)
                              printf("%d atup ", j);
                         else if (getDynamicStatus(j) == atLowerBound)
                              printf("%d atlo ", j);
                         else if (getDynamicStatus(j) == soloKey)
                              printf("%d solo ", j);
                         else
                              abort();
                         j = next_[j]; //onto next in set
                    }
                    printf("\n");
               }
#endif
               if (keyVariable_[savedBestSet_] < maximumGubColumns_) {
                    // slack not key
                    model_->pivotVariable()[newRow] = firstAvailable_;
                    backToPivotRow_[firstAvailable_] = newRow;
                    model->setStatus(iSequence, getStatus(savedBestSet_));
                    model->djRegion()[iSequence] = savedBestGubDual_;
                    solution[iSequence] = valueOfKey;
		    // create variable and pivot in
                    int key = keyVariable_[savedBestSet_];
                    setDynamicStatus(key, inSmall);
                    double * element =  matrix_->getMutableElements();
                    int * row = matrix_->getMutableIndices();
                    CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
                    int * length = matrix_->getMutableVectorLengths();
                    CoinBigIndex numberElements = startColumn[firstAvailable_];
                    int numberThis = startColumn_[key+1] - startColumn_[key] + 1;
                    if (numberElements + numberThis > numberElements_) {
                         // need to redo
                         numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
                         matrix_->reserve(lastDynamic_, 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_[key];
                    CoinBigIndex base = startColumn_[key];
                    for (int j = 0; j < numberThis - 1; j++) {
                         row[numberElements] = row_[base+j];
                         element[numberElements++] = element_[base+j];
                    }
                    row[numberElements] = newRow;
                    element[numberElements++] = 1.0;
                    id_[firstAvailable_-firstDynamic_] = key;
                    model->setObjectiveOffset(model->objectiveOffset() + cost_[key]*
                                              valueOfKey);
                    model->solutionRegion()[firstAvailable_] = valueOfKey;
                    model->setStatus(firstAvailable_, ClpSimplex::basic);
                    // ***** need to adjust effective rhs
                    if (!columnLower_ && !columnUpper_) {
                         columnLower[firstAvailable_] = 0.0;
                         columnUpper[firstAvailable_] = COIN_DBL_MAX;
                    }  else {
                         if (columnLower_)
                              columnLower[firstAvailable_] = columnLower_[key];
                         else
                              columnLower[firstAvailable_] = 0.0;
                         if (columnUpper_)
                              columnUpper[firstAvailable_] = columnUpper_[key];
                         else
                              columnUpper[firstAvailable_] = COIN_DBL_MAX;
                    }
                    model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
                                                   columnLower[firstAvailable_],
                                                   columnUpper[firstAvailable_], cost_[key]);
                    startColumn[firstAvailable_+1] = numberElements;
                    reducedCost[firstAvailable_] = 0.0;
                    modifyOffset(key, valueOfKey);
                    rhsOffset_[newRow] = -shift; // sign?
#ifdef CLP_DEBUG
                    model->rowArray(1)->checkClear();
#endif
                    // now pivot in
                    unpack(model, model->rowArray(1), firstAvailable_);
                    model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(1));
                    double alpha = model->rowArray(1)->denseVector()[newRow];
                    int updateStatus =
                         model->factorization()->replaceColumn(model,
                                   model->rowArray(2),
                                   model->rowArray(1),
                                   newRow, alpha);
                    model->rowArray(1)->clear();
                    if (updateStatus) {
                         if (updateStatus == 3) {
                              // out of memory
                              // increase space if not many iterations
                              if (model->factorization()->pivots() <
                                        0.5 * model->factorization()->maximumPivots() &&
                                        model->factorization()->pivots() < 400)
                                   model->factorization()->areaFactor(
                                        model->factorization()->areaFactor() * 1.1);
                         } else {
                              printf("Bad returncode %d from replaceColumn\n", updateStatus);
                         }
                         bestSequence = -1;
                         return;
                    }
                    // firstAvailable_ only finally updated if good pivot (in updatePivot)
                    // otherwise it reverts to firstAvailableBefore_
                    firstAvailable_++;
               } else {
                    // slack key
                    model->setStatus(iSequence, ClpSimplex::basic);
                    model->djRegion()[iSequence] = 0.0;
                    solution[iSequence] = valueOfKey+shift;
                    rhsOffset_[newRow] = -shift; // sign?
               }
               // correct slack
               model->costRegion()[iSequence] = 0.0;
               model->nonLinearCost()->setOne(iSequence, solution[iSequence], columnLower[iSequence],
                                              columnUpper[iSequence], 0.0);
          }
          if (savedBestSequence_ >= structuralOffset) {
               // recompute dj and create
               double value = cost_[bestSequence2] - savedBestGubDual_;
               for (CoinBigIndex jBigIndex = startColumn_[bestSequence2];
                         jBigIndex < startColumn_[bestSequence2+1]; jBigIndex++) {
                    int jRow = row_[jBigIndex];
                    value -= duals[jRow] * element_[jBigIndex];
               }
               int gubRow = toIndex_[savedBestSet_] + numberStaticRows_;
               double * element =  matrix_->getMutableElements();
               int * row = matrix_->getMutableIndices();
               CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
               int * length = matrix_->getMutableVectorLengths();
               CoinBigIndex numberElements = startColumn[firstAvailable_];
               int numberThis = startColumn_[bestSequence2+1] - startColumn_[bestSequence2] + 1;
               if (numberElements + numberThis > numberElements_) {
                    // need to redo
                    numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
                    matrix_->reserve(lastDynamic_, 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_[bestSequence2];
               CoinBigIndex base = startColumn_[bestSequence2];
               for (int j = 0; j < numberThis - 1; j++) {
                    row[numberElements] = row_[base+j];
                    element[numberElements++] = element_[base+j];
               }
               row[numberElements] = gubRow;
               element[numberElements++] = 1.0;
               id_[firstAvailable_-firstDynamic_] = bestSequence2;
               //printf("best %d\n",bestSequence2);
               model->solutionRegion()[firstAvailable_] = 0.0;
	       model->clearFlagged(firstAvailable_);
               if (!columnLower_ && !columnUpper_) {
                    model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
                    columnLower[firstAvailable_] = 0.0;
                    columnUpper[firstAvailable_] = COIN_DBL_MAX;
               }  else {
                    DynamicStatus status = getDynamicStatus(bestSequence2);
                    if (columnLower_)
                         columnLower[firstAvailable_] = columnLower_[bestSequence2];
                    else
                         columnLower[firstAvailable_] = 0.0;
                    if (columnUpper_)
                         columnUpper[firstAvailable_] = columnUpper_[bestSequence2];
                    else
                         columnUpper[firstAvailable_] = COIN_DBL_MAX;
                    if (status == atLowerBound) {
                         solution[firstAvailable_] = columnLower[firstAvailable_];
                         model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
                    } else {
                         solution[firstAvailable_] = columnUpper[firstAvailable_];
                         model->setStatus(firstAvailable_, ClpSimplex::atUpperBound);
                    }
               }
               model->setObjectiveOffset(model->objectiveOffset() + cost_[bestSequence2]*
                                         solution[firstAvailable_]);
               model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
                                              columnLower[firstAvailable_],
                                              columnUpper[firstAvailable_], cost_[bestSequence2]);
               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,savedBestGubDual_);
               reducedCost[bestSequence] = value;
          } else {
               // slack - row must just have been created
               assert (toIndex_[savedBestSet_] == numberActiveSets_ - 1);
               int newRow = numberStaticRows_ + numberActiveSets_ - 1;
               bestSequence = lastDynamic_ + newRow;
               reducedCost[bestSequence] = savedBestGubDual_;
          }
     }
     // clear for next iteration
     savedBestSequence_ = -1;
}
// Returns reduced cost of a variable
double
ClpDynamicMatrix::reducedCost(ClpSimplex * model, int sequence) const
{
     int numberRows = model->numberRows();
     int slackOffset = lastDynamic_ + numberRows;
     if (sequence < slackOffset)
          return model->djRegion()[sequence];
     else
          return savedBestDj_;
}
// Does gub crash
void
ClpDynamicMatrix::gubCrash()
{
     // Do basis - cheapest or slack if feasible
     int longestSet = 0;
     int iSet;
     for (iSet = 0; iSet < numberSets_; iSet++) {
          int n = 0;
          int j = startSet_[iSet];
          while (j >= 0) {
               n++;
               j = next_[j];
          }
          longestSet = CoinMax(longestSet, n);
     }
     double * upper = new double[longestSet+1];
     double * cost = new double[longestSet+1];
     double * lower = new double[longestSet+1];
     double * solution = new double[longestSet+1];
     int * back = new int[longestSet+1];
     double tolerance = model_->primalTolerance();
     double objectiveOffset = 0.0;
     for (iSet = 0; iSet < numberSets_; iSet++) {
          int iBasic = -1;
          double value = 0.0;
          // find cheapest
          int numberInSet = 0;
          int j = startSet_[iSet];
          while (j >= 0) {
               if (!columnLower_)
                    lower[numberInSet] = 0.0;
               else
                    lower[numberInSet] = columnLower_[j];
               if (!columnUpper_)
                    upper[numberInSet] = COIN_DBL_MAX;
               else
                    upper[numberInSet] = columnUpper_[j];
               back[numberInSet++] = j;
               j = next_[j];
          }
          CoinFillN(solution, numberInSet, 0.0);
          // and slack
          iBasic = numberInSet;
          solution[iBasic] = -value;
          lower[iBasic] = -upperSet_[iSet];
          upper[iBasic] = -lowerSet_[iSet];
          int kphase;
          if (value < lowerSet_[iSet] - tolerance || value > upperSet_[iSet] + tolerance) {
               // infeasible
               kphase = 0;
               // remember bounds are flipped so opposite to natural
               if (value < lowerSet_[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_[back[j]];
                    }
                    // 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 (columnLower_ || columnUpper_) {
               for (int j = 0; j < numberInSet; j++) {
                    if (j != iBasic) {
                         objectiveOffset += solution[j] * cost[j];
                         if (columnLower_ && columnUpper_) {
                              if (fabs(solution[j] - columnLower_[back[j]]) >
                                        fabs(solution[j] - columnUpper_[back[j]]))
                                   setDynamicStatus(back[j], atUpperBound);
                         } else if (columnUpper_ && solution[j] > 0.0) {
                              setDynamicStatus(back[j], atUpperBound);
                         } else {
                              setDynamicStatus(back[j], atLowerBound);
                              assert(!solution[j]);
                         }
                    }
               }
          }
          // convert iBasic back and do bounds
          if (iBasic == numberInSet) {
               // slack basic
               setStatus(iSet, ClpSimplex::basic);
               iBasic = iSet + maximumGubColumns_;
          } else {
               iBasic = back[iBasic];
               setDynamicStatus(iBasic, soloKey);
               // 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();
          }
          keyVariable_[iSet] = iBasic;
     }
     model_->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
     delete [] lower;
     delete [] solution;
     delete [] upper;
     delete [] cost;
     delete [] back;
     // make sure matrix is in good shape
     matrix_->orderMatrix();
}
// Populates initial matrix from dynamic status
void
ClpDynamicMatrix::initialProblem()
{
     int iSet;
     double * element =  matrix_->getMutableElements();
     int * row = matrix_->getMutableIndices();
     CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
     int * length = matrix_->getMutableVectorLengths();
     double * cost = model_->objective();
     double * solution = model_->primalColumnSolution();
     double * columnLower = model_->columnLower();
     double * columnUpper = model_->columnUpper();
     double * rowSolution = model_->primalRowSolution();
     double * rowLower = model_->rowLower();
     double * rowUpper = model_->rowUpper();
     CoinBigIndex numberElements = startColumn[firstDynamic_];

     firstAvailable_ = firstDynamic_;
     numberActiveSets_ = 0;
     for (iSet = 0; iSet < numberSets_; iSet++) {
          toIndex_[iSet] = -1;
          int numberActive = 0;
          int whichKey = -1;
          if (getStatus(iSet) == ClpSimplex::basic) {
               whichKey = maximumGubColumns_ + iSet;
	       numberActive = 1;
          } else {
               whichKey = -1;
	  }
          int j = startSet_[iSet];
          while (j >= 0) {
               assert (getDynamicStatus(j) != soloKey || whichKey == -1);
               if (getDynamicStatus(j) == inSmall) {
                    numberActive++;
               } else if (getDynamicStatus(j) == soloKey) {
                    whichKey = j;
                    numberActive++;
	       }
               j = next_[j]; //onto next in set
          }
          if (numberActive > 1) {
               int iRow = numberActiveSets_ + numberStaticRows_;
               rowSolution[iRow] = 0.0;
               double lowerValue;
               if (lowerSet_[iSet] > -1.0e20)
                    lowerValue = lowerSet_[iSet];
               else
                    lowerValue = -COIN_DBL_MAX;
               double upperValue;
               if (upperSet_[iSet] < 1.0e20)
                    upperValue = upperSet_[iSet];
               else
                    upperValue = COIN_DBL_MAX;
               rowLower[iRow] = lowerValue;
               rowUpper[iRow] = upperValue;
               if (getStatus(iSet) == ClpSimplex::basic) {
                    model_->setRowStatus(iRow, ClpSimplex::basic);
                    rowSolution[iRow] = 0.0;
               } else if (getStatus(iSet) == ClpSimplex::atLowerBound) {
                    model_->setRowStatus(iRow, ClpSimplex::atLowerBound);
                    rowSolution[iRow] = lowerValue;
               } else {
                    model_->setRowStatus(iRow, ClpSimplex::atUpperBound);
                    rowSolution[iRow] = upperValue;
               }
               j = startSet_[iSet];
               while (j >= 0) {
                    DynamicStatus status = getDynamicStatus(j);
                    if (status == inSmall) {
                         int numberThis = startColumn_[j+1] - startColumn_[j] + 1;
                         if (numberElements + numberThis > numberElements_) {
                              // need to redo
                              numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
                              matrix_->reserve(lastDynamic_, numberElements_);
                              element =  matrix_->getMutableElements();
                              row = matrix_->getMutableIndices();
                              // these probably okay but be safe
                              startColumn = matrix_->getMutableVectorStarts();
                              length = matrix_->getMutableVectorLengths();
                         }
                         length[firstAvailable_] = numberThis;
                         cost[firstAvailable_] = cost_[j];
                         CoinBigIndex base = startColumn_[j];
                         for (int k = 0; k < numberThis - 1; k++) {
                              row[numberElements] = row_[base+k];
                              element[numberElements++] = element_[base+k];
                         }
                         row[numberElements] = iRow;
                         element[numberElements++] = 1.0;
                         id_[firstAvailable_-firstDynamic_] = j;
                         solution[firstAvailable_] = 0.0;
                         model_->setStatus(firstAvailable_, ClpSimplex::basic);
                         if (!columnLower_ && !columnUpper_) {
                              columnLower[firstAvailable_] = 0.0;
                              columnUpper[firstAvailable_] = COIN_DBL_MAX;
                         }  else {
                              if (columnLower_)
                                   columnLower[firstAvailable_] = columnLower_[j];
                              else
                                   columnLower[firstAvailable_] = 0.0;
                              if (columnUpper_)
                                   columnUpper[firstAvailable_] = columnUpper_[j];
                              else
                                   columnUpper[firstAvailable_] = COIN_DBL_MAX;
                              if (status != atUpperBound) {
                                   solution[firstAvailable_] = columnLower[firstAvailable_];
                              } else {
                                   solution[firstAvailable_] = columnUpper[firstAvailable_];
                              }
                         }
                         firstAvailable_++;
                         startColumn[firstAvailable_] = numberElements;
                    }
                    j = next_[j]; //onto next in set
               }
               model_->setRowStatus(numberActiveSets_ + numberStaticRows_, getStatus(iSet));
               toIndex_[iSet] = numberActiveSets_;
               fromIndex_[numberActiveSets_++] = iSet;
	  } else {
	    // solo key
	    bool needKey=false;
	    if (numberActive) {
	      if (whichKey<maximumGubColumns_) {
		// structural - assume ok
		needKey = false;
	      } else {
		// slack
		keyVariable_[iSet] = maximumGubColumns_ + iSet;
		double value = keyValue(iSet);
		if (value<lowerSet_[iSet]-1.0e-8||
		    value>upperSet_[iSet]+1.0e-8)
		  needKey=true;
	      }
	    } else {
	      needKey = true;
	    }
	    if (needKey) {
	      // all to lb then up some (slack/null if possible)
	      int length=99999999;
	      int which=-1;
	      double sum=0.0;
	      for (int iColumn=startSet_[iSet];iColumn<startSet_[iSet+1];iColumn++) {
		setDynamicStatus(iColumn,atLowerBound);
		sum += columnLower_[iColumn];
		if (length>startColumn_[iColumn+1]-startColumn_[iColumn]) {
		  which=iColumn;
		  length=startColumn_[iColumn+1]-startColumn_[iColumn];
		}
	      }
	      if (sum>lowerSet_[iSet]-1.0e-8) {
		// slack can be basic
		setStatus(iSet,ClpSimplex::basic);
		keyVariable_[iSet] = maximumGubColumns_ + iSet;
	      } else {
		// use shortest
		setDynamicStatus(which,soloKey);
		keyVariable_[iSet] = which;
		setStatus(iSet,ClpSimplex::atLowerBound);
	      }
	    }
          }
          assert (toIndex_[iSet] >= 0 || whichKey >= 0);
          keyVariable_[iSet] = whichKey;
     }
     // clean up pivotVariable
     int numberColumns = model_->numberColumns();
     int numberRows = model_->numberRows();
     int * pivotVariable = model_->pivotVariable();
     if (pivotVariable) {
       for (int i=0; i<numberStaticRows_+numberActiveSets_;i++) {
	 if (model_->getRowStatus(i)!=ClpSimplex::basic)
	   pivotVariable[i]=-1;
	 else
	   pivotVariable[i]=numberColumns+i;
       }
       for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
	 pivotVariable[i]=i+numberColumns;
       }
       int put=-1;
       for (int i=0;i<numberColumns;i++) {
	 if (model_->getColumnStatus(i)==ClpSimplex::basic) {
	   while(put<numberRows) {
	     put++;
	     if (pivotVariable[put]==-1) {
	       pivotVariable[put]=i;
	       break;
	     }
	   }
	 }
       }
       for (int i=CoinMax(put,0);i<numberRows;i++) {
	 if (pivotVariable[i]==-1) 
	   pivotVariable[i]=i+numberColumns;
       }
     }
     if (rhsOffset_) {
       double * cost = model_->costRegion();
       double * columnLower = model_->lowerRegion();
       double * columnUpper = model_->upperRegion();
       double * solution = model_->solutionRegion();
       int numberRows = model_->numberRows();
       for (int i = numberActiveSets_; i < numberRows-numberStaticRows_; i++) {
	 int iSequence = i + numberStaticRows_ + numberColumns;
	 solution[iSequence] = 0.0;
	 columnLower[iSequence] = -COIN_DBL_MAX;
	 columnUpper[iSequence] = COIN_DBL_MAX;
	 cost[iSequence] = 0.0;
	 model_->nonLinearCost()->setOne(iSequence, solution[iSequence],
					columnLower[iSequence],
					columnUpper[iSequence], 0.0);
	 model_->setStatus(iSequence, ClpSimplex::basic);
	 rhsOffset_[i+numberStaticRows_] = 0.0;
       }
#if 0
       for (int i=0;i<numberStaticRows_;i++)
	 printf("%d offset %g\n",
		i,rhsOffset_[i]);
#endif
     }
     numberActiveColumns_ = firstAvailable_;
#if 0
     for (iSet = 0; iSet < numberSets_; iSet++) {
       for (int j=startSet_[iSet];j<startSet_[iSet+1];j++) {
	 if (getDynamicStatus(j)==soloKey)
	   printf("%d in set %d is solo key\n",j,iSet);
	 else if (getDynamicStatus(j)==inSmall)
	   printf("%d in set %d is in small\n",j,iSet);
       }
     }
#endif
     return;
}
// Writes out model (without names)
void 
ClpDynamicMatrix::writeMps(const char * name)
{
  int numberTotalRows = numberStaticRows_+numberSets_;
  int numberTotalColumns = firstDynamic_+numberGubColumns_;
  // over estimate
  int numberElements = getNumElements()+startColumn_[numberGubColumns_]
    + numberGubColumns_;
  double * columnLower = new double [numberTotalColumns];
  double * columnUpper = new double [numberTotalColumns];
  double * cost = new double [numberTotalColumns];
  double * rowLower = new double [numberTotalRows];
  double * rowUpper = new double [numberTotalRows];
  CoinBigIndex * start = new CoinBigIndex[numberTotalColumns+1];
  int * row = new int [numberElements];
  double * element = new double [numberElements];
  // Fill in
  const CoinBigIndex * startA = getVectorStarts();
  const int * lengthA = getVectorLengths();
  const int * rowA = getIndices();
  const double * elementA = getElements();
  const double * columnLowerA = model_->columnLower();
  const double * columnUpperA = model_->columnUpper();
  const double * costA = model_->objective();
  const double * rowLowerA = model_->rowLower();
  const double * rowUpperA = model_->rowUpper();
  start[0]=0;
  numberElements=0;
  for (int i=0;i<firstDynamic_;i++) {
    columnLower[i] = columnLowerA[i];
    columnUpper[i] = columnUpperA[i];
    cost[i] = costA[i];
    for (CoinBigIndex j = startA[i];j<startA[i]+lengthA[i];j++) {
      row[numberElements] = rowA[j];
      element[numberElements++]=elementA[j];
    }
    start[i+1]=numberElements;
  }
  for (int i=0;i<numberStaticRows_;i++) {
    rowLower[i] = rowLowerA[i];
    rowUpper[i] = rowUpperA[i];
  }
  int putC=firstDynamic_;
  int putR=numberStaticRows_;
  for (int i=0;i<numberSets_;i++) {
    rowLower[putR]=lowerSet_[i];
    rowUpper[putR]=upperSet_[i];
    for (CoinBigIndex k=startSet_[i];k<startSet_[i+1];k++) {
      columnLower[putC]=columnLower_[k];
      columnUpper[putC]=columnUpper_[k];
      cost[putC]=cost_[k];
      putC++;
      for (CoinBigIndex j = startColumn_[k];j<startColumn_[k+1];j++) {
	row[numberElements] = row_[j];
	element[numberElements++]=element_[j];
      }
      row[numberElements] = putR;
      element[numberElements++]=1.0;
      start[putC]=numberElements;
    }
    putR++;
  }

  assert (putR==numberTotalRows);
  assert (putC==numberTotalColumns);
  ClpSimplex modelOut;
  modelOut.loadProblem(numberTotalColumns,numberTotalRows,
		       start,row,element,
		       columnLower,columnUpper,cost,
		       rowLower,rowUpper);
  modelOut.writeMps(name);
  delete [] columnLower;
  delete [] columnUpper;
  delete [] cost;
  delete [] rowLower;
  delete [] rowUpper;
  delete [] start;
  delete [] row; 
  delete [] element;
}
// Adds in a column to gub structure (called from descendant)
int
ClpDynamicMatrix::addColumn(int numberEntries, const int * row, const double * element,
                            double cost, double lower, double upper, int iSet,
                            DynamicStatus status)
{
     // check if already in
     int j = startSet_[iSet];
     while (j >= 0) {
          if (startColumn_[j+1] - startColumn_[j] == numberEntries) {
               const int * row2 = row_ + startColumn_[j];
               const double * element2 = element_ + startColumn_[j];
               bool same = true;
               for (int k = 0; k < numberEntries; k++) {
                    if (row[k] != row2[k] || element[k] != element2[k]) {
                         same = false;
                         break;
                    }
               }
               if (same) {
                    bool odd = false;
                    if (cost != cost_[j])
                         odd = true;
                    if (columnLower_ && lower != columnLower_[j])
                         odd = true;
                    if (columnUpper_ && upper != columnUpper_[j])
                         odd = true;
                    if (odd) {
                         printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n",
                                cost, lower, upper, cost_[j],
                                columnLower_ ? columnLower_[j] : 0.0,
                                columnUpper_ ? columnUpper_[j] : 1.0e100);
                    } else {
		         setDynamicStatus(j, status);
                         return j;
		    }
               }
          }
          j = next_[j];
     }

     if (numberGubColumns_ == maximumGubColumns_ ||
               startColumn_[numberGubColumns_] + numberEntries > maximumElements_) {
          CoinBigIndex j;
          int i;
          int put = 0;
          int numberElements = 0;
          CoinBigIndex start = 0;
          // compress - leave ones at ub and basic
          int * which = new int [numberGubColumns_];
          for (i = 0; i < numberGubColumns_; i++) {
               CoinBigIndex end = startColumn_[i+1];
	       // what about ubs if column generation?
               if (getDynamicStatus(i) != atLowerBound) {
                    // keep in
                    for (j = start; j < end; j++) {
                         row_[numberElements] = row_[j];
                         element_[numberElements++] = element_[j];
                    }
                    startColumn_[put+1] = numberElements;
                    cost_[put] = cost_[i];
                    if (columnLower_)
                         columnLower_[put] = columnLower_[i];
                    if (columnUpper_)
                         columnUpper_[put] = columnUpper_[i];
                    dynamicStatus_[put] = dynamicStatus_[i];
                    id_[put] = id_[i];
                    which[i] = put;
                    put++;
               } else {
                    // out
                    which[i] = -1;
               }
               start = end;
          }
          // now redo startSet_ and next_
          int * newNext = new int [maximumGubColumns_];
          for (int jSet = 0; jSet < numberSets_; jSet++) {
               int sequence = startSet_[jSet];
               while (which[sequence] < 0) {
                    // out
                    assert (next_[sequence] >= 0);
                    sequence = next_[sequence];
               }
               startSet_[jSet] = which[sequence];
               int last = which[sequence];
               while (next_[sequence] >= 0) {
                    sequence = next_[sequence];
                    if(which[sequence] >= 0) {
                         // keep
                         newNext[last] = which[sequence];
                         last = which[sequence];
                    }
               }
               newNext[last] = -jSet - 1;
          }
          delete [] next_;
          next_ = newNext;
          delete [] which;
          abort();
     }
     CoinBigIndex start = startColumn_[numberGubColumns_];
     CoinMemcpyN(row, numberEntries, row_ + start);
     CoinMemcpyN(element, numberEntries, element_ + start);
     startColumn_[numberGubColumns_+1] = start + numberEntries;
     cost_[numberGubColumns_] = cost;
     if (columnLower_)
          columnLower_[numberGubColumns_] = lower;
     else
          assert (!lower);
     if (columnUpper_)
          columnUpper_[numberGubColumns_] = upper;
     else
          assert (upper > 1.0e20);
     setDynamicStatus(numberGubColumns_, status);
     // Do next_
     j = startSet_[iSet];
     startSet_[iSet] = numberGubColumns_;
     next_[numberGubColumns_] = j;
     numberGubColumns_++;
     return numberGubColumns_ - 1;
}
// Returns which set a variable is in
int
ClpDynamicMatrix::whichSet (int sequence) const
{
     while (next_[sequence] >= 0)
          sequence = next_[sequence];
     int iSet = - next_[sequence] - 1;
     return iSet;
}
