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

#if defined(_MSC_VER)
// Turn off compiler warning about long names
#  pragma warning(disable:4786)
#endif

#include "CoinUtilsConfig.h"

#include <algorithm>
#include <numeric>
#include <cassert>
#include <cstdio>
#include <iostream>

#include "CoinSort.hpp"
#include "CoinHelperFunctions.hpp"
#ifndef CLP_NO_VECTOR
#include "CoinPackedVectorBase.hpp"
#endif
#include "CoinPackedMatrix.hpp"
#include <stdio.h>

//#############################################################################
static inline int
CoinLengthWithExtra(int len, double extraGap)
{
  return static_cast<int>(ceil(len * (1 + extraGap)));
}

//#############################################################################

static inline void
CoinTestSortedIndexSet(const int num, const int * sorted, const int maxEntry,
		      const char * testingMethod)
{
   if (sorted[0] < 0 || sorted[num-1] >= maxEntry)
      throw CoinError("bad index", testingMethod, "CoinPackedMatrix");
   if (std::adjacent_find(sorted, sorted + num) != sorted + num)
      throw CoinError("duplicate index", testingMethod, "CoinPackedMatrix");
}

//-----------------------------------------------------------------------------

static inline int *
CoinTestIndexSet(const int numDel, const int * indDel, const int maxEntry,
		const char * testingMethod)
{
   if (! CoinIsSorted(indDel, indDel + numDel)) {
      // if not sorted then sort it, test for consistency and return a pointer
      // to the sorted array
      int * sorted = new int[numDel];
      CoinMemcpyN(indDel, numDel, sorted);
      std::sort(sorted, sorted + numDel);
      CoinTestSortedIndexSet(numDel, sorted, maxEntry, testingMethod);
      return sorted;
   }

   // Otherwise it's already sorted, so just test for consistency and return a
   // 0 pointer.
   CoinTestSortedIndexSet(numDel, indDel, maxEntry, testingMethod);
   return 0;
}

//#############################################################################

void
CoinPackedMatrix::reserve(const int newMaxMajorDim, const CoinBigIndex newMaxSize,
			  bool create)
{
   if (newMaxMajorDim > maxMajorDim_) {
      maxMajorDim_ = newMaxMajorDim;
      int * oldlength = length_;
      CoinBigIndex * oldstart = start_;
      length_ = new int[newMaxMajorDim];
      start_ = new CoinBigIndex[newMaxMajorDim+1];
      start_[0]=0;
      if (majorDim_ > 0) {
	 CoinMemcpyN(oldlength, majorDim_, length_);
	 CoinMemcpyN(oldstart, majorDim_ + 1, start_);
      }
      if (create) {
	// create empty vectors
	CoinFillN(length_+majorDim_,maxMajorDim_-majorDim_,0);
	CoinFillN(start_+majorDim_+1,maxMajorDim_-majorDim_,0);
	majorDim_=maxMajorDim_;
      }
      delete[] oldlength;
      delete[] oldstart;
   }
   if (newMaxSize > maxSize_) {
      maxSize_ = newMaxSize;
      int * oldind = index_;
      double * oldelem = element_;
      index_ = new int[newMaxSize];
      element_ = new double[newMaxSize];
      for (int i = majorDim_ - 1; i >= 0; --i) {
	 CoinMemcpyN(oldind+start_[i], length_[i], index_+start_[i]);
	 CoinMemcpyN(oldelem+start_[i], length_[i], element_+start_[i]);
      }
      delete[] oldind;
      delete[] oldelem;
   }
}

//-----------------------------------------------------------------------------

void
CoinPackedMatrix::clear()
{
   majorDim_ = 0;
   minorDim_ = 0;
   size_ = 0;
}

//#############################################################################
//#############################################################################

void
CoinPackedMatrix::setDimensions(int newnumrows, int newnumcols)
{
  const int numrows = getNumRows();
  if (newnumrows < 0)
    newnumrows = numrows;
  if (newnumrows < numrows)
    throw CoinError("Bad new rownum (less than current)",
		    "setDimensions", "CoinPackedMatrix");

  const int numcols = getNumCols();
  if (newnumcols < 0)
    newnumcols = numcols;
  if (newnumcols < numcols)
    throw CoinError("Bad new colnum (less than current)",
		    "setDimensions", "CoinPackedMatrix");

  int numplus = 0;
  if (isColOrdered()) {
    minorDim_ = newnumrows;
    numplus = newnumcols - numcols;
  } else {
    minorDim_ = newnumcols;
    numplus = newnumrows - numrows;
  }
  if (numplus > 0) {
    int* lengths = new int[numplus];
    CoinZeroN(lengths, numplus);
    resizeForAddingMajorVectors(numplus, lengths);
    delete[] lengths;
    majorDim_ += numplus; //forgot to change majorDim_
  }

}

//-----------------------------------------------------------------------------

void
CoinPackedMatrix::setExtraGap(const double newGap)
{
   if (newGap < 0)
      throw CoinError("negative new extra gap",
		     "setExtraGap", "CoinPackedMatrix");
   extraGap_ = newGap;
}

//-----------------------------------------------------------------------------

void
CoinPackedMatrix::setExtraMajor(const double newMajor)
{
   if (newMajor < 0)
      throw CoinError("negative new extra major",
		     "setExtraMajor", "CoinPackedMatrix");
   extraMajor_ = newMajor;
}

//#############################################################################
#ifndef CLP_NO_VECTOR
void
CoinPackedMatrix::appendCol(const CoinPackedVectorBase& vec)
{
   if (colOrdered_)
      appendMajorVector(vec);
   else
      appendMinorVector(vec);
}
#endif
//-----------------------------------------------------------------------------

void
CoinPackedMatrix::appendCol(const int vecsize,
			   const int *vecind,
			   const double *vecelem)
{
   if (colOrdered_)
      appendMajorVector(vecsize, vecind, vecelem);
   else
      appendMinorVector(vecsize, vecind, vecelem);
}

//-----------------------------------------------------------------------------
#ifndef CLP_NO_VECTOR
void
CoinPackedMatrix::appendCols(const int numcols,
			    const CoinPackedVectorBase * const * cols)
{
   if (colOrdered_)
      appendMajorVectors(numcols, cols);
   else
      appendMinorVectors(numcols, cols);
}
#endif
//-----------------------------------------------------------------------------

int 
CoinPackedMatrix::appendCols(const int numcols,
                             const CoinBigIndex * columnStarts, const int * row,
                             const double * element, int numberRows)
{
  int numberErrors;
  if (colOrdered_) {
    numberErrors=appendMajor(numcols, columnStarts, row, element, numberRows);
  } else {
    numberErrors=appendMinor(numcols, columnStarts, row, element, numberRows);
  }
  return numberErrors;
}
//-----------------------------------------------------------------------------
#ifndef CLP_NO_VECTOR
void
CoinPackedMatrix::appendRow(const CoinPackedVectorBase& vec)
{
   if (colOrdered_)
      appendMinorVector(vec);
   else
      appendMajorVector(vec);
}
#endif
//-----------------------------------------------------------------------------

void
CoinPackedMatrix::appendRow(const int vecsize,
			   const int *vecind,
			   const double *vecelem)
{
   if (colOrdered_)
      appendMinorVector(vecsize, vecind, vecelem);
   else
      appendMajorVector(vecsize, vecind, vecelem);
}

//-----------------------------------------------------------------------------
#ifndef CLP_NO_VECTOR
void
CoinPackedMatrix::appendRows(const int numrows,
			    const CoinPackedVectorBase * const * rows)
{
  if (colOrdered_) {
    // make sure enough columns
    if (numrows == 0)
      return;

    int i;
    int maxDim=-1;
    for (i = numrows - 1; i >= 0; --i) {
      const int vecsize = rows[i]->getNumElements();
      const int* vecind = rows[i]->getIndices();
      for (int j = vecsize - 1; j >= 0; --j) 
	maxDim = CoinMax(maxDim,vecind[j]);
    }
    maxDim++;
    if (maxDim>majorDim_) {
      setDimensions(minorDim_,maxDim);
      //int nAdd=maxDim-majorDim_;
      //int * length = new int[nAdd];
      //memset(length,0,nAdd*sizeof(int));
      //resizeForAddingMajorVectors(nAdd,length);
      //delete [] length;
    }
    appendMinorVectors(numrows, rows);
  } else {
    appendMajorVectors(numrows, rows);
  }
}
#endif
//-----------------------------------------------------------------------------

int 
CoinPackedMatrix::appendRows(const int numrows,
                             const CoinBigIndex * rowStarts, const int * column,
                             const double * element, int numberColumns)
{
  int numberErrors;
  if (colOrdered_) {
    numberErrors=appendMinor(numrows, rowStarts, column, element, numberColumns);
  } else {
    numberErrors=appendMajor(numrows, rowStarts, column, element, numberColumns);
  }
  return numberErrors;
}

//#############################################################################

void
CoinPackedMatrix::rightAppendPackedMatrix(const CoinPackedMatrix& matrix)
{
   if (colOrdered_) {
      if (matrix.colOrdered_) {
	 majorAppendSameOrdered(matrix);
      } else {
	 majorAppendOrthoOrdered(matrix);
      }
   } else {
      if (matrix.colOrdered_) {
	 minorAppendOrthoOrdered(matrix);
      } else {
	 minorAppendSameOrdered(matrix);
      }
   }
}

//-----------------------------------------------------------------------------

void
CoinPackedMatrix::bottomAppendPackedMatrix(const CoinPackedMatrix& matrix)
{
   if (colOrdered_) {
      if (matrix.colOrdered_) {
	 minorAppendSameOrdered(matrix);
      } else {
	 minorAppendOrthoOrdered(matrix);
      }
   } else {
      if (matrix.colOrdered_) {
	 majorAppendOrthoOrdered(matrix);
      } else {
	 majorAppendSameOrdered(matrix);
      }
   }
}

//#############################################################################

void
CoinPackedMatrix::deleteCols(const int numDel, const int * indDel)
{
  if (numDel) {
    if (colOrdered_)
      deleteMajorVectors(numDel, indDel);
    else
      deleteMinorVectors(numDel, indDel);
  }
}

//-----------------------------------------------------------------------------

void
CoinPackedMatrix::deleteRows(const int numDel, const int * indDel)
{
  if (numDel) {
    if (colOrdered_)
      deleteMinorVectors(numDel, indDel);
    else
      deleteMajorVectors(numDel, indDel);
  }
}

//#############################################################################
/* Replace the elements of a vector.  The indices remain the same.
   At most the number specified will be replaced.
   The index is between 0 and major dimension of matrix */
void 
CoinPackedMatrix::replaceVector(const int index,
			       const int numReplace, 
			       const double * newElements)
{
  if (index >= 0 && index < majorDim_) {
    int length = (length_[index] < numReplace) ? length_[index] : numReplace;
    CoinMemcpyN(newElements, length, element_ + start_[index]);
  } else {
#ifdef COIN_DEBUG
    throw CoinError("bad index", "replaceVector", "CoinPackedMatrix");
#endif
  }
}
/* Modify one element of packed matrix.  An element may be added.
   If the new element is zero it will be deleted unless
   keepZero true */
