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

#include "CoinPragma.hpp"
#include "ClpFactorization.hpp"
#ifndef SLIM_CLP
#include "ClpQuadraticObjective.hpp"
#endif
#include "CoinHelperFunctions.hpp"
#include "CoinIndexedVector.hpp"
#include "ClpSimplex.hpp"
#include "ClpMatrixBase.hpp"
#ifndef SLIM_CLP
#include "ClpNetworkBasis.hpp"
#include "ClpNetworkMatrix.hpp"
//#define CHECK_NETWORK
#ifdef CHECK_NETWORK
const static bool doCheck=true;
#else
const static bool doCheck=false;
#endif
#endif
//#define CLP_FACTORIZATION_INSTRUMENT
#ifdef CLP_FACTORIZATION_INSTRUMENT
#include "CoinTime.hpp"
double factorization_instrument(int type)
{
  static int times[10]={0,0,0,0,0,0,0,0,0,0};
  static double startTime=0.0;
  static double totalTimes [10]={0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
  if (type<0) {
    assert (!startTime);
    startTime=CoinCpuTime();
    return 0.0;
  } else if (type>0) {
    times[type]++;
    double difference = CoinCpuTime()-startTime;
    totalTimes[type] += difference;
    startTime=0.0;
    return difference;
  } else {
    // report
    const char * types[10]=
      {"","fac=rhs_etc","factorize","replace","update_FT",
       "update","update_transpose","gosparse","getWeights!","update2_FT"};
    double total=0.0;
    for (int i=1;i<10;i++) {
      if (times[i]) {
	printf("%s was called %d times taking %g seconds\n",
	       types[i],times[i],totalTimes[i]);
	total += totalTimes[i];
	times[i]=0;
	totalTimes[i]=0.0;
      }
    }
    return total;
  }
}
#endif
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
#ifndef CLP_MULTIPLE_FACTORIZATIONS    

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
ClpFactorization::ClpFactorization () :
   CoinFactorization() 
{
#ifndef SLIM_CLP
  networkBasis_ = NULL;
#endif
}

//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
ClpFactorization::ClpFactorization (const ClpFactorization & rhs) :
   CoinFactorization(rhs) 
{
#ifndef SLIM_CLP
  if (rhs.networkBasis_)
    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
  else
    networkBasis_=NULL;
#endif
}

ClpFactorization::ClpFactorization (const CoinFactorization & rhs) :
   CoinFactorization(rhs) 
{
#ifndef SLIM_CLP
  networkBasis_=NULL;
#endif
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
ClpFactorization::~ClpFactorization () 
{
#ifndef SLIM_CLP
  delete networkBasis_;
#endif
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpFactorization &
ClpFactorization::operator=(const ClpFactorization& rhs)
{
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(-1);
#endif
  if (this != &rhs) {
    CoinFactorization::operator=(rhs);
#ifndef SLIM_CLP
    delete networkBasis_;
    if (rhs.networkBasis_)
      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
    else
      networkBasis_=NULL;
#endif
  }
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(1);
#endif
  return *this;
}
#if 0
static unsigned int saveList[10000];
int numberSave=-1;
inline bool isDense(int i) {
  return ((saveList[i>>5]>>(i&31))&1)!=0;
}
inline void setDense(int i) {
  unsigned int & value = saveList[i>>5];
  int bit = i&31;
  value |= (1<<bit);
}
#endif
int 
ClpFactorization::factorize ( ClpSimplex * model,
			      int solveType, bool valuesPass)
{
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(-1);
#endif
  ClpMatrixBase * matrix = model->clpMatrix(); 
  int numberRows = model->numberRows();
  int numberColumns = model->numberColumns();
  if (!numberRows)
    return 0;
  // If too many compressions increase area
  if (numberPivots_>1&&numberCompressions_*10 > numberPivots_+10) {
    areaFactor_ *= 1.1;
  }
  //int numberPivots=numberPivots_;
#if 0
  if (model->algorithm()>0)
    numberSave=-1;
#endif
#ifndef SLIM_CLP
  if (!networkBasis_||doCheck) {
#endif
    status_=-99;
    int * pivotVariable = model->pivotVariable();
    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
    while (status_<-98) {
      
      int i;
      int numberBasic=0;
      int numberRowBasic;
      // Move pivot variables across if they look good
      int * pivotTemp = model->rowArray(0)->getIndices();
      assert (!model->rowArray(0)->getNumElements());
      if (!matrix->rhsOffset(model)) {
#if 0
        if (numberSave>0) {
          int nStill=0;
          int nAtBound=0;
          int nZeroDual=0;
          CoinIndexedVector * array = model->rowArray(3);
          CoinIndexedVector * objArray = model->columnArray(1);
          array->clear();
          objArray->clear();
          double * cost = model->costRegion();
          double tolerance = model->primalTolerance();
          double offset=0.0;
          for (i=0;i<numberRows;i++) {
            int iPivot = pivotVariable[i];
            if (iPivot<numberColumns&&isDense(iPivot)) {
              if (model->getColumnStatus(iPivot)==ClpSimplex::basic) {
                nStill++;
                double value=model->solutionRegion()[iPivot];
                double dual = model->dualRowSolution()[i];
                double lower=model->lowerRegion()[iPivot];
                double upper=model->upperRegion()[iPivot];
                ClpSimplex::Status status;
                if (fabs(value-lower)<tolerance) {
                  status=ClpSimplex::atLowerBound;
                  nAtBound++;
                } else if (fabs(value-upper)<tolerance) {
                  nAtBound++;
                  status=ClpSimplex::atUpperBound;
                } else if (value>lower&&value<upper) {
                  status=ClpSimplex::superBasic;
                } else {
                  status=ClpSimplex::basic;
                }
                if (status!=ClpSimplex::basic) {
                  if (model->getRowStatus(i)!=ClpSimplex::basic) {
                    model->setColumnStatus(iPivot,ClpSimplex::atLowerBound);
                    model->setRowStatus(i,ClpSimplex::basic);
                    pivotVariable[i]=i+numberColumns;
                    model->dualRowSolution()[i]=0.0;
                    model->djRegion(0)[i]=0.0;
                    array->add(i,dual);
                    offset += dual*model->solutionRegion(0)[i];
                  }
                }
                if (fabs(dual)<1.0e-5)
                  nZeroDual++;
              }
            }
          }
          printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
                 numberSave,nStill,nAtBound,nZeroDual,offset);
          if (array->getNumElements()) {
            // modify costs
            model->clpMatrix()->transposeTimes(model,1.0,array,model->columnArray(0),
                                               objArray);
            array->clear();
            int n=objArray->getNumElements();
            int * indices = objArray->getIndices();
            double * elements = objArray->denseVector();
            for (i=0;i<n;i++) {
              int iColumn = indices[i];
              cost[iColumn] -= elements[iColumn];
              elements[iColumn]=0.0;
            }
            objArray->setNumElements(0);
          }
        }
#endif
	// Seems to prefer things in order so quickest
	// way is to go though like this
	for (i=0;i<numberRows;i++) {
	  if (model->getRowStatus(i) == ClpSimplex::basic) 
	    pivotTemp[numberBasic++]=i;
	}
	numberRowBasic=numberBasic;
	/* Put column basic variables into pivotVariable
	   This is done by ClpMatrixBase to allow override for gub
	*/
	matrix->generalExpanded(model,0,numberBasic);
      } else {
	// Long matrix - do a different way
	bool fullSearch=false;
	for (i=0;i<numberRows;i++) {
	  int iPivot = pivotVariable[i];
	  if (iPivot>=numberColumns) {
	    pivotTemp[numberBasic++]=iPivot-numberColumns;
	  }
	}
	numberRowBasic=numberBasic;
	for (i=0;i<numberRows;i++) {
	  int iPivot = pivotVariable[i];
	  if (iPivot<numberColumns) {
	    if (iPivot>=0) {
	      pivotTemp[numberBasic++]=iPivot;
	    } else {
	      // not full basis
	      fullSearch=true;
	      break;
	    }
	  }
	}
	if (fullSearch) {
	  // do slow way
	  numberBasic=0;
	  for (i=0;i<numberRows;i++) {
	    if (model->getRowStatus(i) == ClpSimplex::basic) 
	      pivotTemp[numberBasic++]=i;
	  }
	  numberRowBasic=numberBasic;
	  /* Put column basic variables into pivotVariable
	     This is done by ClpMatrixBase to allow override for gub
	  */
	  matrix->generalExpanded(model,0,numberBasic);
	}
      }
      if (numberBasic>model->maximumBasic()) {
#if 0 // ndef NDEBUG
        printf("%d basic - should only be %d\n",
               numberBasic,numberRows);
#endif
        // Take out some
        numberBasic=numberRowBasic;
        for (int i=0;i<numberColumns;i++) {
          if (model->getColumnStatus(i) == ClpSimplex::basic) {
            if (numberBasic<numberRows)
              numberBasic++;
            else
              model->setColumnStatus(i,ClpSimplex::superBasic);
          }
        }
        numberBasic=numberRowBasic;
        matrix->generalExpanded(model,0,numberBasic);
      }
#ifndef SLIM_CLP
      // see if matrix a network
#ifndef NO_RTTI
      ClpNetworkMatrix* networkMatrix =
	dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
#else
      ClpNetworkMatrix* networkMatrix = NULL;
      if (model->clpMatrix()->type()==11)
	networkMatrix = 
	static_cast< ClpNetworkMatrix*>(model->clpMatrix());
#endif
      // If network - still allow ordinary factorization first time for laziness
      if (networkMatrix)
	biasLU_=0; // All to U if network
      //int saveMaximumPivots = maximumPivots();
      delete networkBasis_;
      networkBasis_ = NULL;
      if (networkMatrix&&!doCheck)
	maximumPivots(1);
#endif
      //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
      while (status_==-99) {
	// maybe for speed will be better to leave as many regions as possible
	gutsOfDestructor();
	gutsOfInitialize(2);
	CoinBigIndex numberElements=numberRowBasic;

	// compute how much in basis

	int i;
	// can change for gub
	int numberColumnBasic = numberBasic-numberRowBasic;

	numberElements +=matrix->countBasis(model,
					   pivotTemp+numberRowBasic, 
					   numberRowBasic,
					    numberColumnBasic);
	// and recompute as network side say different
	if (model->numberIterations())
	  numberRowBasic = numberBasic - numberColumnBasic;
	numberElements = 3 * numberBasic + 3 * numberElements + 20000;
#if 0
	// If iteration not zero then may be compressed
	getAreas ( !model->numberIterations() ? numberRows : numberBasic, 
		   numberRowBasic+numberColumnBasic, numberElements,
		   2 * numberElements );
#else
	getAreas ( numberRows,
		   numberRowBasic+numberColumnBasic, numberElements,
		   2 * numberElements );
#endif
	//fill
	// Fill in counts so we can skip part of preProcess
	int * numberInRow = numberInRow_.array();
	int * numberInColumn = numberInColumn_.array();
	CoinZeroN ( numberInRow, numberRows_ + 1 );
	CoinZeroN ( numberInColumn, maximumColumnsExtra_ + 1 );
	double * elementU = elementU_.array();
	int * indexRowU = indexRowU_.array();
	CoinBigIndex * startColumnU = startColumnU_.array();
	for (i=0;i<numberRowBasic;i++) {
	  int iRow = pivotTemp[i];
	  indexRowU[i]=iRow;
	  startColumnU[i]=i;
	  elementU[i]=slackValue_;
	  numberInRow[iRow]=1;
	  numberInColumn[i]=1;
	}
	startColumnU[numberRowBasic]=numberRowBasic;
	// can change for gub so redo
	numberColumnBasic = numberBasic-numberRowBasic;
	matrix->fillBasis(model, 
			  pivotTemp+numberRowBasic, 
			  numberColumnBasic,
			  indexRowU, 
			  startColumnU+numberRowBasic,
			  numberInRow,
			  numberInColumn+numberRowBasic,
			  elementU);
#if 0
	{
	  printf("%d row basic, %d column basic\n",numberRowBasic,numberColumnBasic);
	  for (int i=0;i<numberElements;i++) 
	    printf("row %d col %d value %g\n",indexRowU_.array()[i],indexColumnU_[i],
		   elementU_.array()[i]);
	}
#endif
	// recompute number basic
        numberBasic = numberRowBasic+numberColumnBasic;
	if (numberBasic) 
	  numberElements = startColumnU[numberBasic-1]
	    +numberInColumn[numberBasic-1];
	else
	  numberElements=0;
	lengthU_ = numberElements;
        //saveFactorization("dump.d");
	if (biasLU_>=3||numberRows_!=numberColumns_)
	  preProcess ( 2 );
	else
	  preProcess ( 3 ); // no row copy
	factor (  );
	if (status_==-99) {
	  // get more memory
	  areaFactor(2.0*areaFactor());
	} else if (status_==-1&&model->numberIterations()==0&&
                   denseThreshold_) {
          // Round again without dense
          denseThreshold_=0;
          status_=-99;
        }
      }
      // If we get here status is 0 or -1
      if (status_ == 0) {
	// We may need to tamper with order and redo - e.g. network with side
	int useNumberRows = numberRows;
	// **** we will also need to add test in dual steepest to do
	// as we do for network
        matrix->generalExpanded(model,12,useNumberRows);
	const int * permuteBack = permuteBack_.array();
	const int * back = pivotColumnBack_.array();
	//int * pivotTemp = pivotColumn_.array();
	//ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
	// Redo pivot order
	for (i=0;i<numberRowBasic;i++) {
	  int k = pivotTemp[i];
	  // so rowIsBasic[k] would be permuteBack[back[i]]
	  pivotVariable[permuteBack[back[i]]]=k+numberColumns;
	}
	for (;i<useNumberRows;i++) {
	  int k = pivotTemp[i];
	  // so rowIsBasic[k] would be permuteBack[back[i]]
	  pivotVariable[permuteBack[back[i]]]=k;
	}
#if 0
        if (numberSave>=0) {
          numberSave=numberDense_;
          memset(saveList,0,((numberRows_+31)>>5)*sizeof(int));
          for (i=numberRows_-numberSave;i<numberRows_;i++) {
            int k=pivotTemp[pivotColumn_.array()[i]];
            setDense(k);
          }
        }
#endif
	// Set up permutation vector
	// these arrays start off as copies of permute
	// (and we could use permute_ instead of pivotColumn (not back though))
	ClpDisjointCopyN ( permute_.array(), useNumberRows , pivotColumn_.array()  );
	ClpDisjointCopyN ( permuteBack_.array(), useNumberRows , pivotColumnBack_.array()  );
#ifndef SLIM_CLP
	if (networkMatrix) {
	  maximumPivots(CoinMax(2000,maximumPivots()));
	  // redo arrays
	  for (int iRow=0;iRow<4;iRow++) {
	    int length =model->numberRows()+maximumPivots();
	    if (iRow==3||model->objectiveAsObject()->type()>1)
	      length += model->numberColumns();
	    model->rowArray(iRow)->reserve(length);
	  }
	  // create network factorization
	  if (doCheck)
	    delete networkBasis_; // temp
	  networkBasis_ = new ClpNetworkBasis(model,numberRows_,
					      pivotRegion_.array(),
					      permuteBack_.array(),
					      startColumnU_.array(),
					      numberInColumn_.array(),
					      indexRowU_.array(),
					      elementU_.array());
	  // kill off arrays in ordinary factorization
	  if (!doCheck) {
	    gutsOfDestructor();
	    // but make sure numberRows_ set
	    numberRows_ = model->numberRows();
	    status_=0;
#if 0
	    // but put back permute arrays so odd things will work
	    int numberRows = model->numberRows();
	    pivotColumnBack_ = new int [numberRows];
	    permute_ = new int [numberRows];
	    int i;
	    for (i=0;i<numberRows;i++) {
	      pivotColumnBack_[i]=i;
	      permute_[i]=i;
	    }
#endif
	  }
	} else {
#endif
	  // See if worth going sparse and when
	  if (numberFtranCounts_>100) {
	    ftranCountInput_= CoinMax(ftranCountInput_,1.0);
	    ftranAverageAfterL_ = CoinMax(ftranCountAfterL_/ftranCountInput_,1.0);
	    ftranAverageAfterR_ = CoinMax(ftranCountAfterR_/ftranCountAfterL_,1.0);
	    ftranAverageAfterU_ = CoinMax(ftranCountAfterU_/ftranCountAfterR_,1.0);
	    if (btranCountInput_&&btranCountAfterU_&&btranCountAfterR_) {
	      btranAverageAfterU_ = CoinMax(btranCountAfterU_/btranCountInput_,1.0);
	      btranAverageAfterR_ = CoinMax(btranCountAfterR_/btranCountAfterU_,1.0);
	      btranAverageAfterL_ = CoinMax(btranCountAfterL_/btranCountAfterR_,1.0);
	    } else {
	      // we have not done any useful btrans (values pass?)
	      btranAverageAfterU_ = 1.0;
	      btranAverageAfterR_ = 1.0;
	      btranAverageAfterL_ = 1.0;
	    }
	  }
	  // scale back
	  
	  ftranCountInput_ *= 0.8;
	  ftranCountAfterL_ *= 0.8;
	  ftranCountAfterR_ *= 0.8;
	  ftranCountAfterU_ *= 0.8;
	  btranCountInput_ *= 0.8;
	  btranCountAfterU_ *= 0.8;
	  btranCountAfterR_ *= 0.8;
	  btranCountAfterL_ *= 0.8;
#ifndef SLIM_CLP
	}
#endif
      } else if (status_ == -1&&(solveType==0||solveType==2)) {
	// This needs redoing as it was merged coding - does not need array
	int numberTotal = numberRows+numberColumns;
	int * isBasic = new int [numberTotal];
	int * rowIsBasic = isBasic+numberColumns;
	int * columnIsBasic = isBasic;
	for (i=0;i<numberTotal;i++) 
	  isBasic[i]=-1;
	for (i=0;i<numberRowBasic;i++) {
	  int iRow = pivotTemp[i];
	  rowIsBasic[iRow]=1;
	}
	for (;i<numberBasic;i++) {
	  int iColumn = pivotTemp[i];
	  columnIsBasic[iColumn]=1;
	}
	numberBasic=0;
	for (i=0;i<numberRows;i++) 
	  pivotVariable[i]=-1;
	// mark as basic or non basic
	const int * pivotColumn = pivotColumn_.array();
	for (i=0;i<numberRows;i++) {
	  if (rowIsBasic[i]>=0) {
	    if (pivotColumn[numberBasic]>=0) {
	      rowIsBasic[i]=pivotColumn[numberBasic];
	    } else {
	      rowIsBasic[i]=-1;
              model->setRowStatus(i,ClpSimplex::superBasic);
            }
	    numberBasic++;
	  }
	}
	for (i=0;i<numberColumns;i++) {
	  if (columnIsBasic[i]>=0) {
	    if (pivotColumn[numberBasic]>=0) 
	      columnIsBasic[i]=pivotColumn[numberBasic];
	    else
	      columnIsBasic[i]=-1;
	    numberBasic++;
	  }
	}
	// leave pivotVariable in useful form for cleaning basis
	int * pivotVariable = model->pivotVariable();
	for (i=0;i<numberRows;i++) {
	  pivotVariable[i]=-1;
	}

	for (i=0;i<numberRows;i++) {
	  if (model->getRowStatus(i) == ClpSimplex::basic) {
	    int iPivot = rowIsBasic[i];
	    if (iPivot>=0) 
	      pivotVariable[iPivot]=i+numberColumns;
	  }
	}
	for (i=0;i<numberColumns;i++) {
	  if (model->getColumnStatus(i) == ClpSimplex::basic) {
	    int iPivot = columnIsBasic[i];
	    if (iPivot>=0) 
	      pivotVariable[iPivot]=i;
	  }
	}
	delete [] isBasic;
	double * columnLower = model->lowerRegion();
	double * columnUpper = model->upperRegion();
	double * columnActivity = model->solutionRegion();
	double * rowLower = model->lowerRegion(0);
	double * rowUpper = model->upperRegion(0);
	double * rowActivity = model->solutionRegion(0);
	//redo basis - first take ALL columns out
	int iColumn;
	double largeValue = model->largeValue();
	for (iColumn=0;iColumn<numberColumns;iColumn++) {
	  if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
	    // take out
	    if (!valuesPass) {
	      double lower=columnLower[iColumn];
	      double upper=columnUpper[iColumn];
	      double value=columnActivity[iColumn];
	      if (lower>-largeValue||upper<largeValue) {
		if (fabs(value-lower)<fabs(value-upper)) {
		  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
		  columnActivity[iColumn]=lower;
		} else {
		  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
		  columnActivity[iColumn]=upper;
		}
	      } else {
		model->setColumnStatus(iColumn,ClpSimplex::isFree);
	      }
	    } else {
	      model->setColumnStatus(iColumn,ClpSimplex::superBasic);
	    }
	  }
	}
	int iRow;
	for (iRow=0;iRow<numberRows;iRow++) {
	  int iSequence=pivotVariable[iRow];
	  if (iSequence>=0) {
	    // basic
	    if (iSequence>=numberColumns) {
	      // slack in - leave
	      //assert (iSequence-numberColumns==iRow);
	    } else {
              assert(model->getRowStatus(iRow)!=ClpSimplex::basic);
	      // put back structural
	      model->setColumnStatus(iSequence,ClpSimplex::basic);
	    }
	  } else {
	    // put in slack
	    model->setRowStatus(iRow,ClpSimplex::basic);
	  }
	}
	// Put back any key variables for gub (status_ not touched)
	matrix->generalExpanded(model,1,status_);
	// signal repeat
	status_=-99;
	// set fixed if they are
	for (iRow=0;iRow<numberRows;iRow++) {
	  if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
	    if (rowLower[iRow]==rowUpper[iRow]) {
	      rowActivity[iRow]=rowLower[iRow];
	      model->setRowStatus(iRow,ClpSimplex::isFixed);
	    }
	  }
	}
	for (iColumn=0;iColumn<numberColumns;iColumn++) {
	  if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
	    if (columnLower[iColumn]==columnUpper[iColumn]) {
	      columnActivity[iColumn]=columnLower[iColumn];
	      model->setColumnStatus(iColumn,ClpSimplex::isFixed);
	    }
	  }
	}
      } 
    }
#ifndef SLIM_CLP
  } else {
    // network - fake factorization - do nothing
    status_=0;
    numberPivots_ = 0;
  }
