/* $Id: AbcMatrix.hpp 1910 2013-01-27 02:00:13Z stefan $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#ifndef AbcMatrix_H
#define AbcMatrix_H

#include "CoinPragma.hpp"

#include "ClpMatrixBase.hpp"
#include "AbcSimplex.hpp"
#include "CoinAbcHelperFunctions.hpp"
/** This implements a scaled version of CoinPackedMatrix
    It may have THREE! copies
    1) scaled CoinPackedMatrix without gaps
    2) row copy non-basic,basic, fixed
    3) vector copy
*/
class AbcMatrix2;
class AbcMatrix3;
class AbcMatrix {
  
public:
  /**@name Useful methods */
  //@{
  /// Return a complete CoinPackedMatrix
  inline CoinPackedMatrix * getPackedMatrix() const {
    return matrix_;
  }
  /** Whether the packed matrix is column major ordered or not. */
  inline bool isColOrdered() const {
    return true;
  }
  /** Number of entries in the packed matrix. */
  inline CoinBigIndex getNumElements() const {
    return matrix_->getNumElements();
  }
  /** Number of columns. */
  inline int getNumCols() const {
    assert(matrix_->getNumCols()==model_->numberColumns());return matrix_->getNumCols();
  }
  /** Number of rows. */
  inline int getNumRows() const {
    assert(matrix_->getNumRows()==model_->numberRows());return matrix_->getNumRows();
  }
  /// Sets model
  void setModel(AbcSimplex * model);
  /// A vector containing the elements in the packed matrix.
  inline const double * getElements() const {
    return matrix_->getElements();
  }
  /// Mutable elements
  inline double * getMutableElements() const {
    return matrix_->getMutableElements();
  }
  /// A vector containing the minor indices of the elements in the packed matrix. 
  inline const int * getIndices() const {
    return matrix_->getIndices();
  }
  /// A vector containing the minor indices of the elements in the packed matrix. 
  inline int * getMutableIndices() const {
    return matrix_->getMutableIndices();
  }
  /// Starts
  inline const CoinBigIndex * getVectorStarts() const {
    return matrix_->getVectorStarts();
  }
  inline  CoinBigIndex * getMutableVectorStarts() const {
    return matrix_->getMutableVectorStarts();
  }
  /** The lengths of the major-dimension vectors. */
  inline const int * getVectorLengths() const {
    return matrix_->getVectorLengths();
  }
  /** The lengths of the major-dimension vectors. */
  inline int * getMutableVectorLengths() const {
    return matrix_->getMutableVectorLengths();
  }
  /// Row starts
  CoinBigIndex * rowStart() const;
  /// Row ends
  CoinBigIndex * rowEnd() const;
  /// Row elements
  double * rowElements() const;
  /// Row columns
  CoinSimplexInt * rowColumns() const;
  /** Returns a new matrix in reverse order without gaps */
  CoinPackedMatrix * reverseOrderedCopy() const;
  /// Returns number of elements in column part of basis
  CoinBigIndex countBasis(const int * whichColumn,
			  int & numberColumnBasic);
  /// Fills in column part of basis
  void fillBasis(const int * whichColumn,
		 int & numberColumnBasic,
		 int * row, int * start,
		 int * rowCount, int * columnCount,
		 CoinSimplexDouble * element);
  /// Fills in column part of basis
  void fillBasis(const int * whichColumn,
		 int & numberColumnBasic,
		 int * row, int * start,
		 int * rowCount, int * columnCount,
		 long double * element);
  /** Scales and creates row copy
  */
  void scale(int numberRowsAlreadyScaled);
  /// Creates row copy
  void createRowCopy();
  /// Take out of useful
  void takeOutOfUseful(int sequence,CoinIndexedVector & spare);
  /// Put into useful
  void putIntofUseful(int sequence,CoinIndexedVector & spare);
  /// Put in and out for useful
  void inOutUseful(int sequenceIn,int sequenceOut);
  /// Make all useful
  void makeAllUseful(CoinIndexedVector & spare);
  /// Sort into useful
  void sortUseful(CoinIndexedVector & spare);
  /// Move largest in column to beginning (not used as doesn't help factorization)
  void moveLargestToStart();
  
  /** Unpacks a column into an CoinIndexedVector
   */
  void unpack(CoinIndexedVector & rowArray,
	      int column) const ;
  /** Adds multiple of a column (or slack) into an CoinIndexedvector
      You can use quickAdd to add to vector */
  void add(CoinIndexedVector & rowArray, int column, double multiplier) const ;
  //@}
  