void 
CoinPackedMatrix::modifyCoefficient(int row, int column, double newElement,
				    bool keepZero)
{
  int minorIndex,majorIndex;
  if (colOrdered_) {
    majorIndex=column;
    minorIndex=row;
  } else {
    minorIndex=column;
    majorIndex=row;
  }
  if (majorIndex >= 0 && majorIndex < majorDim_) {
    if (minorIndex >= 0 && minorIndex < minorDim_) {
      CoinBigIndex j;
      CoinBigIndex end=start_[majorIndex]+length_[majorIndex];;
      for (j=start_[majorIndex];j<end;j++) {
	if (minorIndex==index_[j]) {
	  // replacement
	  if (newElement||keepZero) {
	    element_[j]=newElement;
	  } else {
	    // pack down and return
	    length_[majorIndex]--;
	    end--;
	    size_--;
	    for (;j<end;j++) {
	      element_[j]=element_[j+1];
	      index_[j]=index_[j+1];
	    }
	  }
	  return;
	}
      }
      if (j==end&&(newElement||keepZero)) {
	// we need to insert - keep in minor order if possible
	if (end>=start_[majorIndex+1]) {
	   int * addedEntries = new int[majorDim_];
	   memset(addedEntries, 0, majorDim_ * sizeof(int));
	   addedEntries[majorIndex] = 1;
	   resizeForAddingMinorVectors(addedEntries);
	   delete[] addedEntries;
        }
        // So where to insert? We're just going to assume that the entries
        // in the major vector are in increasing order, so we'll insert the
        // new entry to the last place we can
        const CoinBigIndex start = start_[majorIndex];
        end = start_[majorIndex]+length_[majorIndex]; // recalculate end
        for (j = end - 1; j >= start; --j) {
          if (index_[j] < minorIndex)
            break;
          index_[j+1] = index_[j];
          element_[j+1] = element_[j];
        }
        ++j;
        index_[j] = minorIndex;
        element_[j] = newElement;
        size_++;
        length_[majorIndex]++;
      }
    } else {
#ifdef COIN_DEBUG
      throw CoinError("bad minor index", "modifyCoefficient",
		      "CoinPackedMatrix");
#endif
    }
  } else {
#ifdef COIN_DEBUG
    throw CoinError("bad major index", "modifyCoefficient",
		    "CoinPackedMatrix");
#endif
  }
}
/* Return one element of packed matrix.
   This works for either ordering
   If it is not present will return 0.0 */
double 
CoinPackedMatrix::getCoefficient(int row, int column) const
{
  int minorIndex,majorIndex;
  if (colOrdered_) {
    majorIndex=column;
    minorIndex=row;
  } else {
    minorIndex=column;
    majorIndex=row;
  }
  double value=0.0;
  if (majorIndex >= 0 && majorIndex < majorDim_) {
    if (minorIndex >= 0 && minorIndex < minorDim_) {
      CoinBigIndex j;
      CoinBigIndex end=start_[majorIndex]+length_[majorIndex];;
      for (j=start_[majorIndex];j<end;j++) {
	if (minorIndex==index_[j]) {
          value = element_[j];
          break;
	}
      }
    } else {
#ifdef COIN_DEBUG
      throw CoinError("bad minor index", "modifyCoefficient",
		      "CoinPackedMatrix");
#endif
    }
  } else {
#ifdef COIN_DEBUG
    throw CoinError("bad major index", "modifyCoefficient",
		    "CoinPackedMatrix");
#endif
  }
  return value;
}

//#############################################################################
/* Eliminate all elements in matrix whose 
   absolute value is less than threshold.
   The column starts are not affected.  Returns number of elements
   eliminated.  Elements eliminated are at end of each vector
*/
int 
CoinPackedMatrix::compress(double threshold)
{
  CoinBigIndex numberEliminated =0;
  // space for eliminated
  int * eliminatedIndex = new int[minorDim_];
  double * eliminatedElement = new double[minorDim_];
  int i;
  for (i=0;i<majorDim_;i++) {
    int length = length_[i];
    CoinBigIndex k=start_[i];
    int kbad=0;
    CoinBigIndex j;
    for (j=start_[i];j<start_[i]+length;j++) {
      if (fabs(element_[j])>=threshold) {
	element_[k]=element_[j];
	index_[k++]=index_[j];
      } else {
	eliminatedElement[kbad]=element_[j];
	eliminatedIndex[kbad++]=index_[j];
      }
    }
    if (kbad) {
      numberEliminated += kbad;
      length_[i] = k-start_[i];
      memcpy(index_+k,eliminatedIndex,kbad*sizeof(int));
      memcpy(element_+k,eliminatedElement,kbad*sizeof(double));
    }
  }
  size_ -= numberEliminated;
  delete [] eliminatedIndex;
  delete [] eliminatedElement;
  return numberEliminated;
}
//#############################################################################
/* Eliminate all elements in matrix whose 
   absolute value is less than threshold.ALSO removes duplicates
   The column starts are not affected.  Returns number of elements
   eliminated. 
*/
int 
CoinPackedMatrix::eliminateDuplicates(double threshold)
{
  CoinBigIndex numberEliminated =0;
  // space for eliminated
  int * mark = new int [minorDim_];
  int i;
  for (i=0;i<minorDim_;i++)
    mark[i]=-1;
  for (i=0;i<majorDim_;i++) {
    CoinBigIndex k=start_[i];
    CoinBigIndex end = k+length_[i];
    CoinBigIndex j;
    for (j=k;j<end;j++) {
      int index = index_[j];
      if (mark[index]==-1) {
	mark[index]=j;
      } else {
	// duplicate
	int jj = mark[index];
	element_[jj] += element_[j];
	element_[j]=0.0;
      }
    }
    for (j=k;j<end;j++) {
      int index = index_[j];
      mark[index]=-1;
      if (fabs(element_[j])>=threshold) {
	element_[k]=element_[j];
	index_[k++]=index_[j];
      }
    }
    numberEliminated += end-k;
    length_[i] = k-start_[i];
  }
  size_ -= numberEliminated;
  delete [] mark;
  return numberEliminated;
}
//#############################################################################

void
CoinPackedMatrix::removeGaps(double removeValue)
{
  if (removeValue<0.0) {
    if (extraGap_) {
      for (int i = 1; i < majorDim_; ++i) {
	const CoinBigIndex si = start_[i];
	const int li = length_[i];
	start_[i] = start_[i-1] + length_[i-1];
	CoinCopy(index_ + si, index_ + (si + li), index_ + start_[i]);
	CoinCopy(element_ + si, element_ + (si + li), element_ + start_[i]);
      }
      start_[majorDim_] = size_;
    } else {
#ifndef NDEBUG
      for (int i = 1; i < majorDim_; ++i) {
	assert (start_[i] == start_[i-1] + length_[i-1]);
      }
      assert(start_[majorDim_] == size_);
#endif
    }
  } else {
    CoinBigIndex put=0;
    assert (!start_[0]);
    CoinBigIndex start = 0;
    for (int i = 0; i < majorDim_; ++i) {
      const CoinBigIndex si = start;
      start = start_[i+1];
      const int li = length_[i];
      for (CoinBigIndex j = si;j<si+li;j++) {
	double value = element_[j];
	if (fabs(value)>removeValue) {
	  index_[put]=index_[j];
	  element_[put++]=value;
	}
      }
      length_[i]=put-start_[i];
      start_[i+1] = put;
    }
    size_ = put;
  }
}

//#############################################################################

/* Really clean up matrix.
   a) eliminate all duplicate AND small elements in matrix 
   b) remove all gaps and set extraGap_ and extraMajor_ to 0.0
   c) reallocate arrays and make max lengths equal to lengths
   d) orders elements
   returns number of elements eliminated
*/
int 
CoinPackedMatrix::cleanMatrix(double threshold)
{
  if (!majorDim_) {
    extraGap_=0.0;
    extraMajor_=0.0;
    return 0;
  }
  CoinBigIndex numberEliminated =0;
  // space for eliminated
  int * mark = new int [minorDim_];
  int i;
  for (i=0;i<minorDim_;i++)
    mark[i]=-1;
  CoinBigIndex n = 0;
  for (i=0;i<majorDim_;i++) {
    CoinBigIndex k=start_[i];
    start_[i]=n;
    CoinBigIndex end = k+length_[i];
    CoinBigIndex j;
    for (j=k;j<end;j++) {
      int index = index_[j];
      if (mark[index]==-1) {
	mark[index]=j;
      } else {
	// duplicate
	int jj = mark[index];
	element_[jj] += element_[j];
	element_[j]=0.0;
      }
    }
    for (j=k;j<end;j++) {
      int index = index_[j];
      mark[index]=-1;
      if (fabs(element_[j])>=threshold) {
	element_[n]=element_[j];
	index_[n++]=index_[j];
	k++;
      }
    }
    numberEliminated += end-k;
    length_[i] = n-start_[i];
    // sort
    CoinSort_2(index_+start_[i],index_+n,element_+start_[i]);
  }
  start_[majorDim_]=n;
  size_ -= numberEliminated;
  assert (n==size_);
  delete [] mark;
  extraGap_=0.0;
  extraMajor_=0.0;
  maxMajorDim_=majorDim_;
  maxSize_=size_;
  // Now reallocate - do smallest ones first
  int * temp = CoinCopyOfArray(length_,majorDim_);
  delete [] length_;
  length_ = temp;
  CoinBigIndex * temp2 = CoinCopyOfArray(start_,majorDim_+1);
  delete [] start_;
  start_ = temp2;
  temp = CoinCopyOfArray(index_,size_);
  delete [] index_;
  index_ = temp;
  double * temp3 = CoinCopyOfArray(element_,size_);
  delete [] element_;
  element_ = temp3;
  return numberEliminated;
}

//#############################################################################

void
CoinPackedMatrix::submatrixOf(const CoinPackedMatrix& matrix,
			     const int numMajor, const int * indMajor)
{
   int i;
   int* sortedIndPtr = CoinTestIndexSet(numMajor, indMajor, matrix.majorDim_,
				       "submatrixOf");
   const int * sortedInd = sortedIndPtr == 0 ? indMajor : sortedIndPtr;

   gutsOfDestructor();

   // Count how many nonzeros there'll be
   CoinBigIndex nzcnt = 0;
   const int* length = matrix.getVectorLengths();
   for (i = 0; i < numMajor; ++i) {
      nzcnt += length[sortedInd[i]];
   }

   colOrdered_ = matrix.colOrdered_;
   maxMajorDim_ = int(numMajor * (1+extraMajor_) + 1);
   maxSize_ = (CoinBigIndex) (nzcnt * (1+extraMajor_) * (1+extraGap_) + 100);
   length_ = new int[maxMajorDim_];
   start_ = new CoinBigIndex[maxMajorDim_+1];
   start_[0]=0;
   index_ = new int[maxSize_];
   element_ = new double[maxSize_];
   majorDim_ = 0;
   minorDim_ = matrix.minorDim_;
   size_ = 0;
#ifdef CLP_NO_VECTOR
   for (i = 0; i < numMajor; ++i) {
     int j = sortedInd[i];
     CoinBigIndex start = matrix.start_[j];
     appendMajorVector(matrix.length_[j],matrix.index_+start,matrix.element_+start);
   }
#else
   for (i = 0; i < numMajor; ++i) {
      const CoinShallowPackedVector reqdBySunCC = matrix.getVector(sortedInd[i]) ;
      appendMajorVector(reqdBySunCC);
   }
#endif

   delete[] sortedIndPtr;
}

//#############################################################################