#endif
#ifndef SLIM_CLP
  if (!status_) {
    // take out part if quadratic
    if (model->algorithm()==2) {
      ClpObjective * obj = model->objectiveAsObject();
      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
      assert (quadraticObj);
      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
      int numberXColumns = quadratic->getNumCols();
      assert (numberXColumns<numberColumns);
      int base = numberColumns-numberXColumns;
      int * which = new int [numberXColumns];
      int * pivotVariable = model->pivotVariable();
      int * permute = pivotColumn();
      int i;
      int n=0;
      for (i=0;i<numberRows;i++) {
	int iSj = pivotVariable[i]-base;
	if (iSj>=0&&iSj<numberXColumns) 
	  which[n++]=permute[i];
      }
      if (n)
	emptyRows(n,which);
      delete [] which;
    }
  }
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(2);
#endif
  return status_;
}
/* Replaces one Column to basis,
   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
   If checkBeforeModifying is true will do all accuracy checks
   before modifying factorization.  Whether to set this depends on
   speed considerations.  You could just do this on first iteration
   after factorization and thereafter re-factorize
   partial update already in U */
int 
ClpFactorization::replaceColumn ( const ClpSimplex * model, 
				  CoinIndexedVector * regionSparse,
				  CoinIndexedVector * tableauColumn,
				  int pivotRow,
				  double pivotCheck ,
				  bool checkBeforeModifying)
{
  int returnCode;
#ifndef SLIM_CLP
  if (!networkBasis_) {
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(-1);
#endif
    // see if FT
    if (doForrestTomlin_) {
      returnCode= CoinFactorization::replaceColumn(regionSparse,
					      pivotRow,
					      pivotCheck,
					      checkBeforeModifying);
    } else {
      returnCode= CoinFactorization::replaceColumnPFI(tableauColumn,
					      pivotRow,pivotCheck); // Note array
    }
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(3);
#endif

#ifndef SLIM_CLP
  } else {
    if (doCheck) {
      returnCode = CoinFactorization::replaceColumn(regionSparse,
							pivotRow,
							pivotCheck,
							checkBeforeModifying);
      networkBasis_->replaceColumn(regionSparse,
				   pivotRow);
    } else {
      // increase number of pivots
      numberPivots_++;
      returnCode = networkBasis_->replaceColumn(regionSparse,
				   pivotRow);
    }
  }
#endif
  return returnCode;
}

