// Copyright (C) 2004, International Business Machines
// Corporation and others.  All Rights Reserved.


#include <cstdio>

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

#include "ClpSimplex.hpp"
#include "ClpFactorization.hpp"
#include "ClpQuadraticObjective.hpp"
#include "ClpNonLinearCost.hpp"
// at end to get min/max!
#include "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_,numberSets_);
  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_);
  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_,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_];
  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;
    }
  }
  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;
  // 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_) {
    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 /= (double) numberColumns;
  guess *= 2*numberGubInSmall;
  numberElements_ = (int) guess;
  numberElements_ = CoinMin(numberElements_,numberElements)+originalMatrix->getNumElements();
  matrix_ = originalMatrix;
  flags_ &= ~1;
  // resize model (matrix stays same)
  int newRowSize = numberRows+CoinMin(numberSets_,CoinMax(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,numberSets_);
    assert (dynamicStatus);
    dynamicStatus_ = ClpCopyOfArray(dynamicStatus,numberGubColumns_);
  } else {
    status_= new unsigned char [numberSets_];
    memset(status_,0,numberSets_);
    int i;
    for (i=0;i<numberSets_;i++) {
      // make slack key
      setStatus(i,ClpSimplex::basic);
    }
    dynamicStatus_ = new unsigned char [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_,numberSets_);
    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_);
    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_,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 = (int) (startFraction*numberSets_);
    int endG2 = (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)) {
		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 check)
{
  //forceRefresh=true;
  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]);
      }
      memcpy(saveE,rhs,numberRows*sizeof(double));
      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];
	      } 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];
		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 -(double) lowerSet_[iSet])<
	fabs(valueOut -(double) upperSet_[iSet]))
      setStatus(iSet,ClpSimplex::atLowerBound);
    else
      setStatus(iSet,ClpSimplex::atUpperBound);
    if (lowerSet_[iSet]==upperSet_[iSet])
      setStatus(iSet,ClpSimplex::isFixed);
    if (getStatus(iSet)!=model->getStatus(sequenceOut))
      printf("** set %d status %d, var status %d\n",iSet,
	     getStatus(iSet),model->getStatus(sequenceOut));
  }
  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
  return 0;
}
/*
     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;
  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:
    {
    }
    break;
    // restore status
  case 6:
    {
      // maybe mismatch - redo all in small model
      //for (int i=firstDynamic_;i<firstAvailable_;i++) {
      //int jColumn = id_[i-firstDynamic_];
      //setDynamicStatus(jColumn,inSmall);
      //}
    }
    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_) {
	assert (number==model->sequenceIn());
	// id will be sitting at firstAvailable
	int sequence = id_[firstAvailable_-firstDynamic_];
	setFlagged(sequence);
	model->clearFlagged(firstAvailable_);
      }
    }
  case 11:
    {
      int sequenceIn = model->sequenceIn();
      if (sequenceIn>=firstDynamic_&&sequenceIn<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[sequenceIn]) {
	  int * row = matrix_->getMutableIndices();
	  CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
	  int iRow = row[startColumn[sequenceIn]+length[sequenceIn]-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_[sequenceIn-firstDynamic_];
	if (model->getStatus(sequenceIn)==ClpSimplex::atLowerBound) {
	  setDynamicStatus(iBig,atLowerBound);
	  if (columnLower_)
	    modifyOffset(sequenceIn,columnLower_[iBig]);
	} else {
	  setDynamicStatus(iBig,atUpperBound);
	  modifyOffset(sequenceIn,columnUpper_[iBig]);
	}
      }
    }
    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;
  }
  // 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;
      if (columnUpper_) {
	if (!columnLower_) {
	  if (fabs(value-columnUpper_[jColumn])<fabs(value)) 
	  toLowerBound=false;
	} else if (fabs(value-columnUpper_[jColumn])<fabs(value-columnLower_[iColumn])) {
	  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;
      }
    }
    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++;
    }
  }
  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
      }
    }
  }
  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_;
	if (getStatus(savedBestSet_)==ClpSimplex::atUpperBound)
	  solution[iSequence] = columnUpper[iSequence];
	else
	  solution[iSequence] = columnLower[iSequence];
	// 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]=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;
      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;
	    double distance=0.0;
	    if (iphase||i<numberInSet)
	      assert (solution[i]>=lower[i]&&solution[i]<=upper[i]);
	    if (dj>dualTolerance)
	      improvement = dj*(solution[i]-lower[i]);
	    else if (dj<-dualTolerance)
	      improvement = dj*(solution[i]-upper[i]);
	    if (improvement>best) {
	      best=improvement;
	      chosen=i;
	      if (dj<0.0) {
		way = 1;
		distance = upper[i]-solution[i];
	      } else {
		way = -1;
		distance = solution[i]-lower[i];
	      }
	    }
	  }
	  if (chosen>=0) {
	    improve=true;
	    // now see how far
	    if (way>0) {
	      // incoming increasing so basic decreasing
	      // if phase 0 then go to nearest bound
	      double distance=upper[chosen]-solution[chosen];
	      double basicDistance;
	      if (!iphase) {
		assert (iBasic==numberInSet);
		assert (solution[iBasic]>upper[iBasic]);
		basicDistance = solution[iBasic]-upper[iBasic];
	      } else {
		basicDistance = solution[iBasic]-lower[iBasic];
	      }
	      // need extra coding for unbounded
	      assert (CoinMin(distance,basicDistance)<1.0e20);
	      if (distance>basicDistance) {
		// incoming becomes basic
		solution[chosen] += basicDistance;
		if (!iphase) 
		  solution[iBasic]=upper[iBasic];
		else 
		  solution[iBasic]=lower[iBasic];
		iBasic = chosen;
	      } else {
		// flip
		solution[chosen]=upper[chosen];
		solution[iBasic] -= distance;
	      }
	    } else {
	      // incoming decreasing so basic increasing
	      // if phase 0 then go to nearest bound
	      double distance=solution[chosen]-lower[chosen];
	      double basicDistance;
	      if (!iphase) {
		assert (iBasic==numberInSet);
		assert (solution[iBasic]<lower[iBasic]);
		basicDistance = lower[iBasic]-solution[iBasic];
	      } else {
		basicDistance = upper[iBasic]-solution[iBasic];
	      }
	      // need extra coding for unbounded - for now just exit
	      if (CoinMin(distance,basicDistance)>1.0e20) {
		printf("unbounded on set %d\n",iSet);
		iphase=1;
		iBasic=numberInSet;
		break;
	      }
	      if (distance>basicDistance) {
		// incoming becomes basic
		solution[chosen] -= basicDistance;
		if (!iphase) 
		  solution[iBasic]=lower[iBasic];
		else 
		  solution[iBasic]=upper[iBasic];
		iBasic = chosen;
	      } else {
		// flip
		solution[chosen]=lower[chosen];
		solution[iBasic] += distance;
	      }
	    }
	    if (!iphase) {
	      if(iBasic<numberInSet)
		break; // feasible
	      else if (solution[iBasic]>=lower[iBasic]&&
		       solution[iBasic]<=upper[iBasic])
		break; // feasible (on flip)
	    }
	  }
	}
      }
    }
    // do solution i.e. bounds
    if (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;
    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;
      j = next_[j]; //onto next in set
    }
    if (getStatus(iSet)==ClpSimplex::basic&&numberActive)
      numberActive++;
    if (numberActive) {
      assert (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==atLowerBound) {
	      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;
    }
    assert (toIndex_[iSet]>=0||whichKey>=0);
    keyVariable_[iSet]=whichKey;
  }
  return;
}
// 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
	  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];
      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_];
  memcpy(row_+start,row,numberEntries*sizeof(int));
  memcpy(element_+start,element,numberEntries*sizeof(float));
  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;
}