void
CoinPackedMatrix::submatrixOfWithDuplicates(const CoinPackedMatrix& matrix,
			     const int numMajor, const int * indMajor)
{
   int i;
   // we should allow duplicates - can be useful
   for (i=0; i<numMajor;i++) {
     if (indMajor[i]<0||indMajor[i]>=matrix.majorDim_)
       throw CoinError("bad index", "submatrixOfWithDuplicates", "CoinPackedMatrix");
   }
   gutsOfDestructor();

   // Count how many nonzeros there'll be
   CoinBigIndex nzcnt = 0;
   const int* length = matrix.getVectorLengths();
   for (i = 0; i < numMajor; ++i) {
      nzcnt += length[indMajor[i]];
   }

   colOrdered_ = matrix.colOrdered_;
   maxMajorDim_ = int(numMajor * (1+extraMajor_) + 1);
   maxSize_ = (CoinBigIndex) (nzcnt * (1+extraMajor_) * (1+extraGap_) + 100);
   length_ = new int[maxMajorDim_];
   start_ = new CoinBigIndex[maxMajorDim_+1];
   start_[0]=0;
   index_ = new int[maxSize_];
   element_ = new double[maxSize_];
   majorDim_ = 0;
   minorDim_ = matrix.minorDim_;
   size_ = 0;
#ifdef CLP_NO_VECTOR
   for (i = 0; i < numMajor; ++i) {
     int j = indMajor[i];
     CoinBigIndex start = matrix.start_[j];
     appendMajorVector(matrix.length_[j],matrix.index_+start,matrix.element_+start);
   }
#else
   for (i = 0; i < numMajor; ++i) {
      const CoinShallowPackedVector reqdBySunCC = matrix.getVector(indMajor[i]) ;
      appendMajorVector(reqdBySunCC);
   }
#endif

}

//#############################################################################

void
CoinPackedMatrix::copyOf(const CoinPackedMatrix& rhs)
{
   if (this != &rhs) {
      gutsOfDestructor();
      gutsOfCopyOf(rhs.colOrdered_,
		   rhs.minorDim_, rhs.majorDim_, rhs.size_,
		   rhs.element_, rhs.index_, rhs.start_, rhs.length_,
		   rhs.extraMajor_, rhs.extraGap_);
   }
}

//-----------------------------------------------------------------------------

void
CoinPackedMatrix::copyOf(const bool colordered,
			const int minor, const int major,
			const CoinBigIndex numels,
			const double * elem, const int * ind,
			const CoinBigIndex * start, const int * len,
			const double extraMajor, const double extraGap)
{
   gutsOfDestructor();
   gutsOfCopyOf(colordered, minor, major, numels, elem, ind, start, len,
		extraMajor, extraGap);
}

//#############################################################################

// This method is essentially the same as minorAppendOrthoOrdered(). However,
// since we start from an empty matrix, lots of fluff can be avoided.

void
CoinPackedMatrix::reverseOrderedCopyOf(const CoinPackedMatrix& rhs)
{
   if (this == &rhs) {
      reverseOrdering();
      return;
   }

   int i;
   colOrdered_ = !rhs.colOrdered_;
   majorDim_ = rhs.minorDim_;
   minorDim_ = rhs.majorDim_;
   size_ = rhs.size_;

   if (size_ == 0) {
     // we still need to allocate starts and lengths
     maxMajorDim_=majorDim_;
     delete[] start_;
     delete[] length_;
     delete[] index_;
     delete[] element_;
     start_ = new CoinBigIndex[maxMajorDim_ + 1];
     length_ = new int[maxMajorDim_];
     for (i = 0; i < majorDim_; ++i) {
       start_[i] = 0;
       length_[i]=0;
     }
     start_[majorDim_]=0;
     index_ = new int[maxSize_];
     element_ = new double[maxSize_];
     return;
   }

   // first compute how long each major-dimension vector will be
   // this trickery is needed because MSVC++ is not willing to delete[] a
   // 'const int *'
   int * orthoLengthPtr = rhs.countOrthoLength();
   const int * orthoLength = orthoLengthPtr;

   // Allocate sufficient space (resizeForAddingMinorVectors())
   
   const int newMaxMajorDim_ =
     CoinMax(maxMajorDim_, CoinLengthWithExtra(majorDim_, extraMajor_));

   if (newMaxMajorDim_ > maxMajorDim_) {
      maxMajorDim_ = newMaxMajorDim_;
      delete[] start_;
      delete[] length_;
      start_ = new CoinBigIndex[maxMajorDim_ + 1];
      length_ = new int[maxMajorDim_];
   }

   start_[0] = 0;
   if (extraGap_ == 0) {
      for (i = 0; i < majorDim_; ++i)
	 start_[i+1] = start_[i] + orthoLength[i];
   } else {
      const double eg = extraGap_;
      for (i = 0; i < majorDim_; ++i)
	start_[i+1] = start_[i] + CoinLengthWithExtra(orthoLength[i], eg);
   }

   const CoinBigIndex newMaxSize =
      CoinMax(maxSize_, getLastStart()+ (CoinBigIndex) extraMajor_);

   if (newMaxSize > maxSize_) {
      maxSize_ = newMaxSize;
      delete[] index_;
      delete[] element_;
      index_ = new int[maxSize_];
      element_ = new double[maxSize_];
#     ifdef ZEROFAULT
      memset(index_,0,(maxSize_*sizeof(int))) ;
      memset(element_,0,(maxSize_*sizeof(double))) ;
#     endif
   }

   // now insert the entries of matrix
   
   minorDim_ = 0;
   // Copy lengths
   CoinCopyN(orthoLength, majorDim_, length_);
   delete[] orthoLengthPtr;

   for (i = 0; i < rhs.majorDim_; ++i) {
      const CoinBigIndex last = rhs.getVectorLast(i);
      for (CoinBigIndex j = rhs.getVectorFirst(i); j != last; ++j) {
	 const int ind = rhs.index_[j];
         CoinBigIndex put = start_[ind];
         start_[ind] = put +1;
	 element_[put] = rhs.element_[j];
	 index_[put] = minorDim_;
      }
      ++minorDim_;
   }
   // and re-adjust start_
   for (i = 0; i < majorDim_; ++i) {
     start_[i] -= length_[i];
   }
}
   
//#############################################################################

void
CoinPackedMatrix::assignMatrix(const bool colordered,
			      const int minor, const int major,
			      const CoinBigIndex numels,
			      double *& elem, int *& ind,
			      CoinBigIndex *& start, int *& len,
			      const int maxmajor, const CoinBigIndex maxsize)
{
   gutsOfDestructor();
   colOrdered_ = colordered;
   element_ = elem;
   index_ = ind;
   start_ = start;
   majorDim_ = major;
   minorDim_ = minor;
   size_ = numels;
   maxMajorDim_ = maxmajor != -1 ? maxmajor : major;
   maxSize_ = maxsize != -1 ? maxsize : numels;
   if (len == NULL) {
     delete [] length_;
     length_ = new int[maxMajorDim_];
     std::adjacent_difference(start + 1, start + (major + 1), length_);
   } else {
      length_ = len;
   }
   elem = NULL;
   ind = NULL;
   start = NULL;
   len = NULL;
}

//#############################################################################

CoinPackedMatrix &
CoinPackedMatrix::operator=(const CoinPackedMatrix& rhs)
{
   if (this != &rhs) {
      gutsOfDestructor();
      extraGap_=rhs.extraGap_;
      extraMajor_=rhs.extraMajor_;
      gutsOfOpEqual(rhs.colOrdered_,
		    rhs.minorDim_,  rhs.majorDim_, rhs.size_,
		    rhs.element_, rhs.index_, rhs.start_, rhs.length_);
   }
   return *this;
}

//#############################################################################

void
CoinPackedMatrix::reverseOrdering()
{
   CoinPackedMatrix m;
   m.extraGap_ = extraMajor_;
   m.extraMajor_ = extraGap_;
   m.reverseOrderedCopyOf(*this);
   swap(m);
}

//-----------------------------------------------------------------------------

void
CoinPackedMatrix::transpose()
{
   colOrdered_ = ! colOrdered_;
}

//-----------------------------------------------------------------------------

void
CoinPackedMatrix::swap(CoinPackedMatrix& m)
{
   std::swap(colOrdered_,  m.colOrdered_);
   std::swap(extraGap_,	   m.extraGap_);
   std::swap(extraMajor_,  m.extraMajor_);
   std::swap(element_, 	   m.element_);
   std::swap(index_,	   m.index_);
   std::swap(start_,	   m.start_);
   std::swap(length_,	   m.length_);
   std::swap(majorDim_,	   m.majorDim_);
   std::swap(minorDim_,	   m.minorDim_);
   std::swap(size_,	   m.size_);
   std::swap(maxMajorDim_, m.maxMajorDim_);
   std::swap(maxSize_,     m.maxSize_);
}

//#############################################################################
//#############################################################################

void
CoinPackedMatrix::times(const double * x, double * y) const 
{
   if (colOrdered_)
      timesMajor(x, y);
   else
      timesMinor(x, y);
}

//-----------------------------------------------------------------------------
#ifndef CLP_NO_VECTOR
void
CoinPackedMatrix::times(const CoinPackedVectorBase& x, double * y) const 
{
   if (colOrdered_)
      timesMajor(x, y);
   else
      timesMinor(x, y);
}
#endif
//-----------------------------------------------------------------------------

void
CoinPackedMatrix::transposeTimes(const double * x, double * y) const 
{
   if (colOrdered_)
      timesMinor(x, y);
   else
      timesMajor(x, y);
}

//-----------------------------------------------------------------------------
#ifndef CLP_NO_VECTOR
void
CoinPackedMatrix::transposeTimes(const CoinPackedVectorBase& x, double * y) const
{
   if (colOrdered_)
      timesMinor(x, y);
   else
      timesMajor(x, y);
}
#endif
//#############################################################################
//#############################################################################

int *
CoinPackedMatrix::countOrthoLength() const
{
   int * orthoLength = new int[minorDim_];
   CoinZeroN(orthoLength, minorDim_);
   if (size_!=start_[majorDim_]) {
     // has gaps
     for (int i = 0; i <majorDim_ ; ++i) {
       const CoinBigIndex first = start_[i];
       const CoinBigIndex last = first + length_[i];
       for (CoinBigIndex j = first; j < last; ++j) {
         assert( index_[j] < minorDim_ && index_[j]>=0);
         ++orthoLength[index_[j]];
       }
     }
   } else {
     // no gaps 
     const CoinBigIndex last = start_[majorDim_];
     for (CoinBigIndex j = 0; j < last; ++j) {
       assert( index_[j] < minorDim_ && index_[j]>=0);
       ++orthoLength[index_[j]];
     }
   }
   return orthoLength;
}

//#############################################################################
/* Returns an array containing major indices.  The array is
    getNumElements long and if getVectorStarts() is 0,2,5 then
    the array would start 0,0,1,1,1,2...
    This method is provided to go back from a packed format
    to a triple format.
    The returned array is allocated with <code>new int[]</code>,
    free it with  <code>delete[]</code>. */
int * 
CoinPackedMatrix::getMajorIndices() const
{
  // Check valid
  if (!majorDim_||start_[majorDim_]!=size_)
    return NULL;
  int * array = new int [size_];
  for (int i=0;i<majorDim_;i++) {
    for (CoinBigIndex k=start_[i];k<start_[i+1];k++)
      array[k]=i;
  }
  return array;
}
//#############################################################################