  /**@name Matrix times vector methods */
  //@{
  /** Return <code>y + A * scalar *x</code> in <code>y</code>.
      @pre <code>x</code> must be of size <code>numColumns()</code>
      @pre <code>y</code> must be of size <code>numRows()</code> */
  void timesModifyExcludingSlacks(double scalar,
	     const double * x, double * y) const;
  /** Return <code>y + A * scalar(+-1) *x</code> in <code>y</code>.
      @pre <code>x</code> must be of size <code>numColumns()+numRows()</code>
      @pre <code>y</code> must be of size <code>numRows()</code> */
  void timesModifyIncludingSlacks(double scalar,
	     const double * x, double * y) const;
  /** Return <code>A * scalar(+-1) *x</code> in <code>y</code>.
      @pre <code>x</code> must be of size <code>numColumns()+numRows()</code>
      @pre <code>y</code> must be of size <code>numRows()</code> */
  void timesIncludingSlacks(double scalar,
	     const double * x, double * y) const;
  /** Return A * scalar(+-1) *x + y</code> in <code>y</code>.
      @pre <code>x</code> must be of size <code>numRows()</code>
      @pre <code>y</code> must be of size <code>numRows()+numColumns()</code> */
  void transposeTimesNonBasic(double scalar,
			    const double * x, double * y) const;
  /** Return y - A * x</code> in <code>y</code>.
      @pre <code>x</code> must be of size <code>numRows()</code>
      @pre <code>y</code> must be of size <code>numRows()+numColumns()</code> */
  void transposeTimesAll(const double * x, double * y) const;
  /** Return y + A * scalar(+-1) *x</code> in <code>y</code>.
      @pre <code>x</code> must be of size <code>numRows()</code>
      @pre <code>y</code> must be of size <code>numRows()</code> */
  void transposeTimesBasic(double scalar,
			    const double * x, double * y) const;
  /** Return <code>x * scalar * A/code> in <code>z</code>.
      Note - x unpacked mode - z packed mode including slacks
      All these return atLo/atUp first then free/superbasic
      number of first set returned
      pivotVariable is extended to have that order
      reversePivotVariable used to update that list
      free/superbasic only stored in normal format
      can use spare array to get this effect
      may put djs alongside atLo/atUp
      Squashes small elements and knows about AbcSimplex */
  int transposeTimesNonBasic(double scalar,
		      const CoinIndexedVector & x,
		      CoinIndexedVector & z) const;
  /// gets sorted tableau row and a possible value of theta
  double dualColumn1(const CoinIndexedVector & update,
		     CoinPartitionedVector & tableauRow,
		     CoinPartitionedVector & candidateList) const;
  /// gets sorted tableau row and a possible value of theta
  double dualColumn1Row(int iBlock, double upperThetaSlack, int & freeSequence,
			const CoinIndexedVector & update,
			CoinPartitionedVector & tableauRow,
			CoinPartitionedVector & candidateList) const;
  /// gets sorted tableau row and a possible value of theta
  double dualColumn1RowFew(int iBlock, double upperThetaSlack, int & freeSequence,
			 const CoinIndexedVector & update,
			 CoinPartitionedVector & tableauRow,
			 CoinPartitionedVector & candidateList) const;
  /// gets sorted tableau row and a possible value of theta
  double dualColumn1Row2(double upperThetaSlack, int & freeSequence,
			 const CoinIndexedVector & update,
			 CoinPartitionedVector & tableauRow,
			 CoinPartitionedVector & candidateList) const;
  /// gets sorted tableau row and a possible value of theta
  double dualColumn1Row1(double upperThetaSlack, int & freeSequence,
			 const CoinIndexedVector & update,
			 CoinPartitionedVector & tableauRow,
			 CoinPartitionedVector & candidateList) const;
  /** gets sorted tableau row and a possible value of theta
      On input first,last give what to scan 
      On output is number in tableauRow and candidateList */
  void dualColumn1Part(int iBlock,int & sequenceIn, double & upperTheta,
			 const CoinIndexedVector & update,
			 CoinPartitionedVector & tableauRow,
			 CoinPartitionedVector & candidateList) const;
  /// rebalance for parallel
  void rebalance() const;
  /// Get sequenceIn when Dantzig
  int pivotColumnDantzig(const CoinIndexedVector & updates,
			 CoinPartitionedVector & spare) const; 
  /// Get sequenceIn when Dantzig (One block)
  int pivotColumnDantzig(int iBlock,bool doByRow,const CoinIndexedVector & updates,
			 CoinPartitionedVector & spare,
			 double & bestValue) const; 
  /// gets tableau row - returns number of slacks in block 
  int primalColumnRow(int iBlock,bool doByRow,const CoinIndexedVector & update,
			  CoinPartitionedVector & tableauRow) const;
  /// gets tableau row and dj row - returns number of slacks in block 
  int primalColumnRowAndDjs(int iBlock,const CoinIndexedVector & updateTableau,
			  const CoinIndexedVector & updateDjs,
			    CoinPartitionedVector & tableauRow) const;
  /** Chooses best weighted dj
   */
  int chooseBestDj(int iBlock,const CoinIndexedVector & infeasibilities, 
		   const double * weights) const;
  /** does steepest edge double or triple update
      If scaleFactor!=0 then use with tableau row to update djs
      otherwise use updateForDjs
      Returns best sequence
   */
  int primalColumnDouble(int iBlock,CoinPartitionedVector & updateForTableauRow,
			  CoinPartitionedVector & updateForDjs,
			  const CoinIndexedVector & updateForWeights,
			  CoinPartitionedVector & spareColumn1,
			  double * infeasibilities, 
			  double referenceIn, double devex,
			  // Array for exact devex to say what is in reference framework
			  unsigned int * reference,
			  double * weights, double scaleFactor) const;
  /** does steepest edge double or triple update
      If scaleFactor!=0 then use with tableau row to update djs
      otherwise use updateForDjs
      Returns best sequence
   */
  int primalColumnSparseDouble(int iBlock,CoinPartitionedVector & updateForTableauRow,
			  CoinPartitionedVector & updateForDjs,
			  const CoinIndexedVector & updateForWeights,
			  CoinPartitionedVector & spareColumn1,
			  double * infeasibilities, 
			  double referenceIn, double devex,
			  // Array for exact devex to say what is in reference framework
			  unsigned int * reference,
			  double * weights, double scaleFactor) const;
  /** does steepest edge double or triple update
      If scaleFactor!=0 then use with tableau row to update djs
      otherwise use updateForDjs
      Returns best sequence
   */
  int primalColumnDouble(CoinPartitionedVector & updateForTableauRow,
			 CoinPartitionedVector & updateForDjs,
			 const CoinIndexedVector & updateForWeights,
			 CoinPartitionedVector & spareColumn1,
			 CoinIndexedVector & infeasible, 
			 double referenceIn, double devex,
			 // Array for exact devex to say what is in reference framework
			 unsigned int * reference,
			 double * weights, double scaleFactor) const;
  /// gets subset updates
  void primalColumnSubset(int iBlock,const CoinIndexedVector & update,
			  const CoinPartitionedVector & tableauRow,
			  CoinPartitionedVector & weights) const;
  /// Partial pricing
  void partialPricing(double startFraction, double endFraction,
		      int & bestSequence, int & numberWanted);
  /** Return <code>x *A</code> in <code>z</code> but
      just for indices Already in z.
      Note - z always packed mode */
  void subsetTransposeTimes(const CoinIndexedVector & x,
			    CoinIndexedVector & z) const;
  /// Return <code>-x *A</code> in <code>z</code>
  void transposeTimes(const CoinIndexedVector & x,
			    CoinIndexedVector & z) const;
  //@}
  
