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


#include <cstdio>

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

#include "ClpSimplex.hpp"
#include "ClpFactorization.hpp"
// at end to get min/max!
#include "ClpNetworkMatrix.hpp"
#include "ClpPlusMinusOneMatrix.hpp"
#include "ClpMessage.hpp"
#include <iostream>
#include <cassert>

//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
ClpNetworkMatrix::ClpNetworkMatrix () 
  : ClpMatrixBase()
{
  setType(11);
  matrix_ = NULL;
  lengths_=NULL;
  indices_=NULL;
  numberRows_=0;
  numberColumns_=0;
  trueNetwork_=false;
}

/* Constructor from two arrays */
ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int * head,
				   const int * tail)
  : ClpMatrixBase()
{
  setType(11);
  matrix_ = NULL;
  lengths_=NULL;
  indices_=new int[2*numberColumns];;
  numberRows_=-1;
  numberColumns_=numberColumns;
  trueNetwork_=true;
  int iColumn;
  CoinBigIndex j=0;
  for (iColumn=0;iColumn<numberColumns_;iColumn++, j+=2) {
    int iRow = head[iColumn];
    numberRows_ = CoinMax(numberRows_,iRow);
    indices_[j]=iRow;
    iRow = tail[iColumn];
    numberRows_ = CoinMax(numberRows_,iRow);
    indices_[j+1]=iRow;
  }
  numberRows_++;
}
//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
ClpNetworkMatrix::ClpNetworkMatrix (const ClpNetworkMatrix & rhs) 
: ClpMatrixBase(rhs)
{  
  matrix_ = NULL;
  lengths_=NULL;
  indices_=NULL;
  numberRows_=rhs.numberRows_;
  numberColumns_=rhs.numberColumns_;
  trueNetwork_=rhs.trueNetwork_;
  if (numberColumns_) {
    indices_ = new int [ 2*numberColumns_];
    memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
  }
  int numberRows = getNumRows();
  if (rhs.rhsOffset_&&numberRows) {
    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
  } else {
    rhsOffset_=NULL;
  }
}