/* Updates one column (FTRAN) from region2
   number returned is negative if no room
   region1 starts as zero and is zero at end */
int 
ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
				   CoinIndexedVector * regionSparse2)
{
#ifdef CLP_DEBUG
  regionSparse->checkClear();
#endif
  if (!numberRows_)
    return 0;
#ifndef SLIM_CLP
  if (!networkBasis_) {
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(-1);
#endif
    collectStatistics_ = true;
    int returnCode= CoinFactorization::updateColumnFT(regionSparse,
						       regionSparse2);
    collectStatistics_ = false;
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(4);
#endif
    return returnCode;
#ifndef SLIM_CLP
  } else {
#ifdef CHECK_NETWORK
    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
    double * check = new double[numberRows_];
    int returnCode = CoinFactorization::updateColumnFT(regionSparse,
						       regionSparse2);
    networkBasis_->updateColumn(regionSparse,save,-1);
    int i;
    double * array = regionSparse2->denseVector();
    int * indices = regionSparse2->getIndices();
    int n=regionSparse2->getNumElements();
    memset(check,0,numberRows_*sizeof(double));
    double * array2 = save->denseVector();
    int * indices2 = save->getIndices();
    int n2=save->getNumElements();
    assert (n==n2);
    if (save->packedMode()) {
      for (i=0;i<n;i++) {
	check[indices[i]]=array[i];
      }
      for (i=0;i<n;i++) {
	double value2 = array2[i];
	assert (check[indices2[i]]==value2);
      }
    } else {
      for (i=0;i<numberRows_;i++) {
	double value1 = array[i];
	double value2 = array2[i];
	assert (value1==value2);
      }
    }
    delete save;
    delete [] check;
    return returnCode;
#else
    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
    return 1;
#endif
  }
#endif
}
/* Updates one column (FTRAN) from region2
   number returned is negative if no room
   region1 starts as zero and is zero at end */
int 
ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
				 CoinIndexedVector * regionSparse2,
				 bool noPermute) const
{
#ifdef CLP_DEBUG
  if (!noPermute)
    regionSparse->checkClear();
#endif
  if (!numberRows_)
    return 0;
#ifndef SLIM_CLP
  if (!networkBasis_) {
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(-1);
#endif
    collectStatistics_ = true;
    int returnCode= CoinFactorization::updateColumn(regionSparse,
						     regionSparse2,
						     noPermute);
    collectStatistics_ = false;
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(5);
#endif
    return returnCode;
#ifndef SLIM_CLP
  } else {
#ifdef CHECK_NETWORK
    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
    double * check = new double[numberRows_];
    int returnCode = CoinFactorization::updateColumn(regionSparse,
						     regionSparse2,
						     noPermute);
    networkBasis_->updateColumn(regionSparse,save,-1);
    int i;
    double * array = regionSparse2->denseVector();
    int * indices = regionSparse2->getIndices();
    int n=regionSparse2->getNumElements();
    memset(check,0,numberRows_*sizeof(double));
    double * array2 = save->denseVector();
    int * indices2 = save->getIndices();
    int n2=save->getNumElements();
    assert (n==n2);
    if (save->packedMode()) {
      for (i=0;i<n;i++) {
	check[indices[i]]=array[i];
      }
      for (i=0;i<n;i++) {
	double value2 = array2[i];
	assert (check[indices2[i]]==value2);
      }
    } else {
      for (i=0;i<numberRows_;i++) {
	double value1 = array[i];
	double value2 = array2[i];
	assert (value1==value2);
      }
    }
    delete save;
    delete [] check;
    return returnCode;
#else
    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
    return 1;
#endif
  }
#endif
}
/* Updates one column (FTRAN) from region2
   Tries to do FT update
   number returned is negative if no room.
   Also updates region3
   region1 starts as zero and is zero at end */
int 
ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
				       CoinIndexedVector * regionSparse2,
				       CoinIndexedVector * regionSparse3,
				       bool noPermuteRegion3) 
{
  int returnCode=updateColumnFT(regionSparse1,regionSparse2);
  updateColumn(regionSparse1,regionSparse3,noPermuteRegion3);
  return returnCode;
}
/* Updates one column (FTRAN) from region2
   number returned is negative if no room
   region1 starts as zero and is zero at end */
int 
ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
				 CoinIndexedVector * regionSparse2,
				 bool noPermute) const
{
  if (!noPermute)
    regionSparse->checkClear();
  if (!numberRows_)
    return 0;
  collectStatistics_ = false;
  int returnCode= CoinFactorization::updateColumn(regionSparse,
                                                   regionSparse2,
                                                   noPermute);
  return returnCode;
}
/* Updates one column (BTRAN) from region2
   region1 starts as zero and is zero at end */
int 
ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
    			  CoinIndexedVector * regionSparse2) const
{
  if (!numberRows_)
    return 0;
#ifndef SLIM_CLP
  if (!networkBasis_) {
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(-1);
#endif
    collectStatistics_ = true;
    int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
						    regionSparse2);
    collectStatistics_ = false;
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(6);
#endif
    return returnCode;
#ifndef SLIM_CLP
  } else {
#ifdef CHECK_NETWORK
    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
    double * check = new double[numberRows_];
    int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
							      regionSparse2);
    networkBasis_->updateColumnTranspose(regionSparse,save);
    int i;
    double * array = regionSparse2->denseVector();
    int * indices = regionSparse2->getIndices();
    int n=regionSparse2->getNumElements();
    memset(check,0,numberRows_*sizeof(double));
    double * array2 = save->denseVector();
    int * indices2 = save->getIndices();
    int n2=save->getNumElements();
    assert (n==n2);
    if (save->packedMode()) {
      for (i=0;i<n;i++) {
	check[indices[i]]=array[i];
      }
      for (i=0;i<n;i++) {
	double value2 = array2[i];
	assert (check[indices2[i]]==value2);
      }
    } else {
      for (i=0;i<numberRows_;i++) {
	double value1 = array[i];
	double value2 = array2[i];
	assert (value1==value2);
      }
    }
    delete save;
    delete [] check;
    return returnCode;