  /**@name Other */
  //@{
  /// Returns CoinPackedMatrix (non const)
  inline CoinPackedMatrix * matrix() const {
    return matrix_;
  }
  /** Partial pricing tuning parameter - minimum number of "objects" to scan.
      e.g. number of Gub sets but could be number of variables */
  inline int minimumObjectsScan() const {
    return minimumObjectsScan_;
  }
  inline void setMinimumObjectsScan(int value) {
    minimumObjectsScan_ = value;
  }
  /// Partial pricing tuning parameter - minimum number of negative reduced costs to get
  inline int minimumGoodReducedCosts() const {
    return minimumGoodReducedCosts_;
  }
  inline void setMinimumGoodReducedCosts(int value) {
    minimumGoodReducedCosts_ = value;
  }
  /// Current start of search space in matrix (as fraction)
  inline double startFraction() const {
    return startFraction_;
  }
  inline void setStartFraction(double value) {
    startFraction_ = value;
  }
  /// Current end of search space in matrix (as fraction)
  inline double endFraction() const {
    return endFraction_;
  }
  inline void setEndFraction(double value) {
    endFraction_ = value;
  }
  /// Current best reduced cost
  inline double savedBestDj() const {
    return savedBestDj_;
  }
  inline void setSavedBestDj(double value) {
    savedBestDj_ = value;
  }
  /// Initial number of negative reduced costs wanted
  inline int originalWanted() const {
    return originalWanted_;
  }
  inline void setOriginalWanted(int value) {
    originalWanted_ = value;
  }
  /// Current number of negative reduced costs which we still need
  inline int currentWanted() const {
    return currentWanted_;
  }
  inline void setCurrentWanted(int value) {
    currentWanted_ = value;
  }
  /// Current best sequence
  inline int savedBestSequence() const {
    return savedBestSequence_;
  }
  inline void setSavedBestSequence(int value) {
    savedBestSequence_ = value;
  }
  /// Start of each column block
  inline int * startColumnBlock() const
  {return startColumnBlock_;}
  /// Start of each block (in stored)
  inline const int * blockStart() const
  { return blockStart_;}
  inline bool gotRowCopy() const
  { return rowStart_!=0;}
  /// Start of each block (in stored)
  inline int blockStart(int block) const
  { return blockStart_[block];}
  /// Number of actual column blocks
  inline int numberColumnBlocks() const
  { return numberColumnBlocks_;}
  /// Number of actual row blocks
  inline int numberRowBlocks() const
  { return numberRowBlocks_;}
  //@}
  
  
  /**@name Constructors, destructor */
  //@{
  /** Default constructor. */
  AbcMatrix();
  /** Destructor */
  ~AbcMatrix();
  //@}
  