ClpNetworkMatrix::ClpNetworkMatrix (const CoinPackedMatrix & rhs) 
  : ClpMatrixBase()
{  
  setType(11);
  matrix_ = NULL;
  lengths_=NULL;
  indices_=NULL;
  int iColumn;
  assert (rhs.isColOrdered());
  // get matrix data pointers
  const int * row = rhs.getIndices();
  const CoinBigIndex * columnStart = rhs.getVectorStarts();
  const int * columnLength = rhs.getVectorLengths(); 
  const double * elementByColumn = rhs.getElements();
  numberColumns_ = rhs.getNumCols();
  int goodNetwork=1;
  numberRows_=-1;
  indices_ = new int[2*numberColumns_];
  CoinBigIndex j=0;
  for (iColumn=0;iColumn<numberColumns_;iColumn++, j+=2) {
    CoinBigIndex k=columnStart[iColumn];
    int iRow;
    switch (columnLength[iColumn]) {
    case 0:
      goodNetwork=-1; // not classic network
      indices_[j]=-1;
      indices_[j+1]=-1;
      break;
      
    case 1:
      goodNetwork=-1; // not classic network
      if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
	indices_[j] = -1;
	iRow = row[k];
	numberRows_ = CoinMax(numberRows_,iRow);
	indices_[j+1]=iRow;
      } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
	indices_[j+1] = -1;
	iRow = row[k];
	numberRows_ = CoinMax(numberRows_,iRow);
	indices_[j]=iRow;
      } else {
	goodNetwork = 0; // not a network
      }
      break;
      
    case 2:
      if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
	if (fabs(elementByColumn[k+1]+1.0)<1.0e-10) {
	  iRow = row[k];
	  numberRows_ = CoinMax(numberRows_,iRow);
	  indices_[j+1]=iRow;
	  iRow = row[k+1];
	  numberRows_ = CoinMax(numberRows_,iRow);
	  indices_[j] = iRow;
	} else {
	  goodNetwork = 0; // not a network
	}
      } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
	if (fabs(elementByColumn[k+1]-1.0)<1.0e-10) {
	  iRow = row[k];
	  numberRows_ = CoinMax(numberRows_,iRow);
	  indices_[j]=iRow;
	  iRow = row[k+1];
	  numberRows_ = CoinMax(numberRows_,iRow);
	  indices_[j+1] = iRow;
	} else {
	  goodNetwork = 0; // not a network
	}
      } else {
	goodNetwork = 0; // not a network
      }
      break;

    default:
      goodNetwork = 0; // not a network
      break;
    }
    if (!goodNetwork)
      break;
  }
  if (!goodNetwork) {
    delete [] indices_;
    // put in message
    printf("Not a network - can test if indices_ null\n");
    indices_=NULL;
    numberRows_=0;
    numberColumns_=0;
  } else {
    numberRows_ ++; //  correct
    trueNetwork_ = goodNetwork>0;
  }
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
ClpNetworkMatrix::~ClpNetworkMatrix ()
{
  delete matrix_;
  delete [] lengths_;
  delete [] indices_;
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpNetworkMatrix &
ClpNetworkMatrix::operator=(const ClpNetworkMatrix& rhs)
{
  if (this != &rhs) {
    ClpMatrixBase::operator=(rhs);
    delete matrix_;
    delete [] lengths_;
    delete [] indices_;
    matrix_ = NULL;
    lengths_=NULL;
    indices_=NULL;
    numberRows_=rhs.numberRows_;
    numberColumns_=rhs.numberColumns_;
    trueNetwork_=rhs.trueNetwork_;
    if (numberColumns_) {
      indices_ = new int [ 2*numberColumns_];
      memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
    }
  }
  return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpMatrixBase * ClpNetworkMatrix::clone() const
{
  return new ClpNetworkMatrix(*this);
}

/* Returns a new matrix in reverse order without gaps */
ClpMatrixBase * 
ClpNetworkMatrix::reverseOrderedCopy() const
{
  // count number in each row
  CoinBigIndex * tempP = new CoinBigIndex [numberRows_];
  CoinBigIndex * tempN = new CoinBigIndex [numberRows_];
  memset(tempP,0,numberRows_*sizeof(CoinBigIndex));
  memset(tempN,0,numberRows_*sizeof(CoinBigIndex));
  CoinBigIndex j=0;
  int i;
  for (i=0;i<numberColumns_;i++,j+=2) {
    int iRow = indices_[j];
    tempN[iRow]++;
    iRow = indices_[j+1];
    tempP[iRow]++;
  }
  int * newIndices = new int [2*numberColumns_];
  CoinBigIndex * newP = new CoinBigIndex [numberRows_+1];
  CoinBigIndex * newN = new CoinBigIndex[numberRows_];
  int iRow;
  j=0;
  // do starts
  for (iRow=0;iRow<numberRows_;iRow++) {
    newP[iRow]=j;
    j += tempP[iRow];
    tempP[iRow]=newP[iRow];
    newN[iRow] = j;
    j += tempN[iRow];
    tempN[iRow]=newN[iRow];
  }
  newP[numberRows_]=j;
  j=0;
  for (i=0;i<numberColumns_;i++,j+=2) {
    int iRow = indices_[j];
    CoinBigIndex put = tempN[iRow];
    newIndices[put++] = i;
    tempN[iRow] = put;
    iRow = indices_[j+1];
    put = tempP[iRow];
    newIndices[put++] = i;
    tempP[iRow] = put;
  }
  delete [] tempP;
  delete [] tempN;
  ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
  newCopy->passInCopy(numberRows_, numberColumns_,
		      false,  newIndices, newP, newN);
  return newCopy;
}
//unscaled versions
void 
ClpNetworkMatrix::times(double scalar,
		   const double * x, double * y) const
{
  int iColumn;
  CoinBigIndex j=0;
  if (trueNetwork_) {
    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
      double value = scalar*x[iColumn];
      if (value) {
	int iRowM = indices_[j];
	int iRowP = indices_[j+1];
	y[iRowM] -= value;
	y[iRowP] += value;
      }
    }
  } else {
    // skip negative rows
    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
      double value = scalar*x[iColumn];
      if (value) {
	int iRowM = indices_[j];
	int iRowP = indices_[j+1];
	if (iRowM>=0)
	  y[iRowM] -= value;
	if (iRowP>=0)
	  y[iRowP] += value;
      }
    }
  }
}
void 
ClpNetworkMatrix::transposeTimes(double scalar,
				const double * x, double * y) const
{
  int iColumn;
  CoinBigIndex j=0;
  if (trueNetwork_) {
    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
      double value = y[iColumn];
      int iRowM = indices_[j];
      int iRowP = indices_[j+1];
      value -= scalar*x[iRowM];
      value += scalar*x[iRowP];
      y[iColumn] = value;
    }
  } else {
    // skip negative rows
    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
      double value = y[iColumn];
      int iRowM = indices_[j];
      int iRowP = indices_[j+1];
      if (iRowM>=0)
	value -= scalar*x[iRowM];
      if (iRowP>=0)
	value += scalar*x[iRowP];
      y[iColumn] = value;
    }
  }
}
void 
ClpNetworkMatrix::times(double scalar,
		       const double * x, double * y,
		       const double * rowScale, 
		       const double * columnScale) const
{
  // we know it is not scaled 
  times(scalar, x, y);
}
void 
ClpNetworkMatrix::transposeTimes( double scalar,
				 const double * x, double * y,
				 const double * rowScale, 
				 const double * columnScale, double * spare) const
{
  // we know it is not scaled 
  transposeTimes(scalar, x, y);
}
/* Return <code>x * A + y</code> in <code>z</code>. 
	Squashes small elements and knows about ClpSimplex */