#else
    return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
#endif
  }
#endif
}
/* makes a row copy of L for speed and to allow very sparse problems */
void 
ClpFactorization::goSparse()
{
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(-1);
#endif
#ifndef SLIM_CLP
  if (!networkBasis_) 
#endif
    CoinFactorization::goSparse();
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(7);
#endif
}
// Cleans up i.e. gets rid of network basis 
void 
ClpFactorization::cleanUp()
{
#ifndef SLIM_CLP
  delete networkBasis_;
  networkBasis_=NULL;
#endif
  resetStatistics();
}
/// Says whether to redo pivot order
bool 
ClpFactorization::needToReorder() const
{
#ifdef CHECK_NETWORK
  return true;
#endif
#ifndef SLIM_CLP
  if (!networkBasis_)
#endif
    return true;
#ifndef SLIM_CLP
  else
    return false;
#endif
}
// Get weighted row list 
void
ClpFactorization::getWeights(int * weights) const
{
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(-1);
#endif
#ifndef SLIM_CLP
  if (networkBasis_) {
    // Network - just unit
    for (int i=0;i<numberRows_;i++) 
      weights[i]=1;
    return;
  }
#endif
  int * numberInRow = numberInRow_.array();
  int * numberInColumn = numberInColumn_.array();
  int * permuteBack = pivotColumnBack_.array();
  int * indexRowU = indexRowU_.array();
  const CoinBigIndex * startColumnU = startColumnU_.array();
  const CoinBigIndex * startRowL = startRowL_.array();
  if (!startRowL||!numberInRow_.array()) {
    int * temp = new int[numberRows_];
    memset(temp,0,numberRows_*sizeof(int));
    int i;
    for (i=0;i<numberRows_;i++) {
      // one for pivot
      temp[i]++;
      CoinBigIndex j;
      for (j=startColumnU[i];j<startColumnU[i]+numberInColumn[i];j++) {
	int iRow=indexRowU[j];
	temp[iRow]++;
      }
    }
    CoinBigIndex * startColumnL = startColumnL_.array();
    int * indexRowL = indexRowL_.array();
    for (i=baseL_;i<baseL_+numberL_;i++) {
      CoinBigIndex j;
      for (j=startColumnL[i];j<startColumnL[i+1];j++) {
	int iRow = indexRowL[j];
	temp[iRow]++;
      }
    }
    for (i=0;i<numberRows_;i++) {
      int number = temp[i];
      int iPermute = permuteBack[i];
      weights[iPermute]=number;
    }
    delete [] temp;
  } else {
    int i;
    for (i=0;i<numberRows_;i++) {
      int number = startRowL[i+1]-startRowL[i]+numberInRow[i]+1;
      //number = startRowL[i+1]-startRowL[i]+1;
      //number = numberInRow[i]+1;
      int iPermute = permuteBack[i];
      weights[iPermute]=number;
    }
  }
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(8);
#endif
}
#else
// This one allows multiple factorizations

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
ClpFactorization::ClpFactorization () 
{
#ifndef SLIM_CLP
  networkBasis_ = NULL;
#endif
  //coinFactorizationA_ = NULL;
  coinFactorizationA_ = new CoinFactorization() ;
  coinFactorizationB_ = NULL;
  //coinFactorizationB_ = new CoinDenseFactorization();
  goDenseThreshold_ = -1;
}

//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
ClpFactorization::ClpFactorization (const ClpFactorization & rhs) 
{
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(-1);
#endif
#ifndef SLIM_CLP
  if (rhs.networkBasis_)
    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
  else
    networkBasis_=NULL;
#endif
  if (rhs.coinFactorizationA_)
    coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_));
  else
    coinFactorizationA_=NULL;
  if (rhs.coinFactorizationB_)
    coinFactorizationB_ = new CoinDenseFactorization(*(rhs.coinFactorizationB_));
  else
    coinFactorizationB_=NULL;
  goDenseThreshold_ = rhs.goDenseThreshold_;
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(1);
#endif
}

ClpFactorization::ClpFactorization (const CoinFactorization & rhs) 
{
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(-1);
#endif
#ifndef SLIM_CLP
  networkBasis_=NULL;
#endif
  coinFactorizationA_ = new CoinFactorization(rhs);
  coinFactorizationB_=NULL;
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(1);
  goDenseThreshold_ = -1;
#endif
}

ClpFactorization::ClpFactorization (const CoinDenseFactorization & rhs) 
{
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(-1);
#endif
#ifndef SLIM_CLP
  networkBasis_=NULL;
#endif
  coinFactorizationA_ = NULL;
  coinFactorizationB_ = new CoinDenseFactorization(rhs);
  goDenseThreshold_ = -1;
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(1);
#endif
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
ClpFactorization::~ClpFactorization () 
{
#ifndef SLIM_CLP
  delete networkBasis_;
#endif
  delete coinFactorizationA_;
  delete coinFactorizationB_;
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpFactorization &
ClpFactorization::operator=(const ClpFactorization& rhs)
{
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(-1);
#endif
  if (this != &rhs) {
#ifndef SLIM_CLP
    delete networkBasis_;
    if (rhs.networkBasis_)
      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
    else
      networkBasis_=NULL;
#endif
    delete coinFactorizationA_;
    if (rhs.coinFactorizationA_)
      coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_));
    else
      coinFactorizationA_=NULL;
    delete coinFactorizationB_;
    if (rhs.coinFactorizationB_)
      coinFactorizationB_ = new CoinDenseFactorization(*(rhs.coinFactorizationB_));
    else
      coinFactorizationB_=NULL;
  }
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(1);
#endif
  return *this;
}
// Go over to dense code
void 
ClpFactorization::goDense() 
{
  delete coinFactorizationA_;
  delete coinFactorizationB_;
  coinFactorizationA_ = NULL;
  coinFactorizationB_ = new CoinDenseFactorization();
  //printf("going dense\n");
}
int 
ClpFactorization::factorize ( ClpSimplex * model,
			      int solveType, bool valuesPass)
{
  ClpMatrixBase * matrix = model->clpMatrix(); 
  int numberRows = model->numberRows();
  int numberColumns = model->numberColumns();
  if (!numberRows)
    return 0;
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(-1);
#endif
  if (coinFactorizationB_) {
    setStatus(-99);
    int * pivotVariable = model->pivotVariable();
    //returns 0 -okay, -1 singular, -2 too many in basis */
    while (status()<-98) {
      
      int i;
      int numberBasic=0;
      int numberRowBasic;
      // Move pivot variables across if they look good
      int * pivotTemp = model->rowArray(0)->getIndices();
      assert (!model->rowArray(0)->getNumElements());
      if (!matrix->rhsOffset(model)) {
	// Seems to prefer things in order so quickest
	// way is to go though like this
	for (i=0;i<numberRows;i++) {
	  if (model->getRowStatus(i) == ClpSimplex::basic) 
	    pivotTemp[numberBasic++]=i;
	}
	numberRowBasic=numberBasic;
	/* Put column basic variables into pivotVariable
	   This is done by ClpMatrixBase to allow override for gub
	*/
	matrix->generalExpanded(model,0,numberBasic);
      } else {
	// Long matrix - do a different way
	bool fullSearch=false;
	for (i=0;i<numberRows;i++) {
	  int iPivot = pivotVariable[i];
	  if (iPivot>=numberColumns) {
	    pivotTemp[numberBasic++]=iPivot-numberColumns;
	  }
	}
	numberRowBasic=numberBasic;
	for (i=0;i<numberRows;i++) {
	  int iPivot = pivotVariable[i];
	  if (iPivot<numberColumns) {
	    if (iPivot>=0) {
	      pivotTemp[numberBasic++]=iPivot;
	    } else {
	      // not full basis
	      fullSearch=true;
	      break;
	    }
	  }
	}
	if (fullSearch) {
	  // do slow way
	  numberBasic=0;
	  for (i=0;i<numberRows;i++) {
	    if (model->getRowStatus(i) == ClpSimplex::basic) 
	      pivotTemp[numberBasic++]=i;
	  }
	  numberRowBasic=numberBasic;
	  /* Put column basic variables into pivotVariable
	     This is done by ClpMatrixBase to allow override for gub
	  */
	  matrix->generalExpanded(model,0,numberBasic);
	}
      }
      if (numberBasic>model->maximumBasic()) {
        // Take out some
        numberBasic=numberRowBasic;
        for (int i=0;i<numberColumns;i++) {
          if (model->getColumnStatus(i) == ClpSimplex::basic) {
            if (numberBasic<numberRows)
              numberBasic++;
            else
              model->setColumnStatus(i,ClpSimplex::superBasic);
          }
        }
        numberBasic=numberRowBasic;
        matrix->generalExpanded(model,0,numberBasic);
      }
      CoinBigIndex numberElements=numberRowBasic;
      
      // compute how much in basis
      // can change for gub
      int numberColumnBasic = numberBasic-numberRowBasic;
      
      numberElements +=matrix->countBasis(model,
					  pivotTemp+numberRowBasic, 
					  numberRowBasic,
					  numberColumnBasic);
      // Not needed for dense
      numberElements = 3 * numberBasic + 3 * numberElements + 20000;
      coinFactorizationB_->getAreas ( numberRows,
				      numberRowBasic+numberColumnBasic, numberElements,
				      2 * numberElements );
      // Fill in counts so we can skip part of preProcess
      // This is NOT needed for dense but would be needed for later versions
      double * elementU;
      int * indexRowU;
      CoinBigIndex * startColumnU;
      int * numberInRow;
      int * numberInColumn;
      elementU = coinFactorizationB_->elements();
      indexRowU = coinFactorizationB_->indices();
      startColumnU = coinFactorizationB_->starts();
#ifndef COIN_FAST_CODE
      double slackValue;
      slackValue = coinFactorizationB_->slackValue();
#else
#define slackValue -1.0
#endif
      // At present we don't need counts
      numberInRow = coinFactorizationB_->intWorkArea();
      numberInColumn = numberInRow;
      CoinZeroN ( numberInRow, numberRows + 1 ); // just for valgrind
      for (i=0;i<numberRowBasic;i++) {
	int iRow = pivotTemp[i];
	// Change pivotTemp to correct sequence
	pivotTemp[i]=iRow+numberColumns;
	indexRowU[i]=iRow;
	startColumnU[i]=i;
	elementU[i]=slackValue;
      }
      startColumnU[numberRowBasic]=numberRowBasic;
      // can change for gub so redo
      numberColumnBasic = numberBasic-numberRowBasic;
      matrix->fillBasis(model, 
			pivotTemp+numberRowBasic, 
			numberColumnBasic,
			indexRowU, 
			startColumnU+numberRowBasic,
			numberInRow,
			numberInColumn+numberRowBasic,
			elementU);
      // recompute number basic
      numberBasic = numberRowBasic+numberColumnBasic;
      for (i=numberBasic;i<numberRows;i++)
	pivotTemp[i]=-1; // mark not there
      if (numberBasic) 
	numberElements = startColumnU[numberBasic-1]
	  +numberInColumn[numberBasic-1];
      else
	numberElements=0;
      coinFactorizationB_->preProcess ( );
      coinFactorizationB_->factor (  );
      // If we get here status is 0 or -1
      if (status() == 0&&numberBasic==numberRows) {
	coinFactorizationB_->postProcess(pivotTemp,pivotVariable);
      } else {
	// Change pivotTemp to be correct list
	coinFactorizationB_->makeNonSingular(pivotTemp,numberColumns);
	double * columnLower = model->lowerRegion();
	double * columnUpper = model->upperRegion();
	double * columnActivity = model->solutionRegion();
	double * rowLower = model->lowerRegion(0);
	double * rowUpper = model->upperRegion(0);
	double * rowActivity = model->solutionRegion(0);
	//redo basis - first take ALL out
	int iColumn;
	double largeValue = model->largeValue();
	for (iColumn=0;iColumn<numberColumns;iColumn++) {
	  if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
	    // take out
	    if (!valuesPass) {
	      double lower=columnLower[iColumn];
	      double upper=columnUpper[iColumn];
	      double value=columnActivity[iColumn];
	      if (lower>-largeValue||upper<largeValue) {
		if (fabs(value-lower)<fabs(value-upper)) {
		  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
		  columnActivity[iColumn]=lower;
		} else {
		  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
		  columnActivity[iColumn]=upper;
		}
	      } else {
		model->setColumnStatus(iColumn,ClpSimplex::isFree);
	      }
	    } else {
	      model->setColumnStatus(iColumn,ClpSimplex::superBasic);
	    }
	  }
	}
	int iRow;
	for (iRow=0;iRow<numberRows;iRow++) {
	  if (model->getRowStatus(iRow)==ClpSimplex::basic) {
	    // take out
	    if (!valuesPass) {
	      double lower=columnLower[iRow];
	      double upper=columnUpper[iRow];
	      double value=columnActivity[iRow];
	      if (lower>-largeValue||upper<largeValue) {
		if (fabs(value-lower)<fabs(value-upper)) {
		  model->setRowStatus(iRow,ClpSimplex::atLowerBound);
		  columnActivity[iRow]=lower;
		} else {
		  model->setRowStatus(iRow,ClpSimplex::atUpperBound);
		  columnActivity[iRow]=upper;
		}
	      } else {
		model->setRowStatus(iRow,ClpSimplex::isFree);
	      }
	    } else {
	      model->setRowStatus(iRow,ClpSimplex::superBasic);
	    }
	  }
	}
	for (iRow=0;iRow<numberRows;iRow++) {
	  int iSequence=pivotTemp[iRow];
	  assert (iSequence>=0);
	  // basic
	  model->setColumnStatus(iSequence,ClpSimplex::basic);
	}
	// signal repeat
	coinFactorizationB_->setStatus(-99);
	// set fixed if they are
	for (iRow=0;iRow<numberRows;iRow++) {
	  if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
	    if (rowLower[iRow]==rowUpper[iRow]) {
	      rowActivity[iRow]=rowLower[iRow];
	      model->setRowStatus(iRow,ClpSimplex::isFixed);
	    }
	  }
	}
	for (iColumn=0;iColumn<numberColumns;iColumn++) {
	  if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
	    if (columnLower[iColumn]==columnUpper[iColumn]) {
	      columnActivity[iColumn]=columnLower[iColumn];
	      model->setColumnStatus(iColumn,ClpSimplex::isFixed);
	    }
	  }
	}
      } 
    }