  /**@name Copy method */
  //@{
  /** The copy constructor. */
  AbcMatrix(const AbcMatrix&);
  /** The copy constructor from an CoinPackedMatrix. */
  AbcMatrix(const CoinPackedMatrix&);
  /** Subset constructor (without gaps).  Duplicates are allowed
      and order is as given */
  AbcMatrix (const AbcMatrix & wholeModel,
	     int numberRows, const int * whichRows,
	     int numberColumns, const int * whichColumns);
  AbcMatrix (const CoinPackedMatrix & wholeModel,
	     int numberRows, const int * whichRows,
	     int numberColumns, const int * whichColumns);
  
  AbcMatrix& operator=(const AbcMatrix&);
  /// Copy contents - resizing if necessary - otherwise re-use memory
  void copy(const AbcMatrix * from);
  //@}
private:
  
protected:
  /**@name Data members
     The data members are protected to allow access for derived classes. */
  //@{
  /// Data
  CoinPackedMatrix * matrix_;
  /// Model
  mutable AbcSimplex * model_;
#if ABC_PARALLEL==0
#define NUMBER_ROW_BLOCKS 1
#define NUMBER_COLUMN_BLOCKS 1
#elif ABC_PARALLEL==1
#define NUMBER_ROW_BLOCKS 4
#define NUMBER_COLUMN_BLOCKS 4
#else
#define NUMBER_ROW_BLOCKS 8
#define NUMBER_COLUMN_BLOCKS 8
#endif
  /** Start of each row (per block) - last lot are useless
      first all row starts for block 0, then for block2
      so NUMBER_ROW_BLOCKS+2 times number rows */
  CoinBigIndex * rowStart_;
  /// Values by row
  double * element_;
  /// Columns
  int * column_;
  /// Start of each column block
  mutable int startColumnBlock_[NUMBER_COLUMN_BLOCKS+1];
  /// Start of each block (in stored)
  int blockStart_[NUMBER_ROW_BLOCKS+1];
  /// Number of actual column blocks
  mutable int numberColumnBlocks_;
  /// Number of actual row blocks
  int numberRowBlocks_;
  //#define COUNT_COPY
#ifdef COUNT_COPY
#define MAX_COUNT 13
  /// Start in elements etc
  CoinBigIndex countStart_[MAX_COUNT+1];
  /// First column
  int countFirst_[MAX_COUNT+1];
  // later int countEndUseful_[MAX_COUNT+1];
  int * countRealColumn_;
  // later int * countInverseRealColumn_;
  CoinBigIndex * countStartLarge_;
  int * countRow_;
  double * countElement_;
  int smallestCount_;
  int largestCount_;
#endif 
  /// Special row copy
  //AbcMatrix2 * rowCopy_;
  /// Special column copy
  //AbcMatrix3 * columnCopy_;
  /// Current start of search space in matrix (as fraction)
  double startFraction_;
  /// Current end of search space in matrix (as fraction)
  double endFraction_;
  /// Best reduced cost so far
  double savedBestDj_;
  /// Initial number of negative reduced costs wanted
  int originalWanted_;
  /// Current number of negative reduced costs which we still need
  int currentWanted_;
  /// Saved best sequence in pricing
  int savedBestSequence_;
  /// Partial pricing tuning parameter - minimum number of "objects" to scan
  int minimumObjectsScan_;
  /// Partial pricing tuning parameter - minimum number of negative reduced costs to get
  int minimumGoodReducedCosts_;
  //@}
};
#ifdef THREAD
#include <pthread.h>
typedef struct {
  double acceptablePivot;
  const AbcSimplex * model;
  double * spare;
  int * spareIndex;
  double * arrayTemp;
  int * indexTemp;
  int * numberInPtr;
  double * bestPossiblePtr;
  double * upperThetaPtr;
  int * posFreePtr;
  double * freePivotPtr;
  int * numberOutPtr;
  const unsigned short * count;
  const double * pi;
  const CoinBigIndex * rowStart;
  const double * element;
  const unsigned short * column;
  int offset;
  int numberInRowArray;
  int numberLook;
} dualColumn0Struct;
#endif
class AbcMatrix2 {
  
public:
  /**@name Useful methods */
  //@{
  /** Return <code>x * -1 * A in <code>z</code>.
      Note - x packed and z will be packed mode
      Squashes small elements and knows about AbcSimplex */
  void transposeTimes(const AbcSimplex * model,
		      const CoinPackedMatrix * rowCopy,
		      const CoinIndexedVector & x,
		      CoinIndexedVector & spareArray,
		      CoinIndexedVector & z) const;
  /// Returns true if copy has useful information
  inline bool usefulInfo() const {
    return rowStart_ != NULL;
  }
  //@}
  
  
  /**@name Constructors, destructor */
  //@{
  /** Default constructor. */
  AbcMatrix2();
  /** Constructor from copy. */
  AbcMatrix2(AbcSimplex * model, const CoinPackedMatrix * rowCopy);
  /** Destructor */
  ~AbcMatrix2();
  //@}
  