void 
ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar,
			      const CoinIndexedVector * rowArray,
			      CoinIndexedVector * y,
			      CoinIndexedVector * columnArray) const
{
  // we know it is not scaled 
  columnArray->clear();
  double * pi = rowArray->denseVector();
  int numberNonZero=0;
  int * index = columnArray->getIndices();
  double * array = columnArray->denseVector();
  int numberInRowArray = rowArray->getNumElements();
  // maybe I need one in OsiSimplex
  double zeroTolerance = model->factorization()->zeroTolerance();
  int numberRows = model->numberRows();
#ifndef NO_RTTI
  ClpPlusMinusOneMatrix* rowCopy =
    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
#else
  ClpPlusMinusOneMatrix* rowCopy =
    static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
#endif
  bool packed = rowArray->packedMode();
  double factor = 0.3;
  // We may not want to do by row if there may be cache problems
  int numberColumns = model->numberColumns();
  // It would be nice to find L2 cache size - for moment 512K
  // Be slightly optimistic
  if (numberColumns*sizeof(double)>1000000) {
    if (numberRows*10<numberColumns)
      factor=0.1;
    else if (numberRows*4<numberColumns)
      factor=0.15;
    else if (numberRows*2<numberColumns)
      factor=0.2;
    //if (model->numberIterations()%50==0)
    //printf("%d nonzero\n",numberInRowArray);
  }
  if (numberInRowArray>factor*numberRows||!rowCopy) {
    // do by column
    int iColumn;
    assert (!y->getNumElements());
    CoinBigIndex j=0;
    if (packed) {
      // need to expand pi into y
      assert(y->capacity()>=numberRows);
      double * piOld = pi;
      pi = y->denseVector();
      const int * whichRow = rowArray->getIndices();
      int i;
      // modify pi so can collapse to one loop
      for (i=0;i<numberInRowArray;i++) {
	int iRow = whichRow[i];
	pi[iRow]=scalar*piOld[i];
      }
      if (trueNetwork_) {
	for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
	  double value = 0.0;
	  int iRowM = indices_[j];
	  int iRowP = indices_[j+1];
	  value -= pi[iRowM];
	  value += pi[iRowP];
	  if (fabs(value)>zeroTolerance) {
	    array[numberNonZero]=value;
	    index[numberNonZero++]=iColumn;
	  }
	}
      } else {
	// skip negative rows
	for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
	  double value = 0.0;
	  int iRowM = indices_[j];
	  int iRowP = indices_[j+1];
	  if (iRowM>=0)
	    value -= pi[iRowM];
	  if (iRowP>=0)
	    value += pi[iRowP];
	  if (fabs(value)>zeroTolerance) {
	    array[numberNonZero]=value;
	    index[numberNonZero++]=iColumn;
	  }
	}
      }
      for (i=0;i<numberInRowArray;i++) {
	int iRow = whichRow[i];
	pi[iRow]=0.0;
      }
    } else {
      if (trueNetwork_) {
	for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
	  double value = 0.0;
	  int iRowM = indices_[j];
	  int iRowP = indices_[j+1];
	  value -= scalar*pi[iRowM];
	  value += scalar*pi[iRowP];
	  if (fabs(value)>zeroTolerance) {
	    index[numberNonZero++]=iColumn;
	    array[iColumn]=value;
	  }
	}
      } else {
	// skip negative rows
	for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
	  double value = 0.0;
	  int iRowM = indices_[j];
	  int iRowP = indices_[j+1];
	  if (iRowM>=0)
	    value -= scalar*pi[iRowM];
	  if (iRowP>=0)
	    value += scalar*pi[iRowP];
	  if (fabs(value)>zeroTolerance) {
	    index[numberNonZero++]=iColumn;
	    array[iColumn]=value;
	  }
	}
      }
    }
    columnArray->setNumElements(numberNonZero);
  } else {
    // do by row
    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
  }
}
/* Return <code>x *A in <code>z</code> but
   just for indices in y. */