void
CoinPackedMatrix::appendMajorVector(const int vecsize,
				   const int *vecind,
				   const double *vecelem)
{
#ifdef COIN_DEBUG
  for (int i = 0; i < vecsize; ++i) {
    if (vecind[i] < 0 )
      throw CoinError("out of range index",
		     "appendMajorVector", "CoinPackedMatrix");
  }
#if 0   
  if (std::find_if(vecind, vecind + vecsize,
		   compose2(logical_or<bool>(),
			    bind2nd(less<int>(), 0),
			    bind2nd(greater_equal<int>(), minorDim_))) !=
      vecind + vecsize)
    throw CoinError("out of range index",
		   "appendMajorVector", "CoinPackedMatrix");
#endif
#endif
  
  if (majorDim_ == maxMajorDim_ || vecsize > maxSize_ - getLastStart()) {
    resizeForAddingMajorVectors(1, &vecsize);
  }

  // got to get this again since it might change!
  const CoinBigIndex last = getLastStart();

  // OK, now just append the major-dimension vector to the end

  length_[majorDim_] = vecsize;
  CoinMemcpyN(vecind, vecsize, index_ + last);
  CoinMemcpyN(vecelem, vecsize, element_ + last);
  if (majorDim_ == 0)
    start_[0] = 0;
  start_[majorDim_ + 1] =
    CoinMin(last + CoinLengthWithExtra(vecsize, extraGap_), maxSize_ );

   // LL: Do we want to allow appending a vector that has more entries than
   // the current size?
   if (vecsize > 0) {
     minorDim_ = CoinMax(minorDim_,
			(*std::max_element(vecind, vecind+vecsize)) + 1);
   }

   ++majorDim_;
   size_ += vecsize;
}

//-----------------------------------------------------------------------------
#ifndef CLP_NO_VECTOR
void
CoinPackedMatrix::appendMajorVector(const CoinPackedVectorBase& vec)
{
   appendMajorVector(vec.getNumElements(),
		     vec.getIndices(), vec.getElements());
}
//-----------------------------------------------------------------------------

void
CoinPackedMatrix::appendMajorVectors(const int numvecs,
				    const CoinPackedVectorBase * const * vecs)
{
  int i;
  CoinBigIndex nz = 0;
  for (i = 0; i < numvecs; ++i)
    nz += CoinLengthWithExtra(vecs[i]->getNumElements(), extraGap_);
  reserve(majorDim_ + numvecs, getLastStart() + nz);
  for (i = 0; i < numvecs; ++i)
    appendMajorVector(*vecs[i]);
}
#endif

//#############################################################################

void
CoinPackedMatrix::appendMinorVector(const int vecsize,
				   const int *vecind,
				   const double *vecelem)
{
  int i;
#ifdef COIN_DEBUG
  for (i = 0; i < vecsize; ++i) {
    if (vecind[i] < 0 || vecind[i] >= majorDim_)
      throw CoinError("out of range index",
		     "appendMinorVector", "CoinPackedMatrix");
  }
#if 0   
  if (std::find_if(vecind, vecind + vecsize,
		   compose2(logical_or<bool>(),
			    bind2nd(less<int>(), 0),
			    bind2nd(greater_equal<int>(), majorDim_))) !=
      vecind + vecsize)
    throw CoinError("out of range index",
		   "appendMinorVector", "CoinPackedMatrix");
#endif
#endif

  // test that there's a gap at the end of every major-dimension vector where
  // we want to add a new entry
   
  for (i = vecsize - 1; i >= 0; --i) {
    const int j = vecind[i];
    if (start_[j] + length_[j] == start_[j+1])
      break;
  }

  if (i >= 0) {
    int * addedEntries = new int[majorDim_];
    memset(addedEntries, 0, majorDim_ * sizeof(int));
    for (i = vecsize - 1; i >= 0; --i)
      addedEntries[vecind[i]] = 1;
    resizeForAddingMinorVectors(addedEntries);
    delete[] addedEntries;
  }

  // OK, now insert the entries of the minor-dimension vector
  for (i = vecsize - 1; i >= 0; --i) {
    const int j = vecind[i];
    const CoinBigIndex posj = start_[j] + (length_[j]++);
    index_[posj] = minorDim_;
    element_[posj] = vecelem[i];
  }

  ++minorDim_;
  size_ += vecsize;
}

//-----------------------------------------------------------------------------
#ifndef CLP_NO_VECTOR
void
CoinPackedMatrix::appendMinorVector(const CoinPackedVectorBase& vec)
{
   appendMinorVector(vec.getNumElements(),
		     vec.getIndices(), vec.getElements());
}

//-----------------------------------------------------------------------------

void
CoinPackedMatrix::appendMinorVectors(const int numvecs,
				    const CoinPackedVectorBase * const * vecs)
{
  if (numvecs == 0)
    return;

  int i;

  int * addedEntries = new int[majorDim_];
  CoinZeroN(addedEntries, majorDim_);
  for (i = numvecs - 1; i >= 0; --i) {
    const int vecsize = vecs[i]->getNumElements();
    const int* vecind = vecs[i]->getIndices();
    for (int j = vecsize - 1; j >= 0; --j) {
#ifdef COIN_DEBUG
      if (vecind[j] < 0 || vecind[j] >= majorDim_)
	throw CoinError("out of range index", "appendMinorVectors",
		       "CoinPackedMatrix");
#endif
      ++addedEntries[vecind[j]];
    }
  }
 
  for (i = majorDim_ - 1; i >= 0; --i) {
    if (start_[i] + length_[i] + addedEntries[i] > start_[i+1])
      break;
  }
  if (i >= 0)
    resizeForAddingMinorVectors(addedEntries);
  delete[] addedEntries;

  // now insert the entries of the vectors
  for (i = 0; i < numvecs; ++i) {
    const int vecsize = vecs[i]->getNumElements();
    const int* vecind = vecs[i]->getIndices();
    const double* vecelem = vecs[i]->getElements();
    for (int j = vecsize - 1; j >= 0; --j) {
      const int ind = vecind[j];
      element_[start_[ind] + length_[ind]] = vecelem[j];
      index_[start_[ind] + (length_[ind]++)] = minorDim_;
    }
    ++minorDim_;
    size_ += vecsize;
  }
}
#endif

//#############################################################################
//#############################################################################

void
CoinPackedMatrix::majorAppendSameOrdered(const CoinPackedMatrix& matrix)
{
   if (minorDim_ != matrix.minorDim_) {
      throw CoinError("dimension mismatch", "rightAppendSameOrdered",
		     "CoinPackedMatrix");
   }
   if (matrix.majorDim_ == 0)
      return;

   int i;
   if (majorDim_ + matrix.majorDim_ > maxMajorDim_ ||
       getLastStart() + matrix.getLastStart() > maxSize_) {
      // we got to resize before we add. note that the resizing method
      // properly fills out start_ and length_ for the major-dimension
      // vectors to be added!
      resizeForAddingMajorVectors(matrix.majorDim_, matrix.length_);
      start_ += majorDim_;
      for (i = 0; i < matrix.majorDim_; ++i) {
	 const int l = matrix.length_[i];
	 CoinMemcpyN(matrix.index_ + matrix.start_[i], l,
			   index_ + start_[i]);
	 CoinMemcpyN(matrix.element_ + matrix.start_[i], l,
			   element_ + start_[i]);
      }
      start_ -= majorDim_;
   } else {
      start_ += majorDim_;
      length_ += majorDim_;
      for (i = 0; i < matrix.majorDim_; ++i) {
	 const int l = matrix.length_[i];
	 CoinMemcpyN(matrix.index_ + matrix.start_[i], l,
			   index_ + start_[i]);
	 CoinMemcpyN(matrix.element_ + matrix.start_[i], l,
			   element_ + start_[i]);
	 start_[i+1] = start_[i] + matrix.start_[i+1] - matrix.start_[i];
	 length_[i] = l;
      }
      start_ -= majorDim_;
      length_ -= majorDim_;
   }
   majorDim_ += matrix.majorDim_;
   size_ += matrix.size_;
}

//-----------------------------------------------------------------------------
   
void
CoinPackedMatrix::minorAppendSameOrdered(const CoinPackedMatrix& matrix)
{
   if (majorDim_ != matrix.majorDim_) {
      throw CoinError("dimension mismatch", "bottomAppendSameOrdered",
		      "CoinPackedMatrix");
   }
   if (matrix.minorDim_ == 0)
      return;

   int i;
   for (i = majorDim_ - 1; i >= 0; --i) {
      if (start_[i] + length_[i] + matrix.length_[i] > start_[i+1])
	 break;
   }
   if (i >= 0)
      resizeForAddingMinorVectors(matrix.length_);

   // now insert the entries of matrix
   for (i = majorDim_ - 1; i >= 0; --i) {
      const int l = matrix.length_[i];
      std::transform(matrix.index_ + matrix.start_[i],
		matrix.index_ + (matrix.start_[i] + l),
		index_ + (start_[i] + length_[i]),
		std::bind2nd(std::plus<int>(), minorDim_));
      CoinMemcpyN(matrix.element_ + matrix.start_[i], l,
		       element_ + (start_[i] + length_[i]));
      length_[i] += l;
   }
   minorDim_ += matrix.minorDim_;
   size_ += matrix.size_;
}
   
//-----------------------------------------------------------------------------
   
void
CoinPackedMatrix::majorAppendOrthoOrdered(const CoinPackedMatrix& matrix)
{
   if (minorDim_ != matrix.majorDim_) {
      throw CoinError("dimension mismatch", "majorAppendOrthoOrdered",
		     "CoinPackedMatrix");
      }
   if (matrix.majorDim_ == 0)
      return;

   int i;
   CoinBigIndex j;
   // this trickery is needed because MSVC++ is not willing to delete[] a
   // 'const int *'
   int * orthoLengthPtr = matrix.countOrthoLength();
   const int * orthoLength = orthoLengthPtr;

   if (majorDim_ + matrix.minorDim_ > maxMajorDim_) {
      resizeForAddingMajorVectors(matrix.minorDim_, orthoLength);
   } else {
     const double extra_gap = extraGap_;
     start_ += majorDim_;
     for (i = 0; i < matrix.minorDim_ ; ++i) {
       start_[i+1] = start_[i] + CoinLengthWithExtra(orthoLength[i], extra_gap);
     }
     start_ -= majorDim_;
     if (start_[majorDim_ + matrix.minorDim_] > maxSize_) {
       resizeForAddingMajorVectors(matrix.minorDim_, orthoLength);
     }
   }
   // At this point everything is big enough to accommodate the new entries.
   // Also, start_ is set to the correct starting points for all the new
   // major-dimension vectors. The length of the new major-dimension vectors
   // may or may not be correctly set. Hence we just zero them out and they'll
   // be set when the entries are actually added below.

   start_ += majorDim_;
   length_ += majorDim_;

   CoinZeroN(length_, matrix.minorDim_);

   for (i = 0; i < matrix.majorDim_; ++i) {
      const CoinBigIndex last = matrix.getVectorLast(i);
      for (j = matrix.getVectorFirst(i); j < last; ++j) {
	 const int ind = matrix.index_[j];
	 element_[start_[ind] + length_[ind]] = matrix.element_[j];
	 index_[start_[ind] + (length_[ind]++)] = i;
      }
   }
   
   length_ -= majorDim_;
   start_ -= majorDim_;

   // We need to update majorDim_ and size_.  We can just add in from matrix
   majorDim_ += matrix.minorDim_;
   size_ += matrix.size_;

   delete[] orthoLengthPtr;
}

//-----------------------------------------------------------------------------
   