#ifdef CLP_DEBUG
    // check basic
    CoinIndexedVector region1(2*numberRows);
    CoinIndexedVector region2B(2*numberRows);
    int iPivot;
    double * arrayB = region2B.denseVector();
    int i;
    for (iPivot=0;iPivot<numberRows;iPivot++) {
      int iSequence = pivotVariable[iPivot];
      model->unpack(&region2B,iSequence);
      coinFactorizationB_->updateColumn(&region1,&region2B);
      if (fabs(arrayB[iPivot]-1.0)<1.0e-4) {
	// OK?
	arrayB[iPivot]=0.0;
      } else {
	assert (fabs(arrayB[iPivot])<1.0e-4);
	for (i=0;i<numberRows;i++) {
	  if (fabs(arrayB[i]-1.0)<1.0e-4)
	    break;
	}
	assert (i<numberRows);
	printf("variable on row %d landed up on row %d\n",iPivot,i);
	arrayB[i]=0.0;
      }
      for (i=0;i<numberRows;i++)
	assert (fabs(arrayB[i])<1.0e-4);
      region2B.clear();
    }
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(2);
#endif
    return coinFactorizationB_->status();
  }
  // If too many compressions increase area
  if (coinFactorizationA_->pivots()>1&&coinFactorizationA_->numberCompressions()*10 > coinFactorizationA_->pivots()+10) {
    coinFactorizationA_->areaFactor( coinFactorizationA_->areaFactor()* 1.1);
  }
  //int numberPivots=coinFactorizationA_->pivots();
#if 0
  if (model->algorithm()>0)
    numberSave=-1;