void 
ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * model,
			      const CoinIndexedVector * rowArray,
			      const CoinIndexedVector * y,
			      CoinIndexedVector * columnArray) const
{
  columnArray->clear();
  double * pi = rowArray->denseVector();
  double * array = columnArray->denseVector();
  int jColumn;
  int numberToDo = y->getNumElements();
  const int * which = y->getIndices();
  assert (!rowArray->packedMode());
  columnArray->setPacked();
  if (trueNetwork_) {
    for (jColumn=0;jColumn<numberToDo;jColumn++) {
      int iColumn = which[jColumn];
      double value = 0.0;
      CoinBigIndex j=iColumn<<1;
      int iRowM = indices_[j];
      int iRowP = indices_[j+1];
      value -= pi[iRowM];
      value += pi[iRowP];
      array[jColumn]=value;
    }
  } else {
    // skip negative rows
    for (jColumn=0;jColumn<numberToDo;jColumn++) {
      int iColumn = which[jColumn];
      double value = 0.0;
      CoinBigIndex j=iColumn<<1;
      int iRowM = indices_[j];
      int iRowP = indices_[j+1];
      if (iRowM>=0)
	value -= pi[iRowM];
      if (iRowP>=0)
	value += pi[iRowP];
      array[jColumn]=value;
    }
  }
}
/// returns number of elements in column part of basis,
CoinBigIndex 
ClpNetworkMatrix::countBasis(ClpSimplex * model,
				 const int * whichColumn, 
				 int numberBasic,
			     int & numberColumnBasic)
{
  int i;
  CoinBigIndex numberElements=0;
  if (trueNetwork_) {
    numberElements = 2*numberColumnBasic;
  } else {
    for (i=0;i<numberColumnBasic;i++) {
      int iColumn = whichColumn[i];
      CoinBigIndex j=iColumn<<1;
      int iRowM = indices_[j];
      int iRowP = indices_[j+1];
      if (iRowM>=0)
	numberElements ++;
      if (iRowP>=0)
	numberElements ++;
    }
  }
  return numberElements;
}
void
ClpNetworkMatrix::fillBasis(ClpSimplex * model,
			 const int * whichColumn, 
			 int & numberColumnBasic,
			 int * indexRowU, int * start,
			 int * rowCount, int * columnCount,
			 double * elementU)
{
  int i;
  CoinBigIndex numberElements=start[0];
  if (trueNetwork_) {
    for (i=0;i<numberColumnBasic;i++) {
      int iColumn = whichColumn[i];
      CoinBigIndex j=iColumn<<1;
      int iRowM = indices_[j];
      int iRowP = indices_[j+1];
      indexRowU[numberElements]=iRowM;
      rowCount[iRowM]++;
      elementU[numberElements]=-1.0;
      indexRowU[numberElements+1]=iRowP;
      rowCount[iRowP]++;
      elementU[numberElements+1]=1.0;
      numberElements+=2;
      start[i+1]=numberElements;
      columnCount[i]=2;
    }
  } else {
    for (i=0;i<numberColumnBasic;i++) {
      int iColumn = whichColumn[i];
      CoinBigIndex j=iColumn<<1;
      int iRowM = indices_[j];
      int iRowP = indices_[j+1];
      if (iRowM>=0) {
	indexRowU[numberElements]=iRowM;
	rowCount[iRowM]++;
	elementU[numberElements++]=-1.0;
      }
      if (iRowP>=0) {
	indexRowU[numberElements]=iRowP;
	rowCount[iRowP]++;
	elementU[numberElements++]=1.0;
      }
      start[i+1]=numberElements;
      columnCount[i]=numberElements-start[i];
    }
  }
}
/* Unpacks a column into an CoinIndexedvector
 */
void 
ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
		   int iColumn) const 
{
  CoinBigIndex j=iColumn<<1;
  int iRowM = indices_[j];
  int iRowP = indices_[j+1];
  if (iRowM>=0) 
    rowArray->add(iRowM,-1.0);
  if (iRowP>=0) 
    rowArray->add(iRowP,1.0);
}
/* Unpacks a column into an CoinIndexedvector
** in packed foramt
Note that model is NOT const.  Bounds and objective could
be modified if doing column generation (just for this variable) */
void 
ClpNetworkMatrix::unpackPacked(ClpSimplex * model,
			    CoinIndexedVector * rowArray,
			    int iColumn) const
{
  int * index = rowArray->getIndices();
  double * array = rowArray->denseVector();
  int number = 0;
  CoinBigIndex j=iColumn<<1;
  int iRowM = indices_[j];
  int iRowP = indices_[j+1];
  if (iRowM>=0) {
    array[number]=-1.0;
    index[number++]=iRowM;
  }
  if (iRowP>=0) {
    array[number]=1.0;
    index[number++]=iRowP;
  }
  rowArray->setNumElements(number);
  rowArray->setPackedMode(true);
}
/* Adds multiple of a column into an CoinIndexedvector
      You can use quickAdd to add to vector */