void
CoinPackedMatrix::minorAppendOrthoOrdered(const CoinPackedMatrix& matrix)
{
   if (majorDim_ != matrix.minorDim_) {
      throw CoinError("dimension mismatch", "bottomAppendOrthoOrdered",
		     "CoinPackedMatrix");
      }
   if (matrix.majorDim_ == 0)
      return;

   int i;
   // first compute how many entries will be added to each major-dimension
   // vector, and if needed, resize the matrix to accommodate all
   // this trickery is needed because MSVC++ is not willing to delete[] a
   // 'const int *'
   int * addedEntriesPtr = matrix.countOrthoLength();
   const int * addedEntries = addedEntriesPtr;
   for (i = majorDim_ - 1; i >= 0; --i) {
      if (start_[i] + length_[i] + addedEntries[i] > start_[i+1])
	 break;
   }
   if (i >= 0)
      resizeForAddingMinorVectors(addedEntries);
   delete[] addedEntriesPtr;

   // now insert the entries of matrix
   for (i = 0; i < matrix.majorDim_; ++i) {
      const CoinBigIndex last = matrix.getVectorLast(i);
      for (CoinBigIndex j = matrix.getVectorFirst(i); j != last; ++j) {
	 const int ind = matrix.index_[j];
	 element_[start_[ind] + length_[ind]] = matrix.element_[j];
	 index_[start_[ind] + (length_[ind]++)] = minorDim_;
      }
      ++minorDim_;
   }
   size_ += matrix.size_;
}

//#############################################################################
//#############################################################################

void
CoinPackedMatrix::deleteMajorVectors(const int numDel,
				    const int * indDel)
{
   int *sortedDelPtr = CoinTestIndexSet(numDel, indDel, majorDim_,
				       "deleteMajorVectors");
   const int * sortedDel = sortedDelPtr == 0 ? indDel : sortedDelPtr;

   if (numDel == majorDim_) {
      // everything is deleted
      majorDim_ = 0;
      minorDim_ = 0;
      size_ = 0;
      if (sortedDelPtr)
	 delete[] sortedDelPtr;
      // Get rid of memory as well
      maxMajorDim_ = 0;
      delete [] length_;
      length_ = NULL;
      delete [] start_;
      start_ = new CoinBigIndex[1];
      start_[0]=0;;
      delete [] element_;
      element_=NULL;
      delete [] index_;
      index_=NULL;
      maxSize_ = 0;
      return;
   }

   CoinBigIndex deleted = 0;
   const int last = numDel - 1;
   for (int i = 0; i < last; ++i) {
      const int ind = sortedDel[i];
      const int ind1 = sortedDel[i+1];
      deleted += length_[ind];
      if (ind1 - ind > 1) {
	 CoinCopy(start_ + (ind + 1), start_ + ind1, start_ + (ind - i));
	 CoinCopy(length_ + (ind + 1), length_ + ind1, length_ + (ind - i));
      }
   }

   // copy the last block of length_ and start_
   const int ind = sortedDel[last];
   deleted += length_[ind];
   if (sortedDel[last] != majorDim_ - 1) {
      const int ind1 = majorDim_;
      CoinCopy(start_ + (ind + 1), start_ + ind1, start_ + (ind - last));
      CoinCopy(length_ + (ind + 1), length_ + ind1, length_ + (ind - last));
   }
   majorDim_ -= numDel;
   const int lastlength = CoinLengthWithExtra(length_[majorDim_-1], extraGap_);
   start_[majorDim_] = CoinMin(start_[majorDim_-1] + lastlength, maxSize_);
   size_ -= deleted;

   // if the very first major vector was deleted then copy the new first major
   // vector to the beginning to make certain that start_[0] is 0. This may
   // not be necessary, but better safe than sorry...
   if (sortedDel[0] == 0) {
      CoinCopyN(index_ + start_[0], length_[0], index_);
      CoinCopyN(element_ + start_[0], length_[0], element_);
      start_[0] = 0;
   }

   delete[] sortedDelPtr;
}

//#############################################################################

void
CoinPackedMatrix::deleteMinorVectors(const int numDel,
				    const int * indDel)
{
   if (numDel == minorDim_) {
     // everything is deleted
     minorDim_ = 0;
     size_ = 0;
     // Get rid of as much memory as possible
     memset(length_,0,majorDim_*sizeof(int));
     memset(start_,0,(majorDim_+1)*sizeof(CoinBigIndex ));
     delete [] element_;
     element_=NULL;
     delete [] index_;
     index_=NULL;
     maxSize_ = 0;
     return;
   }
  int i, j, k;

  // first compute the new index of every row
  int* newindexPtr = new int[minorDim_];
  CoinIotaN(newindexPtr, minorDim_, 0);
  for (j = 0; j < numDel; ++j) {
    const int ind = indDel[j];
#ifdef COIN_DEBUG
    if (ind < 0 || ind >= minorDim_)
      throw CoinError("out of range index",
		     "deleteMinorVectors", "CoinPackedMatrix");
    if (newindexPtr[ind] == -1)
      throw CoinError("duplicate index",
		     "deleteMinorVectors", "CoinPackedMatrix");
#endif
    newindexPtr[ind] = -1;
  }
  for (i = 0, k = 0; i < minorDim_; ++i) {
    if (newindexPtr[i] != -1) {
      newindexPtr[i] = k++;
    }
  }
  // Now crawl through the matrix
  const int * newindex = newindexPtr;
#ifdef TAKEOUT
  int mcount[400];
  memset(mcount,0,400*sizeof(int));
  for (i = 0; i < majorDim_; ++i) {
    int * index = index_ + start_[i];
    double * elem = element_ + start_[i];
    const int length_i = length_[i];
    for (j = 0, k = 0; j < length_i; ++j) {
      mcount[index[j]]++;
    }
  }
  for (i=0;i<minorDim_;i++) {
    if (mcount[i]==10||mcount[i]==15) {
      if (newindex[i]>=0)
	printf("Keeping original row %d (new %d) with count of %d\n",
	       i,newindex[i],mcount[i]);
      else
	printf("deleting row %d with count of %d\n",
	       i,mcount[i]);
    }
  }
#endif
  int deleted = 0;
  for (i = 0; i < majorDim_; ++i) {
    int * index = index_ + start_[i];
    double * elem = element_ + start_[i];
    const int length_i = length_[i];
    for (j = 0, k = 0; j < length_i; ++j) {
      const int ind = newindex[index[j]];
      if (ind != -1) {
	index[k] = ind;
	elem[k++] = elem[j];
      }
    }
    deleted += length_i - k;
    length_[i] = k;
  }

  delete[] newindexPtr;

  minorDim_ -= numDel;
  size_ -= deleted;
}

//#############################################################################
//#############################################################################

void
CoinPackedMatrix::timesMajor(const double * x, double * y) const 
{
   memset(y, 0, minorDim_ * sizeof(double));
   for (int i = majorDim_ - 1; i >= 0; --i) {
      const double x_i = x[i];
      if (x_i != 0.0) {
	 const CoinBigIndex last = getVectorLast(i);
	 for (CoinBigIndex j = getVectorFirst(i); j < last; ++j)
	    y[index_[j]] += x_i * element_[j];
      }
   }
}

//-----------------------------------------------------------------------------
#ifndef CLP_NO_VECTOR
void
CoinPackedMatrix::timesMajor(const CoinPackedVectorBase& x, double * y) const 
{
   memset(y, 0, minorDim_ * sizeof(double));
   for (CoinBigIndex i = x.getNumElements() - 1; i >= 0; --i) {
      const double x_i = x.getElements()[i];
      if (x_i != 0.0) {
	 const int ind = x.getIndices()[i];
	 const CoinBigIndex last = getVectorLast(ind);
	 for (CoinBigIndex j = getVectorFirst(ind); j < last; ++j)
	    y[index_[j]] += x_i * element_[j];
      }
   }
}
#endif
//-----------------------------------------------------------------------------

void
CoinPackedMatrix::timesMinor(const double * x, double * y) const 
{
   memset(y, 0, majorDim_ * sizeof(double));
   for (int i = majorDim_ - 1; i >= 0; --i) {
      double y_i = 0;
      const CoinBigIndex last = getVectorLast(i);
      for (CoinBigIndex j = getVectorFirst(i); j < last; ++j)
	 y_i += x[index_[j]] * element_[j];
      y[i] = y_i;
   }
}

//-----------------------------------------------------------------------------
#ifndef CLP_NO_VECTOR
void
CoinPackedMatrix::timesMinor(const CoinPackedVectorBase& x, double * y) const 
{
   memset(y, 0, majorDim_ * sizeof(double));
   for (int i = majorDim_ - 1; i >= 0; --i) {
      double y_i = 0;
      const CoinBigIndex last = getVectorLast(i);
      for (CoinBigIndex j = getVectorFirst(i); j < last; ++j)
	 y_i += x[index_[j]] * element_[j];
      y[i] = y_i;
   }
}
#endif
//#############################################################################
//#############################################################################

CoinPackedMatrix::CoinPackedMatrix() :
   colOrdered_(true),
   extraGap_(0.25),
   extraMajor_(0.25),
   element_(0), 
   index_(0),
   length_(0),
   majorDim_(0),
   minorDim_(0),
   size_(0),
   maxMajorDim_(0),
   maxSize_(0) 
{
  start_ = new CoinBigIndex[1];
  start_[0] = 0;
}

//-----------------------------------------------------------------------------

CoinPackedMatrix::CoinPackedMatrix(const bool colordered,
				 const double extraMajor,
				 const double extraGap) :
   colOrdered_(colordered),
   extraGap_(extraGap),
   extraMajor_(extraMajor),
   element_(0), 
   index_(0),
   length_(0),
   majorDim_(0),
   minorDim_(0),
   size_(0),
   maxMajorDim_(0),
   maxSize_(0)
{
  start_ = new CoinBigIndex[1];
  start_[0] = 0;
}

//-----------------------------------------------------------------------------

CoinPackedMatrix::CoinPackedMatrix(const bool colordered,
				 const int minor, const int major,
				 const CoinBigIndex numels,
				 const double * elem, const int * ind,
				 const CoinBigIndex * start, const int * len,
				 const double extraMajor,
				 const double extraGap) :
   colOrdered_(colordered),
   extraGap_(extraGap),
   extraMajor_(extraMajor),
   element_(NULL),
   index_(NULL),
   start_(NULL),
   length_(NULL),
   majorDim_(0),
   minorDim_(0),
   size_(0),
   maxMajorDim_(0),
   maxSize_(0)
{
   gutsOfOpEqual(colordered, minor, major, numels, elem, ind, start, len);
}

//-----------------------------------------------------------------------------
   
CoinPackedMatrix::CoinPackedMatrix(const bool colordered,
				 const int minor, const int major,
         const CoinBigIndex numels,
         const double * elem, const int * ind,
         const CoinBigIndex * start, const int * len) :
   colOrdered_(colordered),
   extraGap_(0.0),
   extraMajor_(0.0),
   element_(NULL),
   index_(NULL),
   start_(NULL),
   length_(NULL),
   majorDim_(0),
   minorDim_(0),
   size_(0),
   maxMajorDim_(0),
   maxSize_(0)
{
     gutsOfOpEqual(colordered, minor, major, numels, elem, ind, start, len);
}
  