#endif
#ifndef SLIM_CLP
  if (!networkBasis_||doCheck) {
#endif
    coinFactorizationA_->setStatus(-99);
    int * pivotVariable = model->pivotVariable();
    int nTimesRound=0;
    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
    while (coinFactorizationA_->status()<-98) {
      nTimesRound++;
      
      int i;
      int numberBasic=0;
      int numberRowBasic;
      // Move pivot variables across if they look good
      int * pivotTemp = model->rowArray(0)->getIndices();
      assert (!model->rowArray(0)->getNumElements());
      if (!matrix->rhsOffset(model)) {
#if 0
        if (numberSave>0) {
          int nStill=0;
          int nAtBound=0;
          int nZeroDual=0;
          CoinIndexedVector * array = model->rowArray(3);
          CoinIndexedVector * objArray = model->columnArray(1);
          array->clear();
          objArray->clear();
          double * cost = model->costRegion();
          double tolerance = model->primalTolerance();
          double offset=0.0;
          for (i=0;i<numberRows;i++) {
            int iPivot = pivotVariable[i];
            if (iPivot<numberColumns&&isDense(iPivot)) {
              if (model->getColumnStatus(iPivot)==ClpSimplex::basic) {
                nStill++;
                double value=model->solutionRegion()[iPivot];
                double dual = model->dualRowSolution()[i];
                double lower=model->lowerRegion()[iPivot];
                double upper=model->upperRegion()[iPivot];
                ClpSimplex::Status status;
                if (fabs(value-lower)<tolerance) {
                  status=ClpSimplex::atLowerBound;
                  nAtBound++;
                } else if (fabs(value-upper)<tolerance) {
                  nAtBound++;
                  status=ClpSimplex::atUpperBound;
                } else if (value>lower&&value<upper) {
                  status=ClpSimplex::superBasic;
                } else {
                  status=ClpSimplex::basic;
                }
                if (status!=ClpSimplex::basic) {
                  if (model->getRowStatus(i)!=ClpSimplex::basic) {
                    model->setColumnStatus(iPivot,ClpSimplex::atLowerBound);
                    model->setRowStatus(i,ClpSimplex::basic);
                    pivotVariable[i]=i+numberColumns;
                    model->dualRowSolution()[i]=0.0;
                    model->djRegion(0)[i]=0.0;
                    array->add(i,dual);
                    offset += dual*model->solutionRegion(0)[i];
                  }
                }
                if (fabs(dual)<1.0e-5)
                  nZeroDual++;
              }
            }
          }
          printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
                 numberSave,nStill,nAtBound,nZeroDual,offset);
          if (array->getNumElements()) {
            // modify costs
            model->clpMatrix()->transposeTimes(model,1.0,array,model->columnArray(0),
                                               objArray);
            array->clear();
            int n=objArray->getNumElements();
            int * indices = objArray->getIndices();
            double * elements = objArray->denseVector();
            for (i=0;i<n;i++) {
              int iColumn = indices[i];
              cost[iColumn] -= elements[iColumn];
              elements[iColumn]=0.0;
            }
            objArray->setNumElements(0);
          }
        }
#endif
	// Seems to prefer things in order so quickest
	// way is to go though like this
	for (i=0;i<numberRows;i++) {
	  if (model->getRowStatus(i) == ClpSimplex::basic) 
	    pivotTemp[numberBasic++]=i;
	}
	numberRowBasic=numberBasic;
	/* Put column basic variables into pivotVariable
	   This is done by ClpMatrixBase to allow override for gub
	*/
	matrix->generalExpanded(model,0,numberBasic);
      } else {
	// Long matrix - do a different way
	bool fullSearch=false;
	for (i=0;i<numberRows;i++) {
	  int iPivot = pivotVariable[i];
	  if (iPivot>=numberColumns) {
	    pivotTemp[numberBasic++]=iPivot-numberColumns;
	  }
	}
	numberRowBasic=numberBasic;
	for (i=0;i<numberRows;i++) {
	  int iPivot = pivotVariable[i];
	  if (iPivot<numberColumns) {
	    if (iPivot>=0) {
	      pivotTemp[numberBasic++]=iPivot;
	    } else {
	      // not full basis
	      fullSearch=true;
	      break;
	    }
	  }
	}
	if (fullSearch) {
	  // do slow way
	  numberBasic=0;
	  for (i=0;i<numberRows;i++) {
	    if (model->getRowStatus(i) == ClpSimplex::basic) 
	      pivotTemp[numberBasic++]=i;
	  }
	  numberRowBasic=numberBasic;
	  /* Put column basic variables into pivotVariable
	     This is done by ClpMatrixBase to allow override for gub
	  */
	  matrix->generalExpanded(model,0,numberBasic);
	}
      }
      if (numberBasic>model->maximumBasic()) {
#if 0 // ndef NDEBUG
        printf("%d basic - should only be %d\n",
               numberBasic,numberRows);
#endif
        // Take out some
        numberBasic=numberRowBasic;
        for (int i=0;i<numberColumns;i++) {
          if (model->getColumnStatus(i) == ClpSimplex::basic) {
            if (numberBasic<numberRows)
              numberBasic++;
            else
              model->setColumnStatus(i,ClpSimplex::superBasic);
          }
        }
        numberBasic=numberRowBasic;
        matrix->generalExpanded(model,0,numberBasic);
      }
#ifndef SLIM_CLP
      // see if matrix a network
#ifndef NO_RTTI
      ClpNetworkMatrix* networkMatrix =
	dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
#else
      ClpNetworkMatrix* networkMatrix = NULL;
      if (model->clpMatrix()->type()==11)
	networkMatrix = 
	static_cast< ClpNetworkMatrix*>(model->clpMatrix());
#endif
      // If network - still allow ordinary factorization first time for laziness
      if (networkMatrix)
	coinFactorizationA_->setBiasLU(0); // All to U if network
      //int saveMaximumPivots = maximumPivots();
      delete networkBasis_;
      networkBasis_ = NULL;
      if (networkMatrix&&!doCheck)
	maximumPivots(1);
#endif
      //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
      while (coinFactorizationA_->status()==-99) {
	// maybe for speed will be better to leave as many regions as possible
	coinFactorizationA_->gutsOfDestructor();
	coinFactorizationA_->gutsOfInitialize(2);
	CoinBigIndex numberElements=numberRowBasic;

	// compute how much in basis

	int i;
	// can change for gub
	int numberColumnBasic = numberBasic-numberRowBasic;

	numberElements +=matrix->countBasis(model,
					   pivotTemp+numberRowBasic, 
					   numberRowBasic,
					    numberColumnBasic);
	// and recompute as network side say different
	if (model->numberIterations())
	  numberRowBasic = numberBasic - numberColumnBasic;
	numberElements = 3 * numberBasic + 3 * numberElements + 20000;
	coinFactorizationA_->getAreas ( numberRows,
		   numberRowBasic+numberColumnBasic, numberElements,
		   2 * numberElements );
	//fill
	// Fill in counts so we can skip part of preProcess
	int * numberInRow = coinFactorizationA_->numberInRow();
	int * numberInColumn = coinFactorizationA_->numberInColumn();
	CoinZeroN ( numberInRow, coinFactorizationA_->numberRows() + 1 );
	CoinZeroN ( numberInColumn, coinFactorizationA_->maximumColumnsExtra() + 1 );
	double * elementU = coinFactorizationA_->elementU();
	int * indexRowU = coinFactorizationA_->indexRowU();
	CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
#ifndef COIN_FAST_CODE
	double slackValue = coinFactorizationA_->slackValue();
#endif
	for (i=0;i<numberRowBasic;i++) {
	  int iRow = pivotTemp[i];
	  indexRowU[i]=iRow;
	  startColumnU[i]=i;
	  elementU[i]=slackValue;
	  numberInRow[iRow]=1;
	  numberInColumn[i]=1;
	}
	startColumnU[numberRowBasic]=numberRowBasic;
	// can change for gub so redo
	numberColumnBasic = numberBasic-numberRowBasic;
	matrix->fillBasis(model, 
			  pivotTemp+numberRowBasic, 
			  numberColumnBasic,
			  indexRowU, 
			  startColumnU+numberRowBasic,
			  numberInRow,
			  numberInColumn+numberRowBasic,
			  elementU);
#if 0
	{
	  printf("%d row basic, %d column basic\n",numberRowBasic,numberColumnBasic);
	  for (int i=0;i<numberElements;i++) 
	    printf("row %d col %d value %g\n",indexRowU[i],indexColumnU_[i],
		   elementU[i]);
	}
#endif
	// recompute number basic
        numberBasic = numberRowBasic+numberColumnBasic;
	if (numberBasic) 
	  numberElements = startColumnU[numberBasic-1]
	    +numberInColumn[numberBasic-1];
	else
	  numberElements=0;
	coinFactorizationA_->setNumberElementsU(numberElements);
        //saveFactorization("dump.d");
	if (coinFactorizationA_->biasLU()>=3||coinFactorizationA_->numberRows()!=coinFactorizationA_->numberColumns())
	  coinFactorizationA_->preProcess ( 2 );
	else
	  coinFactorizationA_->preProcess ( 3 ); // no row copy
	coinFactorizationA_->factor (  );
	if (coinFactorizationA_->status()==-99) {
	  // get more memory
	  coinFactorizationA_->areaFactor(2.0*coinFactorizationA_->areaFactor());
	} else if (coinFactorizationA_->status()==-1&&
		   (model->numberIterations()==0||nTimesRound>2)&&
                   coinFactorizationA_->denseThreshold()) {
          // Round again without dense
          coinFactorizationA_->setDenseThreshold(0);
          coinFactorizationA_->setStatus(-99);
        }
      }
      // If we get here status is 0 or -1
      if (coinFactorizationA_->status() == 0) {
	// We may need to tamper with order and redo - e.g. network with side
	int useNumberRows = numberRows;
	// **** we will also need to add test in dual steepest to do
	// as we do for network
        matrix->generalExpanded(model,12,useNumberRows);
	const int * permuteBack = coinFactorizationA_->permuteBack();
	const int * back = coinFactorizationA_->pivotColumnBack();
	//int * pivotTemp = pivotColumn_.array();
	//ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
	// Redo pivot order
	for (i=0;i<numberRowBasic;i++) {
	  int k = pivotTemp[i];
	  // so rowIsBasic[k] would be permuteBack[back[i]]
	  pivotVariable[permuteBack[back[i]]]=k+numberColumns;
	}
	for (;i<useNumberRows;i++) {
	  int k = pivotTemp[i];
	  // so rowIsBasic[k] would be permuteBack[back[i]]
	  pivotVariable[permuteBack[back[i]]]=k;
	}
#if 0
        if (numberSave>=0) {
          numberSave=numberDense_;
          memset(saveList,0,((coinFactorizationA_->numberRows()+31)>>5)*sizeof(int));
          for (i=coinFactorizationA_->numberRows()-numberSave;i<coinFactorizationA_->numberRows();i++) {
            int k=pivotTemp[pivotColumn_.array()[i]];
            setDense(k);
          }
        }
#endif
	// Set up permutation vector
	// these arrays start off as copies of permute
	// (and we could use permute_ instead of pivotColumn (not back though))
	ClpDisjointCopyN ( coinFactorizationA_->permute(), useNumberRows , coinFactorizationA_->pivotColumn()  );
	ClpDisjointCopyN ( coinFactorizationA_->permuteBack(), useNumberRows , coinFactorizationA_->pivotColumnBack()  );
#ifndef SLIM_CLP
	if (networkMatrix) {
	  maximumPivots(CoinMax(2000,maximumPivots()));
	  // redo arrays
	  for (int iRow=0;iRow<4;iRow++) {
	    int length =model->numberRows()+maximumPivots();
	    if (iRow==3||model->objectiveAsObject()->type()>1)
	      length += model->numberColumns();
	    model->rowArray(iRow)->reserve(length);
	  }
	  // create network factorization
	  if (doCheck)
	    delete networkBasis_; // temp
	  networkBasis_ = new ClpNetworkBasis(model,coinFactorizationA_->numberRows(),
					      coinFactorizationA_->pivotRegion(),
					      coinFactorizationA_->permuteBack(),
					      coinFactorizationA_->startColumnU(),
					      coinFactorizationA_->numberInColumn(),
					      coinFactorizationA_->indexRowU(),
					      coinFactorizationA_->elementU());
	  // kill off arrays in ordinary factorization
	  if (!doCheck) {
	    coinFactorizationA_->gutsOfDestructor();
	    // but make sure coinFactorizationA_->numberRows() set
	    coinFactorizationA_->setNumberRows(model->numberRows());
	    coinFactorizationA_->setStatus(0);
#if 0
	    // but put back permute arrays so odd things will work
	    int numberRows = model->numberRows();
	    pivotColumnBack_ = new int [numberRows];
	    permute_ = new int [numberRows];
	    int i;
	    for (i=0;i<numberRows;i++) {
	      pivotColumnBack_[i]=i;
	      permute_[i]=i;
	    }
#endif
	  }
	} else {
#endif
	  // See if worth going sparse and when
	  coinFactorizationA_->checkSparse();
#ifndef SLIM_CLP
	}
#endif
      } else if (coinFactorizationA_->status() == -1&&(solveType==0||solveType==2)) {
	// This needs redoing as it was merged coding - does not need array
	int numberTotal = numberRows+numberColumns;
	int * isBasic = new int [numberTotal];
	int * rowIsBasic = isBasic+numberColumns;
	int * columnIsBasic = isBasic;
	for (i=0;i<numberTotal;i++) 
	  isBasic[i]=-1;
	for (i=0;i<numberRowBasic;i++) {
	  int iRow = pivotTemp[i];
	  rowIsBasic[iRow]=1;
	}
	for (;i<numberBasic;i++) {
	  int iColumn = pivotTemp[i];
	  columnIsBasic[iColumn]=1;
	}
	numberBasic=0;
	for (i=0;i<numberRows;i++) 
	  pivotVariable[i]=-1;
	// mark as basic or non basic
	const int * pivotColumn = coinFactorizationA_->pivotColumn();
	for (i=0;i<numberRows;i++) {
	  if (rowIsBasic[i]>=0) {
	    if (pivotColumn[numberBasic]>=0) {
	      rowIsBasic[i]=pivotColumn[numberBasic];
	    } else {
	      rowIsBasic[i]=-1;
              model->setRowStatus(i,ClpSimplex::superBasic);
            }
	    numberBasic++;
	  }
	}
	for (i=0;i<numberColumns;i++) {
	  if (columnIsBasic[i]>=0) {
	    if (pivotColumn[numberBasic]>=0) 
	      columnIsBasic[i]=pivotColumn[numberBasic];
	    else
	      columnIsBasic[i]=-1;
	    numberBasic++;
	  }
	}
	// leave pivotVariable in useful form for cleaning basis
	int * pivotVariable = model->pivotVariable();
	for (i=0;i<numberRows;i++) {
	  pivotVariable[i]=-1;
	}

	for (i=0;i<numberRows;i++) {
	  if (model->getRowStatus(i) == ClpSimplex::basic) {
	    int iPivot = rowIsBasic[i];
	    if (iPivot>=0) 
	      pivotVariable[iPivot]=i+numberColumns;
	  }
	}
	for (i=0;i<numberColumns;i++) {
	  if (model->getColumnStatus(i) == ClpSimplex::basic) {
	    int iPivot = columnIsBasic[i];
	    if (iPivot>=0) 
	      pivotVariable[iPivot]=i;
	  }
	}
	delete [] isBasic;
	double * columnLower = model->lowerRegion();
	double * columnUpper = model->upperRegion();
	double * columnActivity = model->solutionRegion();
	double * rowLower = model->lowerRegion(0);
	double * rowUpper = model->upperRegion(0);
	double * rowActivity = model->solutionRegion(0);
	//redo basis - first take ALL columns out
	int iColumn;
	double largeValue = model->largeValue();
	for (iColumn=0;iColumn<numberColumns;iColumn++) {
	  if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
	    // take out
	    if (!valuesPass) {
	      double lower=columnLower[iColumn];
	      double upper=columnUpper[iColumn];
	      double value=columnActivity[iColumn];
	      if (lower>-largeValue||upper<largeValue) {
		if (fabs(value-lower)<fabs(value-upper)) {
		  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
		  columnActivity[iColumn]=lower;
		} else {
		  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
		  columnActivity[iColumn]=upper;
		}
	      } else {
		model->setColumnStatus(iColumn,ClpSimplex::isFree);
	      }
	    } else {
	      model->setColumnStatus(iColumn,ClpSimplex::superBasic);
	    }
	  }
	}
	int iRow;
	for (iRow=0;iRow<numberRows;iRow++) {
	  int iSequence=pivotVariable[iRow];
	  if (iSequence>=0) {
	    // basic
	    if (iSequence>=numberColumns) {
	      // slack in - leave
	      //assert (iSequence-numberColumns==iRow);
	    } else {
              assert(model->getRowStatus(iRow)!=ClpSimplex::basic);
	      // put back structural
	      model->setColumnStatus(iSequence,ClpSimplex::basic);
	    }
	  } else {
	    // put in slack
	    model->setRowStatus(iRow,ClpSimplex::basic);
	  }
	}
	// Put back any key variables for gub 
	int dummy;
	matrix->generalExpanded(model,1,dummy);
	// signal repeat
	coinFactorizationA_->setStatus(-99);
	// set fixed if they are
	for (iRow=0;iRow<numberRows;iRow++) {
	  if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
	    if (rowLower[iRow]==rowUpper[iRow]) {
	      rowActivity[iRow]=rowLower[iRow];
	      model->setRowStatus(iRow,ClpSimplex::isFixed);
	    }
	  }
	}
	for (iColumn=0;iColumn<numberColumns;iColumn++) {
	  if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
	    if (columnLower[iColumn]==columnUpper[iColumn]) {
	      columnActivity[iColumn]=columnLower[iColumn];
	      model->setColumnStatus(iColumn,ClpSimplex::isFixed);
	    }
	  }
	}
      } 
    }