void 
ClpNetworkMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
		   int iColumn, double multiplier) const 
{
  CoinBigIndex j=iColumn<<1;
  int iRowM = indices_[j];
  int iRowP = indices_[j+1];
  if (iRowM>=0) 
    rowArray->quickAdd(iRowM,-multiplier);
  if (iRowP>=0) 
    rowArray->quickAdd(iRowP,multiplier);
}
/* Adds multiple of a column into an array */
void 
ClpNetworkMatrix::add(const ClpSimplex * model,double * array,
		    int iColumn, double multiplier) const
{
  CoinBigIndex j=iColumn<<1;
  int iRowM = indices_[j];
  int iRowP = indices_[j+1];
  if (iRowM>=0) 
    array[iRowM] -= multiplier;
  if (iRowP>=0) 
    array[iRowP] += multiplier;
}

// Return a complete CoinPackedMatrix
CoinPackedMatrix * 
ClpNetworkMatrix::getPackedMatrix() const 
{
  if (!matrix_) {
    assert (trueNetwork_); // fix later
    int numberElements = 2*numberColumns_;
    double * elements = new double [numberElements];
    CoinBigIndex i;
    for (i=0;i<2*numberColumns_;i+=2) {
      elements[i]=-1.0;
      elements[i+1]=1.0;
    }
    CoinBigIndex * starts = new CoinBigIndex [numberColumns_+1];
    for (i=0;i<numberColumns_+1;i++) {
      starts[i]=2*i;
    }
    // use assignMatrix to save space
    delete [] lengths_;
    lengths_ = NULL;
    matrix_ = new CoinPackedMatrix();
    int * indices = CoinCopyOfArray(indices_,2*numberColumns_);
    matrix_->assignMatrix(true,numberRows_,numberColumns_,
			  getNumElements(),
			  elements,indices,
			  starts,lengths_);
    assert(!elements);
    assert(!starts);
    assert (!indices);
    assert (!lengths_);
  }
  return matrix_;
}
/* A vector containing the elements in the packed matrix. Note that there
   might be gaps in this list, entries that do not belong to any
   major-dimension vector. To get the actual elements one should look at
   this vector together with vectorStarts and vectorLengths. */
const double * 
ClpNetworkMatrix::getElements() const 
{
  if (!matrix_) 
    getPackedMatrix();
  return matrix_->getElements();
}