  /**@name Copy method */
  //@{
  /** The copy constructor. */
  AbcMatrix2(const AbcMatrix2&);
  AbcMatrix2& operator=(const AbcMatrix2&);
  //@}
  
  
protected:
  /**@name Data members
     The data members are protected to allow access for derived classes. */
  //@{
  /// Number of blocks
  int numberBlocks_;
  /// Number of rows
  int numberRows_;
  /// Column offset for each block (plus one at end)
  int * offset_;
  /// Counts of elements in each part of row
  mutable unsigned short * count_;
  /// Row starts
  mutable CoinBigIndex * rowStart_;
  /// columns within block
  unsigned short * column_;
  /// work arrays
  double * work_;
#ifdef THREAD
  pthread_t * threadId_;
  dualColumn0Struct * info_;
#endif
  //@}
};
typedef struct {
  CoinBigIndex startElements_; // point to data
  int startIndices_; // point to column_
  int numberInBlock_;
  int numberPrice_; // at beginning
  int numberElements_; // number elements per column
} blockStruct3;
class AbcMatrix3 {
  
public:
  /**@name Useful methods */
  //@{
  /** Return <code>x * -1 * A in <code>z</code>.
      Note - x packed and z will be packed mode
      Squashes small elements and knows about AbcSimplex */
  void transposeTimes(const AbcSimplex * model,
		      const double * pi,
		      CoinIndexedVector & output) const;
  /// Updates two arrays for steepest
  void transposeTimes2(const AbcSimplex * model,
		       const double * pi, CoinIndexedVector & dj1,
		       const double * piWeight,
		       double referenceIn, double devex,
		       // Array for exact devex to say what is in reference framework
		       unsigned int * reference,
		       double * weights, double scaleFactor);
  //@}
  
  
  /**@name Constructors, destructor */
  //@{
  /** Default constructor. */
  AbcMatrix3();
  /** Constructor from copy. */
  AbcMatrix3(AbcSimplex * model, const CoinPackedMatrix * columnCopy);
  /** Destructor */
  ~AbcMatrix3();
  //@}
  
  /**@name Copy method */
  //@{
  /** The copy constructor. */
  AbcMatrix3(const AbcMatrix3&);
  AbcMatrix3& operator=(const AbcMatrix3&);
  //@}
  /**@name Sort methods */
  //@{
  /** Sort blocks */
  void sortBlocks(const AbcSimplex * model);
  /// Swap one variable
  void swapOne(const AbcSimplex * model, const AbcMatrix * matrix,
	       int iColumn);
  //@}
  
  
protected:
  /**@name Data members
     The data members are protected to allow access for derived classes. */
  //@{
  /// Number of blocks
  int numberBlocks_;
  /// Number of columns
  int numberColumns_;
  /// Column indices and reverse lookup (within block)
  int * column_;
  /// Starts for odd/long vectors
  CoinBigIndex * start_;
  /// Rows
  int * row_;
  /// Elements
  double * element_;
  /// Blocks (ordinary start at 0 and go to first block)
  blockStruct * block_;
  //@}
};

#endif