#ifndef SLIM_CLP
  } else {
    // network - fake factorization - do nothing
    coinFactorizationA_->setStatus(0);
    coinFactorizationA_->setPivots(0);
  }
#endif
#ifndef SLIM_CLP
  if (!coinFactorizationA_->status()) {
    // take out part if quadratic
    if (model->algorithm()==2) {
      ClpObjective * obj = model->objectiveAsObject();
      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
      assert (quadraticObj);
      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
      int numberXColumns = quadratic->getNumCols();
      assert (numberXColumns<numberColumns);
      int base = numberColumns-numberXColumns;
      int * which = new int [numberXColumns];
      int * pivotVariable = model->pivotVariable();
      int * permute = pivotColumn();
      int i;
      int n=0;
      for (i=0;i<numberRows;i++) {
	int iSj = pivotVariable[i]-base;
	if (iSj>=0&&iSj<numberXColumns) 
	  which[n++]=permute[i];
      }
      if (n)
	coinFactorizationA_->emptyRows(n,which);
      delete [] which;
    }
  }
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(2);
#endif
  return coinFactorizationA_->status();
}
/* Replaces one Column to basis,
   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
   If checkBeforeModifying is true will do all accuracy checks
   before modifying factorization.  Whether to set this depends on
   speed considerations.  You could just do this on first iteration
   after factorization and thereafter re-factorize
   partial update already in U */
int 
ClpFactorization::replaceColumn ( const ClpSimplex * model, 
				  CoinIndexedVector * regionSparse,
				  CoinIndexedVector * tableauColumn,
				  int pivotRow,
				  double pivotCheck ,
				  bool checkBeforeModifying)
{
#ifndef SLIM_CLP
  if (!networkBasis_) {
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(-1);
#endif
    int returnCode;
    // see if FT
    if (!coinFactorizationA_||coinFactorizationA_->forrestTomlin()) {
      if (coinFactorizationA_) {
	returnCode = coinFactorizationA_->replaceColumn(regionSparse,
							pivotRow,
							pivotCheck,
							checkBeforeModifying);
      } else {
	returnCode= coinFactorizationB_->replaceColumn(tableauColumn,
					      pivotRow,
					      pivotCheck,
					      checkBeforeModifying);
#ifdef CLP_DEBUG
	// check basic
	int numberRows = coinFactorizationB_->numberRows();
	CoinIndexedVector region1(2*numberRows);
	CoinIndexedVector region2A(2*numberRows);
	CoinIndexedVector region2B(2*numberRows);
	int iPivot;
	double * arrayB = region2B.denseVector();
	int * pivotVariable = model->pivotVariable();
	int i;
	for (iPivot=0;iPivot<numberRows;iPivot++) {
	  int iSequence = pivotVariable[iPivot];
	  if (iPivot==pivotRow)
	    iSequence = model->sequenceIn();
	  model->unpack(&region2B,iSequence);
	  coinFactorizationB_->updateColumn(&region1,&region2B);
	  assert (fabs(arrayB[iPivot]-1.0)<1.0e-4);
	  arrayB[iPivot]=0.0;
	  for (i=0;i<numberRows;i++)
	    assert (fabs(arrayB[i])<1.0e-4);
	  region2B.clear();
	}
#endif
      }
    } else {
      returnCode = coinFactorizationA_->replaceColumnPFI(tableauColumn,
					      pivotRow,pivotCheck); // Note array
    }
#ifdef CLP_FACTORIZATION_INSTRUMENT
      factorization_instrument(3);
#endif
      return returnCode;

#ifndef SLIM_CLP
  } else {
    if (doCheck) {
      int returnCode = coinFactorizationA_->replaceColumn(regionSparse,
      						pivotRow,
      						pivotCheck,
      						checkBeforeModifying);
      networkBasis_->replaceColumn(regionSparse,
				   pivotRow);
      return returnCode;
    } else {
      // increase number of pivots
      coinFactorizationA_->setPivots(coinFactorizationA_->pivots()+1);
      return networkBasis_->replaceColumn(regionSparse,
				   pivotRow);
    }
  }
#endif
}

/* Updates one column (FTRAN) from region2
   number returned is negative if no room
   region1 starts as zero and is zero at end */
int 
ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
				   CoinIndexedVector * regionSparse2)
{
#ifdef CLP_DEBUG
  regionSparse->checkClear();
#endif
  if (!numberRows())
      return 0;
#ifndef SLIM_CLP
  if (!networkBasis_) {
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(-1);
#endif
    int returnCode;
    if (coinFactorizationA_) {
      coinFactorizationA_->setCollectStatistics(true);
      returnCode= coinFactorizationA_->updateColumnFT(regionSparse,
							   regionSparse2);
      coinFactorizationA_->setCollectStatistics(false);
    } else {
      returnCode= coinFactorizationB_->updateColumnFT(regionSparse,
							 regionSparse2);
    }
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(4);
#endif
    return returnCode;
#ifndef SLIM_CLP
  } else {
#ifdef CHECK_NETWORK
    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
    double * check = new double[coinFactorizationA_->numberRows()];
    int returnCode = coinFactorizationA_->updateColumnFT(regionSparse,
						       regionSparse2);
    networkBasis_->updateColumn(regionSparse,save,-1);
    int i;
    double * array = regionSparse2->denseVector();
    int * indices = regionSparse2->getIndices();
    int n=regionSparse2->getNumElements();
    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
    double * array2 = save->denseVector();
    int * indices2 = save->getIndices();
    int n2=save->getNumElements();
    assert (n==n2);
    if (save->packedMode()) {
      for (i=0;i<n;i++) {
	check[indices[i]]=array[i];
      }
      for (i=0;i<n;i++) {
	double value2 = array2[i];
	assert (check[indices2[i]]==value2);
      }
    } else {
      int numberRows = coinFactorizationA_->numberRows();
      for (i=0;i<numberRows;i++) {
	double value1 = array[i];
	double value2 = array2[i];
	assert (value1==value2);
      }
    }
    delete save;
    delete [] check;
    return returnCode;
#else
    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
    return 1;
#endif
  }
#endif
}
/* Updates one column (FTRAN) from region2
   number returned is negative if no room
   region1 starts as zero and is zero at end */