const CoinBigIndex * 
ClpNetworkMatrix::getVectorStarts() const 
{
  if (!matrix_) 
    getPackedMatrix();
  return matrix_->getVectorStarts();
}
/* The lengths of the major-dimension vectors. */
const int * 
ClpNetworkMatrix::getVectorLengths() const
{
  assert (trueNetwork_); // fix later
  if (!lengths_) {
    lengths_ = new int [numberColumns_];
    int i;
    for (i=0;i<numberColumns_;i++) {
      lengths_[i]=2;
    }
  }
  return lengths_;
}
/* Delete the columns whose indices are listed in <code>indDel</code>. */
void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel) 
{
  assert (trueNetwork_);
  int iColumn;
  int numberBad=0;
  // Use array to make sure we can have duplicates
  char * which = new char[numberColumns_];
  memset(which,0,numberColumns_);
  int nDuplicate=0;
  for (iColumn=0;iColumn<numDel;iColumn++) {
    int jColumn = indDel[iColumn];
    if (jColumn<0||jColumn>=numberColumns_) {
      numberBad++;
    } else {
      if (which[jColumn])
	nDuplicate++;
      else
	which[jColumn]=1;
    }
  }
  if (numberBad)
    throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix");
  int newNumber = numberColumns_-numDel+nDuplicate;
  // Get rid of temporary arrays
  delete [] lengths_;
  lengths_=NULL;
  delete matrix_;
  matrix_= NULL;
  int newSize=2*newNumber;
  int * newIndices = new int [newSize];
  newSize=0;
  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    if (!which[iColumn]) {
      CoinBigIndex start;
      CoinBigIndex i;
      start = 2*iColumn;
      for (i=start;i<start+2;i++) 
	newIndices[newSize++]=indices_[i];
    }
  }
  delete [] which;
  delete [] indices_;
  indices_= newIndices;
  numberColumns_ = newNumber;
}
/* Delete the rows whose indices are listed in <code>indDel</code>. */
void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel) 
{
  int iRow;
  int numberBad=0;
  // Use array to make sure we can have duplicates
  int * which = new int [numberRows_];
  memset(which,0,numberRows_*sizeof(int));
  for (iRow=0;iRow<numDel;iRow++) {
    int jRow = indDel[iRow];
    if (jRow<0||jRow>=numberRows_) {
      numberBad++;
    } else {
      which[jRow]=1;
    }
  }
  if (numberBad)
    throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix");
  // Only valid of all columns have 0 entries
  int iColumn;
  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    CoinBigIndex start;
    CoinBigIndex i;
    start = 2*iColumn;
    for (i=start;i<start+2;i++) {
      int iRow = indices_[i];
      if (which[iRow])
	numberBad++;
    }
  }
  if (numberBad)
    throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix");
  int newNumber=0;
  for (iRow=0;iRow<numberRows_;iRow++) {
    if (!which[iRow])
      which[iRow]=newNumber++;
    else
      which[iRow]=-1;
  }
  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    CoinBigIndex start;
    CoinBigIndex i;
    start = 2*iColumn;
    for (i=start;i<start+2;i++) {
      int iRow = indices_[i];
      indices_[i]=which[iRow];
    }
  }
  delete [] which;
  numberRows_ = newNumber;
}
/* Given positive integer weights for each row fills in sum of weights
   for each column (and slack).
   Returns weights vector
*/
CoinBigIndex * 
ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
{
  int numberRows = model->numberRows();
  int numberColumns =model->numberColumns();
  int number = numberRows+numberColumns;
  CoinBigIndex * weights = new CoinBigIndex[number];
  int i;
  for (i=0;i<numberColumns;i++) {
    CoinBigIndex j=i<<1;
    CoinBigIndex count=0;
    int iRowM = indices_[j];
    int iRowP = indices_[j+1];
    if (iRowM>=0) {
      count += inputWeights[iRowM];
    }
    if (iRowP>=0) {
      count += inputWeights[iRowP];
    }
    weights[i]=count;
  }
  for (i=0;i<numberRows;i++) {
    weights[i+numberColumns]=inputWeights[i];
  }
  return weights;
}
/* Returns largest and smallest elements of both signs.
   Largest refers to largest absolute value.
*/
void 
ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
		       double & smallestPositive, double & largestPositive)
{
  smallestNegative=-1.0;
  largestNegative=-1.0;
  smallestPositive=1.0;
  largestPositive=1.0;
}
// Says whether it can do partial pricing
bool 
ClpNetworkMatrix::canDoPartialPricing() const
{
  return true; 
}
// Partial pricing 
void 
ClpNetworkMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
			      int & bestSequence, int & numberWanted)
{
  numberWanted=currentWanted_;
  int j;
  int start = (int) (startFraction*numberColumns_);
  int end = CoinMin((int) (endFraction*numberColumns_+1),numberColumns_);
  double tolerance=model->currentDualTolerance();
  double * reducedCost = model->djRegion();
  const double * duals = model->dualRowSolution();
  const double * cost = model->costRegion();
  double bestDj;
  if (bestSequence>=0)
    bestDj = fabs(reducedCost[bestSequence]);
  else
    bestDj=tolerance;
  int sequenceOut = model->sequenceOut();
  int saveSequence = bestSequence;
  if (!trueNetwork_) {
    // Not true network
    int iSequence;
    for (iSequence=start;iSequence<end;iSequence++) {
      if (iSequence!=sequenceOut) {
	double value;
	int iRowM,iRowP;
	ClpSimplex::Status status = model->getStatus(iSequence);
	
	switch(status) {
	  
	case ClpSimplex::basic:
	case ClpSimplex::isFixed:
	  break;
	case ClpSimplex::isFree:
	case ClpSimplex::superBasic:
	  value=cost[iSequence];
	  j = iSequence<<1;
	  // skip negative rows
	  iRowM = indices_[j];
	  iRowP = indices_[j+1];
	  if (iRowM>=0)
	    value += duals[iRowM];
	  if (iRowP>=0)
	    value -= duals[iRowP];
	  value = fabs(value);
	  if (value>FREE_ACCEPT*tolerance) {
	    numberWanted--;
	    // we are going to bias towards free (but only if reasonable)
	    value *= FREE_BIAS;
	    if (value>bestDj) {
	      // check flagged variable and correct dj
	      if (!model->flagged(iSequence)) {
		bestDj=value;
		bestSequence = iSequence;
	      } else {
		// just to make sure we don't exit before got something
		numberWanted++;
	      }
	    }
	  }
	  break;
	case ClpSimplex::atUpperBound:
	  value=cost[iSequence];
	  j = iSequence<<1;
	  // skip negative rows
	  iRowM = indices_[j];
	  iRowP = indices_[j+1];
	  if (iRowM>=0)
	    value += duals[iRowM];
	  if (iRowP>=0)
	    value -= duals[iRowP];
	  if (value>tolerance) {
	    numberWanted--;
	    if (value>bestDj) {
	      // check flagged variable and correct dj
	      if (!model->flagged(iSequence)) {
		bestDj=value;
		bestSequence = iSequence;
	      } else {
		// just to make sure we don't exit before got something
		numberWanted++;
	      }
	    }
	  }
	  break;
	case ClpSimplex::atLowerBound:
	  value=cost[iSequence];
	  j = iSequence<<1;
	  // skip negative rows
	  iRowM = indices_[j];
	  iRowP = indices_[j+1];
	  if (iRowM>=0)
	    value += duals[iRowM];
	  if (iRowP>=0)
	    value -= duals[iRowP];
	  value = -value;
	  if (value>tolerance) {
	    numberWanted--;
	    if (value>bestDj) {
	      // check flagged variable and correct dj
	      if (!model->flagged(iSequence)) {
		bestDj=value;
		bestSequence = iSequence;
	      } else {
		// just to make sure we don't exit before got something
		numberWanted++;
	      }
	    }
	  }
	  break;
	}
      }
      if (!numberWanted)
	break;
    }
    if (bestSequence!=saveSequence) {
      // recompute dj
      double value=cost[bestSequence];
      j = bestSequence<<1;
      // skip negative rows
      int iRowM = indices_[j];
      int iRowP = indices_[j+1];
      if (iRowM>=0)
	value += duals[iRowM];
      if (iRowP>=0)
	value -= duals[iRowP];
      reducedCost[bestSequence] = value;
      savedBestSequence_ = bestSequence;
      savedBestDj_ = reducedCost[savedBestSequence_];
    }
  } else {
    // true network
    int iSequence;
    for (iSequence=start;iSequence<end;iSequence++) {
      if (iSequence!=sequenceOut) {
	double value;
	int iRowM,iRowP;
	ClpSimplex::Status status = model->getStatus(iSequence);
	
	switch(status) {
	  
	case ClpSimplex::basic:
	case ClpSimplex::isFixed:
	  break;
	case ClpSimplex::isFree:
	case ClpSimplex::superBasic:
	  value=cost[iSequence];
	  j = iSequence<<1;
	  iRowM = indices_[j];
	  iRowP = indices_[j+1];
	  value += duals[iRowM];
	  value -= duals[iRowP];
	  value = fabs(value);
	  if (value>FREE_ACCEPT*tolerance) {
	    numberWanted--;
	    // we are going to bias towards free (but only if reasonable)
	    value *= FREE_BIAS;
	    if (value>bestDj) {
	      // check flagged variable and correct dj
	      if (!model->flagged(iSequence)) {
		bestDj=value;
		bestSequence = iSequence;
	      } else {
		// just to make sure we don't exit before got something
		numberWanted++;
	      }
	    }
	  }
	  break;
	case ClpSimplex::atUpperBound:
	  value=cost[iSequence];
	  j = iSequence<<1;
	  iRowM = indices_[j];
	  iRowP = indices_[j+1];
	  value += duals[iRowM];
	  value -= duals[iRowP];
	  if (value>tolerance) {
	    numberWanted--;
	    if (value>bestDj) {
	      // check flagged variable and correct dj
	      if (!model->flagged(iSequence)) {
		bestDj=value;
		bestSequence = iSequence;
	      } else {
		// just to make sure we don't exit before got something
		numberWanted++;
	      }
	    }
	  }
	  break;
	case ClpSimplex::atLowerBound:
	  value=cost[iSequence];
	  j = iSequence<<1;
	  iRowM = indices_[j];
	  iRowP = indices_[j+1];
	  value += duals[iRowM];
	  value -= duals[iRowP];
	  value = -value;
	  if (value>tolerance) {
	    numberWanted--;
	    if (value>bestDj) {
	      // check flagged variable and correct dj
	      if (!model->flagged(iSequence)) {
		bestDj=value;
		bestSequence = iSequence;
	      } else {
		// just to make sure we don't exit before got something
		numberWanted++;
	      }
	    }
	  }
	  break;
	}
      }
      if (!numberWanted)
	break;
    }
    if (bestSequence!=saveSequence) {
      // recompute dj
      double value=cost[bestSequence];
      j = bestSequence<<1;
      int iRowM = indices_[j];
      int iRowP = indices_[j+1];
      value += duals[iRowM];
      value -= duals[iRowP];
      reducedCost[bestSequence] = value;
      savedBestSequence_ = bestSequence;
      savedBestDj_ = reducedCost[savedBestSequence_];
    }
  }
  currentWanted_=numberWanted;
}
// Allow any parts of a created CoinMatrix to be deleted
void 
ClpNetworkMatrix::releasePackedMatrix() const 
{
  delete matrix_;
  delete [] lengths_;
  matrix_=NULL;
  lengths_=NULL;
}
// Append Columns
void 
ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
{
  int iColumn;
  int numberBad=0;
  for (iColumn=0;iColumn<number;iColumn++) {
    int n=columns[iColumn]->getNumElements();
    const double * element = columns[iColumn]->getElements();
    if (n!=2)
      numberBad++;
    if (fabs(element[0])!=1.0||fabs(element[1])!=1.0)
	numberBad++;
    else if (element[0]*element[1]!=-1.0)
	numberBad++;
  }
  if (numberBad)
    throw CoinError("Not network", "appendCols", "ClpNetworkMatrix");
  // Get rid of temporary arrays
  delete [] lengths_;
  lengths_=NULL;
  delete matrix_;
  matrix_= NULL;
  CoinBigIndex size = 2*number;
  int * temp2 = new int [numberColumns_*2+size];
  memcpy(temp2,indices_,numberColumns_*2*sizeof(int));
  delete [] indices_;
  indices_= temp2;
  // now add
  size=2*numberColumns_;
  for (iColumn=0;iColumn<number;iColumn++) {
    const int * row = columns[iColumn]->getIndices();
    const double * element = columns[iColumn]->getElements();
    if (element[0]==-1.0) {
      indices_[size++]=row[0];
      indices_[size++]=row[1];
    } else {
      indices_[size++]=row[1];
      indices_[size++]=row[0];
    }
  }
  
  numberColumns_ += number;
}
// Append Rows
void 
ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
{
  // must be zero arrays
  int numberBad=0;
  int iRow;
  for (iRow=0;iRow<number;iRow++) {
    numberBad += rows[iRow]->getNumElements();
  }
  if (numberBad)
    throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix");
  numberRows_ += number;
}
#ifndef SLIM_CLP
/* Append a set of rows/columns to the end of the matrix. Returns number of errors
   i.e. if any of the new rows/columns contain an index that's larger than the
   number of columns-1/rows-1 (if numberOther>0) or duplicates
   If 0 then rows, 1 if columns */