//-----------------------------------------------------------------------------
// makes column ordered from triplets and takes out duplicates 
// will be sorted 
//
// This is an interesting in-place sorting algorithm; 
// We have triples, and want to sort them so that triples with the same column
// are adjacent.
// We begin by computing how many entries there are for each column (columnCount)
// and using that to compute where each set of column entries will *end* (startColumn).
// As we drop entries into place, startColumn is decremented until it contains
// the position where the column entries *start*.
// The invalid column index -2 means there's a "hole" in that position;
// the invalid column index -1 means the entry in that spot is "where it wants to go".
// Initially, no one is where they want to go.
// Going back to front,
//    if that entry is where it wants to go
//    then leave it there
//    otherwise pick it up (which leaves a hole), and 
//	      for as long as you have an entry in your right hand,
//	- pick up the entry (with your left hand) in the position where the one in 
//		your right hand wants to go;
//	- pass the entry in your left hand to your right hand;
//	- was that entry really just the "hole"?  If so, stop.
// It could be that all the entries get shuffled in the first loop iteration
// and all the rest just confirm that everyone is happy where they are.
// We never move an entry that is where it wants to go, so entries are moved at
// most once.  They may not change position if they happen to initially be
// where they want to go when the for loop gets to them.
// It depends on how many subpermutations the triples initially defined.
// Each while loop takes care of one permutation.
// The while loop has to stop, because each time around we mark one entry as happy.
// We can't run into a happy entry, because we are decrementing the startColumn
// all the time, so we must be running into new entries.
// Once we've processed all the slots for a column, it cannot be the case that
// there are any others that want to go there.
// This all means that we eventually must run into the hole.
CoinPackedMatrix::CoinPackedMatrix(
     const bool colordered,
     const int * indexRow ,
     const int * indexColumn,
     const double * element, 
     CoinBigIndex numberElements ) 
     :
   colOrdered_(colordered),
     extraGap_(0.0),
     extraMajor_(0.0),
     element_(NULL),
     index_(NULL),
     start_(NULL),
     length_(NULL),
     majorDim_(0),
     minorDim_(0),
     size_(0),
     maxMajorDim_(0),
     maxSize_(0)
{
     CoinAbsFltEq eq;
       int * colIndices = new int[numberElements];
       int * rowIndices = new int[numberElements];
       double * elements = new double[numberElements];
       CoinCopyN(element,numberElements,elements);
     if ( colordered ) {
       CoinCopyN(indexColumn,numberElements,colIndices);
       CoinCopyN(indexRow,numberElements,rowIndices);
     }
     else {
       CoinCopyN(indexColumn,numberElements,rowIndices);
       CoinCopyN(indexRow,numberElements,colIndices);
     }

  int numberRows;
  int numberColumns;
  if (numberElements ) {
    numberRows = *std::max_element(rowIndices,rowIndices+numberElements)+1;
    numberColumns = *std::max_element(colIndices,colIndices+numberElements)+1;
  } else {
    numberRows = 0;
    numberColumns = 0;
  }
  int * rowCount = new int[numberRows];
  int * columnCount = new int[numberColumns];
  CoinBigIndex * startColumn = new CoinBigIndex[numberColumns+1];
  int * lengths = new int[numberColumns+1];

  int iColumn,i;
  CoinBigIndex k;
  for (i=0;i<numberRows;i++) {
    rowCount[i]=0;
  }
  for (i=0;i<numberColumns;i++) {
    columnCount[i]=0;
  }
  for (i=0;i<numberElements;i++) {
    int iRow=rowIndices[i];
    int iColumn=colIndices[i];
    rowCount[iRow]++;
    columnCount[iColumn]++;
  }
  CoinBigIndex iCount=0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    /* position after end of Column */
    iCount+=columnCount[iColumn];
    startColumn[iColumn]=iCount;
  } /* endfor */
  startColumn[iColumn]=iCount;
  for (k=numberElements-1;k>=0;k--) {
    iColumn=colIndices[k];
    if (iColumn>=0) {
      /* pick up the entry with your right hand */
      double value = elements[k];
      int iRow=rowIndices[k];
      colIndices[k]=-2;	/* the hole */

      while (1) {
	/* pick this up with your left */
        CoinBigIndex iLook=startColumn[iColumn]-1;
        double valueSave=elements[iLook];
        int iColumnSave=colIndices[iLook];
        int iRowSave=rowIndices[iLook];

	/* put the right-hand entry where it wanted to go */
        startColumn[iColumn]=iLook;
        elements[iLook]=value;
        rowIndices[iLook]=iRow;
        colIndices[iLook]=-1;	/* mark it as being where it wants to be */

	/* there was something there */
        if (iColumnSave>=0) {
          iColumn=iColumnSave;
          value=valueSave;
          iRow=iRowSave;
	} else if (iColumnSave == -2) {	/* that was the hole */
          break;
	} else {
	  assert(1==0);	/* should never happen */
	}
	/* endif */
      } /* endwhile */
    } /* endif */
  } /* endfor */

  /* now pack the elements and combine entries with the same row and column */
  /* also, drop entries with "small" coefficients */
  numberElements=0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    CoinBigIndex start=startColumn[iColumn];
    CoinBigIndex end =startColumn[iColumn+1];
    lengths[iColumn]=0;
    startColumn[iColumn]=numberElements;
    if (end>start) {
      int lastRow;
      double lastValue;
      // sorts on indices dragging elements with
      CoinSort_2(rowIndices+start,rowIndices+end,elements+start,CoinFirstLess_2<int, double>());
      lastRow=rowIndices[start];
      lastValue=elements[start];
      for (i=start+1;i<end;i++) {
        int iRow=rowIndices[i];
        double value=elements[i];
        if (iRow>lastRow) {
          //if(fabs(lastValue)>tolerance) {
          if(!eq(lastValue,0.0)) {
            rowIndices[numberElements]=lastRow;
            elements[numberElements]=lastValue;
            numberElements++;
            lengths[iColumn]++;
          }
          lastRow=iRow;
          lastValue=value;
        } else {
          lastValue+=value;
        } /* endif */
      } /* endfor */
      //if(fabs(lastValue)>tolerance) {
      if(!eq(lastValue,0.0)) {
        rowIndices[numberElements]=lastRow;
        elements[numberElements]=lastValue;
        numberElements++;
        lengths[iColumn]++;
      }
    }
  } /* endfor */
  startColumn[numberColumns]=numberElements;
#if 0
  gutsOfOpEqual(colordered,numberRows,numberColumns,numberElements,elements,rowIndices,startColumn,lengths);
  
  delete [] rowCount;
  delete [] columnCount;
  delete [] startColumn;
  delete [] lengths;

  delete [] colIndices;
  delete [] rowIndices;
  delete [] elements;
#else
  assignMatrix(colordered,numberRows,numberColumns,numberElements,
    elements,rowIndices,startColumn,lengths); 
  delete [] rowCount;
  delete [] columnCount;
  delete [] lengths;
  delete [] colIndices;
#endif

}

//-----------------------------------------------------------------------------

CoinPackedMatrix::CoinPackedMatrix (const CoinPackedMatrix & rhs) :
   colOrdered_(true),
   extraGap_(0.0),
   extraMajor_(0.0),
   element_(0), 
   index_(0),
   start_(0),
   length_(0),
   majorDim_(0),
   minorDim_(0),
   size_(0),
   maxMajorDim_(0),
   maxSize_(0)
{
  bool hasGaps = rhs.size_<rhs.start_[rhs.majorDim_];
  if (!hasGaps&&!rhs.extraMajor_) {
   gutsOfCopyOfNoGaps(rhs.colOrdered_,
		rhs.minorDim_, rhs.majorDim_,
                      rhs.element_, rhs.index_, rhs.start_);
  } else {
   gutsOfCopyOf(rhs.colOrdered_,
		rhs.minorDim_, rhs.majorDim_, rhs.size_,
		rhs.element_, rhs.index_, rhs.start_, rhs.length_,
		rhs.extraMajor_, rhs.extraGap_);
  }
}
// Subset constructor (without gaps)
CoinPackedMatrix::CoinPackedMatrix (const CoinPackedMatrix & rhs,
				    int numberRows, const int * whichRow,
				    int numberColumns, 
				    const int * whichColumn) :
   colOrdered_(true),
   extraGap_(0.0),
   extraMajor_(0.0),
   element_(NULL), 
   index_(NULL),
   start_(NULL),
   length_(NULL),
   majorDim_(0),
   minorDim_(0),
   size_(0),
   maxMajorDim_(0),
   maxSize_(0)
{
  if (numberRows<=0||numberColumns<=0) {
    start_ = new CoinBigIndex[1];
    start_[0] = 0;
  } else {
    if (!rhs.colOrdered_) {
      // just swap lists
      colOrdered_=false;
      const int * temp = whichRow;
      whichRow = whichColumn;
      whichColumn = temp;
      int n = numberRows;
      numberRows = numberColumns;
      numberColumns = n;
    }
    const double * element1 = rhs.element_;
    const int * index1 = rhs.index_;
    const CoinBigIndex * start1 = rhs.start_;
    const int * length1 = rhs.length_;

    majorDim_ = numberColumns;
    maxMajorDim_ = numberColumns;
    minorDim_ = numberRows;
    // Throw exception if rhs empty
    if (rhs.majorDim_ <= 0 || rhs.minorDim_ <= 0)
      throw CoinError("empty rhs", "subset constructor", "CoinPackedMatrix");
    // Array to say if an old row is in new copy
    int * newRow = new int [rhs.minorDim_];
    int iRow;
    for (iRow=0;iRow<rhs.minorDim_;iRow++) 
      newRow[iRow] = -1;
    // and array for duplicating rows
    int * duplicateRow = new int [numberRows];
    int numberBad=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      duplicateRow[iRow] = -1;
      int kRow = whichRow[iRow];
      if (kRow>=0  && kRow < rhs.minorDim_) {
	if (newRow[kRow]<0) {
	  // first time
	  newRow[kRow]=iRow;
	} else {
	  // duplicate
	  int lastRow = newRow[kRow];
	  newRow[kRow]=iRow;
	  duplicateRow[iRow] = lastRow;
	}
      } else {
	// bad row
	numberBad++;
      }
    }

    if (numberBad)
      throw CoinError("bad minor entries", 
		      "subset constructor", "CoinPackedMatrix");
    // now get size and check columns
    size_ = 0;
    int iColumn;
    numberBad=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      int kColumn = whichColumn[iColumn];
      if (kColumn>=0  && kColumn <rhs.majorDim_) {
	CoinBigIndex i;
	for (i=start1[kColumn];i<start1[kColumn]+length1[kColumn];i++) {
	  int kRow = index1[i];
	  kRow = newRow[kRow];
	  while (kRow>=0) {
	    size_++;
	    kRow = duplicateRow[kRow];
	  }
	}
      } else {
	// bad column
	numberBad++;
      }
    }
    if (numberBad)
      throw CoinError("bad major entries", 
		      "subset constructor", "CoinPackedMatrix");
    // now create arrays
    maxSize_=CoinMax((CoinBigIndex) 1,size_);
    start_ = new CoinBigIndex [numberColumns+1];
    length_ = new int [numberColumns];
    index_ = new int[maxSize_];
    element_ = new double [maxSize_];
    // and fill them
    size_ = 0;
    start_[0]=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      int kColumn = whichColumn[iColumn];
      CoinBigIndex i;
      for (i=start1[kColumn];i<start1[kColumn]+length1[kColumn];i++) {
	int kRow = index1[i];
	double value = element1[i];
	kRow = newRow[kRow];
	while (kRow>=0) {
	  index_[size_] = kRow;
	  element_[size_++] = value;
	  kRow = duplicateRow[kRow];
	}
      }
      start_[iColumn+1] = size_;
      length_[iColumn] = size_ - start_[iColumn];
    }
    delete [] newRow;
    delete [] duplicateRow;
  }
}


//-----------------------------------------------------------------------------

CoinPackedMatrix::~CoinPackedMatrix ()
{
   gutsOfDestructor();
}

//#############################################################################
//#############################################################################
//#############################################################################

void
CoinPackedMatrix::gutsOfDestructor()
{
   delete[] length_;
   delete[] start_;
   delete[] index_;
   delete[] element_;
   length_ = 0;
   start_ = 0;
   index_ = 0;
   element_ = 0;
}

//#############################################################################