int 
ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
				 CoinIndexedVector * regionSparse2,
				 bool noPermute) const
{
#ifdef CLP_DEBUG
  if (!noPermute)
    regionSparse->checkClear();
#endif
  if (!numberRows())
    return 0;
#ifndef SLIM_CLP
  if (!networkBasis_) {
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(-1);
#endif
    int returnCode;
    if (coinFactorizationA_) {
      coinFactorizationA_->setCollectStatistics(true);
      returnCode= coinFactorizationA_->updateColumn(regionSparse,
							 regionSparse2,
							 noPermute);
      coinFactorizationA_->setCollectStatistics(false);
    } else {
      returnCode= coinFactorizationB_->updateColumn(regionSparse,
							 regionSparse2,
							 noPermute);
    }
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(5);
#endif
    return returnCode;
#ifndef SLIM_CLP
  } else {
#ifdef CHECK_NETWORK
    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
    double * check = new double[coinFactorizationA_->numberRows()];
    int returnCode = coinFactorizationA_->updateColumn(regionSparse,
						     regionSparse2,
						     noPermute);
    networkBasis_->updateColumn(regionSparse,save,-1);
    int i;
    double * array = regionSparse2->denseVector();
    int * indices = regionSparse2->getIndices();
    int n=regionSparse2->getNumElements();
    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
    double * array2 = save->denseVector();
    int * indices2 = save->getIndices();
    int n2=save->getNumElements();
    assert (n==n2);
    if (save->packedMode()) {
      for (i=0;i<n;i++) {
	check[indices[i]]=array[i];
      }
      for (i=0;i<n;i++) {
	double value2 = array2[i];
	assert (check[indices2[i]]==value2);
      }
    } else {
      int numberRows = coinFactorizationA_->numberRows();
      for (i=0;i<numberRows;i++) {
	double value1 = array[i];
	double value2 = array2[i];
	assert (value1==value2);
      }
    }
    delete save;
    delete [] check;
    return returnCode;
#else
    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
    return 1;
#endif
  }
#endif
}
/* Updates one column (FTRAN) from region2
   Tries to do FT update
   number returned is negative if no room.
   Also updates region3
   region1 starts as zero and is zero at end */
int 
ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
				       CoinIndexedVector * regionSparse2,
				       CoinIndexedVector * regionSparse3,
				       bool noPermuteRegion3)
{
#ifdef CLP_DEBUG
  regionSparse1->checkClear();
#endif
  if (!numberRows())
    return 0;
  int returnCode=0;
#ifndef SLIM_CLP
  if (!networkBasis_) {
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(-1);
#endif
    if (coinFactorizationA_) {
      coinFactorizationA_->setCollectStatistics(true);
      if (coinFactorizationA_->spaceForForrestTomlin()) {
	assert (regionSparse2->packedMode());
	assert (!regionSparse3->packedMode());
	returnCode= coinFactorizationA_->updateTwoColumnsFT(regionSparse1,
							    regionSparse2,
							    regionSparse3,
							    noPermuteRegion3);
      } else {
	returnCode= coinFactorizationA_->updateColumnFT(regionSparse1,
							regionSparse2);
	coinFactorizationA_->updateColumn(regionSparse1,
					  regionSparse3,
					  noPermuteRegion3);
      }
      coinFactorizationA_->setCollectStatistics(false);
    } else {
      returnCode= coinFactorizationB_->updateColumnFT(regionSparse1,
						     regionSparse2);
      coinFactorizationB_->updateColumn(regionSparse1,
					regionSparse3,
					noPermuteRegion3);
    }
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(9);
#endif
    return returnCode;
#ifndef SLIM_CLP
  } else {
    returnCode=updateColumnFT(regionSparse1,regionSparse2);
    updateColumn(regionSparse1,regionSparse3,noPermuteRegion3);
  }
#endif
  return returnCode;
}
/* Updates one column (FTRAN) from region2
   number returned is negative if no room
   region1 starts as zero and is zero at end */
int 
ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
				 CoinIndexedVector * regionSparse2,
				 bool noPermute) const
{
  if (!noPermute)
    regionSparse->checkClear();
  if (!coinFactorizationA_->numberRows())
    return 0;
  coinFactorizationA_->setCollectStatistics(false);
  int returnCode= coinFactorizationA_->updateColumn(regionSparse,
                                                   regionSparse2,
                                                   noPermute);
  return returnCode;
}
/* Updates one column (BTRAN) from region2
   region1 starts as zero and is zero at end */
int 
ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
    			  CoinIndexedVector * regionSparse2) const
{
  if (!numberRows())
    return 0;
#ifndef SLIM_CLP
  if (!networkBasis_) {
#endif
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(-1);
#endif
    int returnCode;
    if (coinFactorizationA_) {
      coinFactorizationA_->setCollectStatistics(true);
      returnCode =  coinFactorizationA_->updateColumnTranspose(regionSparse,
						    regionSparse2);
      coinFactorizationA_->setCollectStatistics(false);
    } else {
      returnCode= coinFactorizationB_->updateColumnTranspose(regionSparse,
							 regionSparse2);
    }
#ifdef CLP_FACTORIZATION_INSTRUMENT
    factorization_instrument(6);
#endif
    return returnCode;
#ifndef SLIM_CLP
  } else {
#ifdef CHECK_NETWORK
    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
    double * check = new double[coinFactorizationA_->numberRows()];
    int returnCode = coinFactorizationA_->updateColumnTranspose(regionSparse,
							      regionSparse2);
    networkBasis_->updateColumnTranspose(regionSparse,save);
    int i;
    double * array = regionSparse2->denseVector();
    int * indices = regionSparse2->getIndices();
    int n=regionSparse2->getNumElements();
    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
    double * array2 = save->denseVector();
    int * indices2 = save->getIndices();
    int n2=save->getNumElements();
    assert (n==n2);
    if (save->packedMode()) {
      for (i=0;i<n;i++) {
	check[indices[i]]=array[i];
      }
      for (i=0;i<n;i++) {
	double value2 = array2[i];
	assert (check[indices2[i]]==value2);
      }
    } else {
      int numberRows = coinFactorizationA_->numberRows();
      for (i=0;i<numberRows;i++) {
	double value1 = array[i];
	double value2 = array2[i];
	assert (value1==value2);
      }
    }
    delete save;
    delete [] check;
    return returnCode;
#else
    return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
#endif
  }
#endif
}
/* makes a row copy of L for speed and to allow very sparse problems */
void 
ClpFactorization::goSparse()
{
#ifndef SLIM_CLP
  if (!networkBasis_) {
#endif
    if (coinFactorizationA_) {
#ifdef CLP_FACTORIZATION_INSTRUMENT
      factorization_instrument(-1);
#endif
      coinFactorizationA_->goSparse();
#ifdef CLP_FACTORIZATION_INSTRUMENT
      factorization_instrument(7);
#endif
    }
  }
}
// Cleans up i.e. gets rid of network basis 
void 
ClpFactorization::cleanUp()
{
#ifndef SLIM_CLP
  delete networkBasis_;
  networkBasis_=NULL;
#endif
  if (coinFactorizationA_)
  coinFactorizationA_->resetStatistics();
}
/// Says whether to redo pivot order
bool 
ClpFactorization::needToReorder() const
{
#ifdef CHECK_NETWORK
  return true;
#endif
#ifndef SLIM_CLP
  if (!networkBasis_)
#endif
    return true;
#ifndef SLIM_CLP
  else
    return false;
#endif
}
// Get weighted row list 
void
ClpFactorization::getWeights(int * weights) const
{
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(-1);
#endif
#ifndef SLIM_CLP
  if (networkBasis_) {
    // Network - just unit
    int numberRows = coinFactorizationA_->numberRows();
    for (int i=0;i<numberRows;i++) 
      weights[i]=1;
    return;
  }
#endif
  int * numberInRow = coinFactorizationA_->numberInRow();
  int * numberInColumn = coinFactorizationA_->numberInColumn();
  int * permuteBack = coinFactorizationA_->pivotColumnBack();
  int * indexRowU = coinFactorizationA_->indexRowU();
  const CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
  const CoinBigIndex * startRowL = coinFactorizationA_->startRowL();
  int numberRows = coinFactorizationA_->numberRows();
  if (!startRowL||!coinFactorizationA_->numberInRow()) {
    int * temp = new int[numberRows];
    memset(temp,0,numberRows*sizeof(int));
    int i;
    for (i=0;i<numberRows;i++) {
      // one for pivot
      temp[i]++;
      CoinBigIndex j;
      for (j=startColumnU[i];j<startColumnU[i]+numberInColumn[i];j++) {
	int iRow=indexRowU[j];
	temp[iRow]++;
      }
    }
    CoinBigIndex * startColumnL = coinFactorizationA_->startColumnL();
    int * indexRowL = coinFactorizationA_->indexRowL();
    int numberL = coinFactorizationA_->numberL();
    CoinBigIndex baseL = coinFactorizationA_->baseL();
    for (i=baseL;i<baseL+numberL;i++) {
      CoinBigIndex j;
      for (j=startColumnL[i];j<startColumnL[i+1];j++) {
	int iRow = indexRowL[j];
	temp[iRow]++;
      }
    }
    for (i=0;i<numberRows;i++) {
      int number = temp[i];
      int iPermute = permuteBack[i];
      weights[iPermute]=number;
    }
    delete [] temp;
  } else {
    int i;
    for (i=0;i<numberRows;i++) {
      int number = startRowL[i+1]-startRowL[i]+numberInRow[i]+1;
      //number = startRowL[i+1]-startRowL[i]+1;
      //number = numberInRow[i]+1;
      int iPermute = permuteBack[i];
      weights[iPermute]=number;
    }
  }
#ifdef CLP_FACTORIZATION_INSTRUMENT
  factorization_instrument(8);
#endif
}
#endif