int 
ClpNetworkMatrix::appendMatrix(int number, int type,
                                    const CoinBigIndex * starts, const int * index,
                                    const double * element, int numberOther)
{
  int numberErrors=0;
  // make into CoinPackedVector
  CoinPackedVectorBase ** vectors=
    new CoinPackedVectorBase * [number];
  int iVector;
  for (iVector=0;iVector<number;iVector++) {
    int iStart = starts[iVector];
    vectors[iVector] = 
      new CoinPackedVector(starts[iVector+1]-iStart,
                           index+iStart,element+iStart);
  }
  if (type==0) {
    // rows
    appendRows(number,vectors);
  } else {
    // columns
    appendCols(number,vectors);
  }
  for (iVector=0;iVector<number;iVector++) 
    delete vectors[iVector];
  delete [] vectors;
  return numberErrors;
}
#endif
/* Subset clone (without gaps).  Duplicates are allowed
   and order is as given */
ClpMatrixBase * 
ClpNetworkMatrix::subsetClone (int numberRows, const int * whichRows,
			      int numberColumns, 
			      const int * whichColumns) const 
{
  return new ClpNetworkMatrix(*this, numberRows, whichRows,
				   numberColumns, whichColumns);
}
/* Subset constructor (without gaps).  Duplicates are allowed
   and order is as given */