void
CoinPackedMatrix::gutsOfCopyOf(const bool colordered,
			      const int minor, const int major,
			      const CoinBigIndex numels,
			      const double * elem, const int * ind,
			      const CoinBigIndex * start, const int * len,
			      const double extraMajor, const double extraGap)
{
   colOrdered_ = colordered;
   majorDim_ = major;
   minorDim_ = minor;
   size_ = numels;

   extraGap_ = extraGap;
   extraMajor_ = extraMajor;

   maxMajorDim_ = CoinLengthWithExtra(majorDim_, extraMajor_);

   if (maxMajorDim_ > 0) {
     delete [] length_;
     length_ = new int[maxMajorDim_];
     if (len == 0) {
       std::adjacent_difference(start + 1, start + (major + 1), length_);
     } else {
       CoinMemcpyN(len, major, length_);
     }
     delete [] start_;
     start_ = new CoinBigIndex[maxMajorDim_+1];
     start_[0]=0;
     CoinMemcpyN(start, major+1, start_);
   } else {
     // empty but be safe
     delete [] length_;
     length_ = NULL;
     delete [] start_;
     start_ = new CoinBigIndex[1];
     start_[0]=0;
   }

   maxSize_ = maxMajorDim_ > 0 ? start_[major] : 0;
   maxSize_ = CoinLengthWithExtra(maxSize_, extraMajor_);

   if (maxSize_ > 0) {
     delete [] element_;
     delete []index_;
     element_ = new double[maxSize_];
     index_ = new int[maxSize_];
     // we can't just simply memcpy these content over, because that can
     // upset memory debuggers like purify if there were gaps and those gaps
     // were uninitialized memory blocks
     for (int i = majorDim_ - 1; i >= 0; --i) {
       CoinMemcpyN(ind + start[i], length_[i], index_ + start_[i]);
       CoinMemcpyN(elem + start[i], length_[i], element_ + start_[i]);
     }
   }
}

//#############################################################################

void
CoinPackedMatrix::gutsOfCopyOfNoGaps(const bool colordered,
			      const int minor, const int major,
			      const double * elem, const int * ind,
                                     const CoinBigIndex * start)
{
   colOrdered_ = colordered;
   majorDim_ = major;
   minorDim_ = minor;
   size_ = start[majorDim_];

   extraGap_ = 0;
   extraMajor_ = 0;

   maxMajorDim_ = majorDim_;

   // delete all arrays
   delete [] length_;
   delete [] start_;
   delete [] element_;
   delete [] index_;
   
   if (maxMajorDim_ > 0) {
     length_ = new int[maxMajorDim_];
     assert (!start[0]);
     start_ = new CoinBigIndex[maxMajorDim_+1];
     start_[0]=0;
     CoinBigIndex last = 0;
     for (int i=0;i<majorDim_;i++) {
       CoinBigIndex first = last;
       last = start[i+1];
       length_[i] = last-first;
       start_[i+1]=last;
     }
   } else {
     // empty but be safe
     length_ = NULL;
     start_ = new CoinBigIndex[1];
     start_[0]=0;
   }

   maxSize_ = start_[majorDim_];

   if (maxSize_ > 0) {
     element_ = new double[maxSize_];
     index_ = new int[maxSize_];
     CoinMemcpyN(ind , maxSize_, index_);
     CoinMemcpyN(elem , maxSize_, element_);
   } else {
     element_ = NULL;
     index_ = NULL;
   }
}

//#############################################################################

void
CoinPackedMatrix::gutsOfOpEqual(const bool colordered,
			       const int minor, const int major,
			       const CoinBigIndex numels,
			       const double * elem, const int * ind,
			       const CoinBigIndex * start, const int * len)
{
   colOrdered_ = colordered;
   majorDim_ = major;
   minorDim_ = minor;
   size_ = numels;

   maxMajorDim_ = CoinLengthWithExtra(majorDim_, extraMajor_);

   int i;
   if (maxMajorDim_ > 0) {
     delete [] length_;
      length_ = new int[maxMajorDim_];
      if (len == 0) {
	 std::adjacent_difference(start + 1, start + (major + 1), length_);
      } else {
	 CoinMemcpyN(len, major, length_);
      }
      delete [] start_;
      start_ = new CoinBigIndex[maxMajorDim_+1];
      start_[0] = 0;
      if (extraGap_ == 0) {
	 for (i = 0; i < major; ++i)
	    start_[i+1] = start_[i] + length_[i];
      } else {
	const double extra_gap = extraGap_;
	for (i = 0; i < major; ++i)
	  start_[i+1] = start_[i] + CoinLengthWithExtra(length_[i], extra_gap);
      }
   } else {
     // empty matrix
     delete [] start_;
     start_ = new CoinBigIndex[1];
     start_[0] = 0;
   }

   maxSize_ = maxMajorDim_ > 0 ? start_[major] : 0;
   maxSize_ = CoinLengthWithExtra(maxSize_, extraMajor_);

   if (maxSize_ > 0) {
     delete [] element_;
     delete [] index_;
     element_ = new double[maxSize_];
     index_ = new int[maxSize_];
     assert (maxSize_>=start_[majorDim_-1]+length_[majorDim_-1]);
     // we can't just simply memcpy these content over, because that can
     // upset memory debuggers like purify if there were gaps and those gaps
     // were uninitialized memory blocks
     for (i = majorDim_ - 1; i >= 0; --i) {
       CoinMemcpyN(ind + start[i], length_[i], index_ + start_[i]);
       CoinMemcpyN(elem + start[i], length_[i], element_ + start_[i]);
     }
   }
#ifndef NDEBUG
   for ( i = majorDim_ - 1; i >= 0; --i) {
     const CoinBigIndex last = getVectorLast(i);
     for (CoinBigIndex j = getVectorFirst(i); j < last; ++j) {
       int index = index_[j];
       assert (index>=0&&index<minorDim_);
     }
   }
#endif
}

//#############################################################################

// This routine is called only if we MUST resize!
void
CoinPackedMatrix::resizeForAddingMajorVectors(const int numVec,
					     const int * lengthVec)
{
  const double extra_gap = extraGap_;
  int i;

  maxMajorDim_ =
    CoinMax(maxMajorDim_, CoinLengthWithExtra(majorDim_ + numVec, extraMajor_));

  CoinBigIndex * newStart = new CoinBigIndex[maxMajorDim_ + 1];
  int * newLength = new int[maxMajorDim_];

  CoinMemcpyN(length_, majorDim_, newLength);
  // fake that the new vectors are there
  CoinMemcpyN(lengthVec, numVec, newLength + majorDim_);
  majorDim_ += numVec;

  newStart[0] = 0;
  if (extra_gap == 0) {
    for (i = 0; i < majorDim_; ++i)
      newStart[i+1] = newStart[i] + newLength[i];
  } else {
    for (i = 0; i < majorDim_; ++i)
      newStart[i+1] = newStart[i] + CoinLengthWithExtra(newLength[i],extra_gap);
  }

  maxSize_ =
    CoinMax(maxSize_,newStart[majorDim_]+ (int) extraMajor_);
  majorDim_ -= numVec;

  int * newIndex = new int[maxSize_];
  double * newElem = new double[maxSize_];
  for (i = majorDim_ - 1; i >= 0; --i) {
    CoinMemcpyN(index_ + start_[i], length_[i], newIndex + newStart[i]);
    CoinMemcpyN(element_ + start_[i], length_[i], newElem + newStart[i]);
  }

  gutsOfDestructor();
  start_   = newStart;
  length_  = newLength;
  index_   = newIndex;
  element_ = newElem;
}

//#############################################################################

void
CoinPackedMatrix::resizeForAddingMinorVectors(const int * addedEntries)
{
   int i;
   maxMajorDim_ =
     CoinMax(CoinLengthWithExtra(majorDim_, extraMajor_), maxMajorDim_);
   CoinBigIndex * newStart = new CoinBigIndex[maxMajorDim_ + 1];
   int * newLength = new int[maxMajorDim_];
   // increase the lengths temporarily so that the correct new start positions
   // can be easily computed (it's faster to modify the lengths and reset them
   // than do a test for every entry when the start positions are computed.
   for (i = majorDim_ - 1; i >= 0; --i)
      newLength[i] = length_[i] + addedEntries[i];

   newStart[0] = 0;
   if (extraGap_ == 0) {
      for (i = 0; i < majorDim_; ++i)
	 newStart[i+1] = newStart[i] + newLength[i];
   } else {
      const double eg = extraGap_;
      for (i = 0; i < majorDim_; ++i)
	 newStart[i+1] = newStart[i] + CoinLengthWithExtra(newLength[i], eg);
   }

   // reset the lengths
   for (i = majorDim_ - 1; i >= 0; --i)
      newLength[i] -= addedEntries[i];

   maxSize_ =
     CoinMax(newStart[majorDim_]+ (int) extraMajor_, maxSize_);
   int * newIndex = new int[maxSize_];
   double * newElem = new double[maxSize_];
   for (i = majorDim_ - 1; i >= 0; --i) {
      CoinMemcpyN(index_ + start_[i], length_[i],
			newIndex + newStart[i]);
      CoinMemcpyN(element_ + start_[i], length_[i],
			newElem + newStart[i]);
   }

   gutsOfDestructor();
   start_   = newStart;
   length_  = newLength;
   index_   = newIndex;
   element_ = newElem;
}

//#############################################################################
//#############################################################################

void
CoinPackedMatrix::dumpMatrix(const char* fname) const
{
  if (! fname) {
    printf("Dumping matrix...\n\n");
    printf("colordered: %i\n", isColOrdered() ? 1 : 0);
    const int major = getMajorDim();
    const int minor = getMinorDim();
    printf("major: %i   minor: %i\n", major, minor);
    for (int i = 0; i < major; ++i) {
      printf("vec %i has length %i with entries:\n", i, length_[i]);
      for (CoinBigIndex j = start_[i]; j < start_[i] + length_[i]; ++j) {
	printf("        %15i  %40.25f\n", index_[j], element_[j]);
      }
    }
    printf("\nFinished dumping matrix\n");
  } else {
    FILE* out = fopen(fname, "w");
    fprintf(out, "Dumping matrix...\n\n");
    fprintf(out, "colordered: %i\n", isColOrdered() ? 1 : 0);
    const int major = getMajorDim();
    const int minor = getMinorDim();
    fprintf(out, "major: %i   minor: %i\n", major, minor);
    for (int i = 0; i < major; ++i) {
      fprintf(out, "vec %i has length %i with entries:\n", i, length_[i]);
      for (CoinBigIndex j = start_[i]; j < start_[i] + length_[i]; ++j) {
	fprintf(out, "        %15i  %40.25f\n", index_[j], element_[j]);
      }
    }
    fprintf(out, "\nFinished dumping matrix\n");
    fclose(out);
  }
}
void
CoinPackedMatrix::printMatrixElement(const int row_val, const int col_val) const
{
int major_index, minor_index;
if (isColOrdered())
  {
  major_index = col_val;
  minor_index = row_val;
  }
else
  {
  major_index = row_val;
  minor_index = col_val;
  }
if (getMajorDim() > major_index)
  std::cout << "Major index out of range: " << major_index << " vs. " << getMajorDim() << "\n";
CoinBigIndex curr_point = start_[major_index];
CoinBigIndex stop_point = curr_point + length_[major_index];
while (curr_point < stop_point)
  {
  if (index_[curr_point] == minor_index)
    {
      std::cout << element_[curr_point];
    return;
    }
  }
 std::cout << 0.0;
}
#ifndef CLP_NO_VECTOR
bool 
CoinPackedMatrix::isEquivalent2(const CoinPackedMatrix& rhs) const
{
  CoinRelFltEq eq;
  // Both must be column order or both row ordered and must be of same size
  if (isColOrdered() ^ rhs.isColOrdered()) {
    std::cerr<<"Ordering "<<isColOrdered()<<
      " rhs - "<<rhs.isColOrdered()<<std::endl;
    return false;
  }
  if (getNumCols() != rhs.getNumCols()) {
    std::cerr<<"NumCols "<<getNumCols()<<
      " rhs - "<<rhs.getNumCols()<<std::endl;
    return false;
  }
  if  (getNumRows() != rhs.getNumRows()) {
    std::cerr<<"NumRows "<<getNumRows()<<
      " rhs - "<<rhs.getNumRows()<<std::endl;
    return false;
  }
  if  (getNumElements() != rhs.getNumElements()) {
    std::cerr<<"NumElements "<<getNumElements()<<
      " rhs - "<<rhs.getNumElements()<<std::endl;
    return false;
  }
  
  for (int i=getMajorDim()-1; i >= 0; --i) {
    CoinShallowPackedVector pv = getVector(i);
    CoinShallowPackedVector rhsPv = rhs.getVector(i);
    if ( !pv.isEquivalent(rhsPv,eq) ) {
      std::cerr<<"vector # "<<i<<" nel "<<pv.getNumElements()<<
      " rhs - "<<rhsPv.getNumElements()<<std::endl;
      int j;
      const int * inds = pv.getIndices();
      const double * elems = pv.getElements();
      const int * inds2 = rhsPv.getIndices();
      const double * elems2 = rhsPv.getElements();
      for ( j = 0 ;j < pv.getNumElements() ;  ++j) {
	double diff = elems[j]-elems2[j];
	if (diff) {
	  std::cerr<<j<<"( "<<inds[j]<<", "<<elems[j]<<"), rhs ( "<<
	    inds2[j]<<", "<<elems2[j]<<") diff "<<
	    diff<<std::endl;
	  const int * xx = (const int *) (elems+j);
	  printf("%x %x",xx[0],xx[1]);
	  xx = (const int *) (elems2+j);
	  printf(" %x %x\n",xx[0],xx[1]);
	}
      }
      //return false;
    }
  }
  return true;
}
#else
/* Equivalence.
   Two matrices are equivalent if they are both by rows or both by columns,
   they have the same dimensions, and each vector is equivalent. 
   In this method the FloatEqual function operator can be specified. 
*/
bool 
CoinPackedMatrix::isEquivalent(const CoinPackedMatrix& rhs, const CoinRelFltEq& eq) const
{
  // Both must be column order or both row ordered and must be of same size
  if ((isColOrdered() ^ rhs.isColOrdered()) ||
      (getNumCols() != rhs.getNumCols()) ||
      (getNumRows() != rhs.getNumRows()) ||
      (getNumElements() != rhs.getNumElements()))
    return false;
  
  const int major = getMajorDim();
  const int minor = getMinorDim();
  double * values = new double[minor];
  memset(values,0,minor*sizeof(double));
  bool same=true;
  for (int i = 0; i < major; ++i) {
    int length = length_[i];
    if (length!=rhs.length_[i]) {
      same=false;
      break;
    } else {
      CoinBigIndex j;
      for ( j = start_[i]; j < start_[i] + length; ++j) {
        int index = index_[j];
        values[index]=element_[j];
      }
      for ( j = rhs.start_[i]; j < rhs.start_[i] + length; ++j) {
        int index = index_[j];
        double oldValue = values[index];
        values[index]=0.0;
        if (!eq(oldValue,rhs.element_[j])) {
          same=false;
          break;
        }
      }
      if (!same)
        break;
    }
  }
  delete [] values;
  return same;
}
#endif
/* Sort all columns so indices are increasing.in each column */
void 
CoinPackedMatrix::orderMatrix()
{
  for (int i=0;i<majorDim_;i++) {
    CoinBigIndex start = start_[i];
    CoinBigIndex end = start + length_[i];
    CoinSort_2(index_+start,index_+end,element_+start);
  }
}
/* 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 
   This version is easy one i.e. adding columns to column ordered */
int 
CoinPackedMatrix::appendMajor(const int number,
                              const CoinBigIndex * starts, const int * index,
                              const double * element, int numberOther)
{
  int i;
  int numberErrors=0;
  CoinBigIndex numberElements = starts[number];
  if (majorDim_ + number > maxMajorDim_ ||
      getLastStart() + numberElements > maxSize_) {
    // we got to resize before we add. note that the resizing method
    // properly fills out start_ and length_ for the major-dimension
    // vectors to be added!
    int * length = new int[number];
    for (i=0;i<number;i++)
      length[i]=starts[i+1]-starts[i];
    resizeForAddingMajorVectors(number, length);
    delete [] length;
    if (numberOther>0) {
      char * which = new char[numberOther];
      memset(which,0,numberOther);
      for (i = 0; i < number; i++) {
        CoinBigIndex put = start_[majorDim_+i];
        CoinBigIndex j;
        for ( j=starts[i];j<starts[i+1];j++) {
          int iIndex = index[j];
          element_[put]=element[j];
          if (iIndex>=0&&iIndex<numberOther) {
            if (!which[iIndex])
              which[iIndex]=1;
            else
              numberErrors++;
          } else {
            numberErrors++;
          }
          index_[put++]=iIndex;
        }
        for ( j=starts[i];j<starts[i+1];j++) {
          int iIndex = index[j];
          if (iIndex>=0&&iIndex<numberOther) 
            which[iIndex]=0;
        }
      }
      delete [] which;
    } else {
      // easy
      int lastMinor=-1;
      if (!extraGap_) {
        // just one copy
        int * index2 = index_+start_[majorDim_];
        for (CoinBigIndex j=0;j<numberElements;j++) {
          int iIndex = index[j];
          index2[j] = iIndex;
          lastMinor = CoinMax(lastMinor,iIndex);
        }
        CoinMemcpyN(element,numberElements,element_+start_[majorDim_]);
      } else {
        start_ += majorDim_;
        for (i = 0; i < number; i++) {
          int length = starts[i+1]-starts[i];
          int * index2 = index_+start_[i];
          const int * index1 = index+starts[i];
          for (CoinBigIndex j=0;j<length;j++) {
            int iIndex = index1[j];
            index2[j] = iIndex;
            lastMinor = CoinMax(lastMinor,iIndex);
          }
          CoinMemcpyN(element + starts[i], length,
                      element_ + start_[i]);
        }
        start_ -= majorDim_;
      }
      // update minorDim if necessary
      minorDim_ = CoinMax(minorDim_,lastMinor+1);
    }
  } else {
    if (numberOther>0) {
      char * which = new char[numberOther];
      memset(which,0,numberOther);
      for (i = 0; i < number; i++) {
        CoinBigIndex put = start_[majorDim_+i];
        CoinBigIndex j;
        for ( j=starts[i];j<starts[i+1];j++) {
          int iIndex = index[j];
          element_[put]=element[j];
          if (iIndex>=0&&iIndex<numberOther) {
            if (!which[iIndex])
              which[iIndex]=1;
            else
              numberErrors++;
          } else {
            numberErrors++;
          }
          index_[put++]=iIndex;
        }
        start_[majorDim_+i+1] = put;
        length_[majorDim_+i] = put-start_[majorDim_+i];;
        for ( j=starts[i];j<starts[i+1];j++) {
          int iIndex = index[j];
          if (iIndex>=0&&iIndex<numberOther) 
            which[iIndex]=0;
        }
      }
      delete [] which;
    } else {
      // easy
      int lastMinor=-1;
      if (!extraGap_) {
        // just one copy
        // just one copy
        int * index2 = index_+start_[majorDim_];
        for (CoinBigIndex j=0;j<numberElements;j++) {
          int iIndex = index[j];
          index2[j] = iIndex;
          lastMinor = CoinMax(lastMinor,iIndex);
        }
        CoinMemcpyN(element,numberElements,element_+start_[majorDim_]);
        start_ += majorDim_;
        for (i = 0; i < number; i++) {
          int length = starts[i+1]-starts[i];
          start_[i+1] = start_[i] + length;
          length_[majorDim_+i] = length;
        }
        start_ -= majorDim_;
      } else {
        start_ += majorDim_;
        for (i = 0; i < number; i++) {
          int length = starts[i+1]-starts[i];
          int * index2 = index_+start_[i];
          const int * index1 = index+starts[i];
          for (CoinBigIndex j=0;j<length;j++) {
            int iIndex = index1[j];
            index2[j] = iIndex;
            lastMinor = CoinMax(lastMinor,iIndex);
          }
          CoinMemcpyN(element + starts[i], length,
                      element_ + start_[i]);
          start_[i+1] = start_[i] + length;
          length_[majorDim_+i] = length;
        }
        start_ -= majorDim_;
      }
      // update minorDim if necessary
      minorDim_ = CoinMax(minorDim_,lastMinor+1);
    }
  }
  majorDim_ += number;
  size_ += numberElements;
  return numberErrors;
}
/* 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
   This version is harder one i.e. adding columns to row ordered */
int 
CoinPackedMatrix::appendMinor(const int number,
                              const CoinBigIndex * starts, const int * index,
                              const double * element, int numberOther)
{
  int i;
  int numberErrors=0;
  // first compute how many entries will be added to each major-dimension
  // vector, and if needed, resize the matrix to accommodate all
  int * addedEntries = NULL;
  if (numberOther>0) {
    addedEntries = new int[majorDim_];
    CoinZeroN(addedEntries,majorDim_);
    numberOther=majorDim_;
    char * which = new char[numberOther];
    memset(which,0,numberOther);
    for (i = 0; i < number; i++) {
      CoinBigIndex j;
      for ( j=starts[i];j<starts[i+1];j++) {
        int iIndex = index[j];
        if (iIndex>=0&&iIndex<numberOther) {
          addedEntries[iIndex]++;
          if (!which[iIndex])
            which[iIndex]=1;
          else
            numberErrors++;
        } else {
          numberErrors++;
        }
      }
      for ( j=starts[i];j<starts[i+1];j++) {
        int iIndex = index[j];
        if (iIndex>=0&&iIndex<numberOther) 
          which[iIndex]=0;
      }
    }
    delete [] which;
  } else {
    int largest = majorDim_-1;
    for (i = 0; i < number; i++) {
      CoinBigIndex j;
      for ( j=starts[i];j<starts[i+1];j++) {
        int iIndex = index[j];
        largest = CoinMax(largest,iIndex);
      }
    }
    if (largest+1>majorDim_) {
      if (isColOrdered())
        setDimensions(-1,largest+1);
      else 
        setDimensions(largest+1,-1);
    }
    addedEntries = new int[majorDim_];
    CoinZeroN(addedEntries,majorDim_);
    // no checking
    for (i = 0; i < number; i++) {
      CoinBigIndex j;
      for ( j=starts[i];j<starts[i+1];j++) {
        int iIndex = index[j];
        addedEntries[iIndex]++;
      }
    }
  }
  for (i = majorDim_ - 1; i >= 0; i--) {
    if (start_[i] + length_[i] + addedEntries[i] > start_[i+1])
      break;
  }
  if (i >= 0)
    resizeForAddingMinorVectors(addedEntries);
  delete[] addedEntries;

  // now insert the entries of matrix
  for (i = 0; i < number; i++) {
    CoinBigIndex j;
    for ( j=starts[i];j<starts[i+1];j++) {
      int iIndex = index[j];
      element_[start_[iIndex] + length_[iIndex]] = element[j];
      index_[start_[iIndex] + (length_[iIndex]++)] = minorDim_;
    }
    ++minorDim_;
  }
  size_ += starts[number];
  return numberErrors;
}