ClpNetworkMatrix::ClpNetworkMatrix (
		       const ClpNetworkMatrix & rhs,
		       int numberRows, const int * whichRow,
		       int numberColumns, const int * whichColumn)
: ClpMatrixBase(rhs)
{  
  setType(11);
  matrix_ = NULL;
  lengths_=NULL;
  indices_=new int[2*numberColumns];;
  numberRows_=numberRows;
  numberColumns_=numberColumns;
  trueNetwork_=true;
  int iColumn;
  int numberBad=0;
  int * which = new int [rhs.numberRows_];
  int iRow;
  for (iRow=0;iRow<rhs.numberRows_;iRow++) 
    which[iRow]=-1;
  int n=0;
  for (iRow=0;iRow<numberRows;iRow++) {
    int jRow=whichRow[iRow];
    assert (jRow>=0&&jRow<rhs.numberRows_);
    which[jRow]=n++;
  }
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    CoinBigIndex start;
    CoinBigIndex i;
    start = 2*iColumn;
    CoinBigIndex offset = 2*whichColumn[iColumn]-start;
    for (i=start;i<start+2;i++) {
      int iRow = rhs.indices_[i+offset];
      iRow = which[iRow];
      if (iRow<0)
	numberBad++;
      else
	indices_[i]=iRow;
    }
  }
  if (numberBad)
    throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix");
}