/* \$Id: AbcMatrix.hpp 1910 2013-01-27 02:00:13Z forrest \$ */ // 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 `y + A * scalar *x` in `y`. @pre `x` must be of size `numColumns()` @pre `y` must be of size `numRows()` */ void timesModifyExcludingSlacks(double scalar, const double * x, double * y) const; /** Return `y + A * scalar(+-1) *x` in `y`. @pre `x` must be of size `numColumns()+numRows()` @pre `y` must be of size `numRows()` */ void timesModifyIncludingSlacks(double scalar, const double * x, double * y) const; /** Return `A * scalar(+-1) *x` in `y`. @pre `x` must be of size `numColumns()+numRows()` @pre `y` must be of size `numRows()` */ void timesIncludingSlacks(double scalar, const double * x, double * y) const; /** Return A * scalar(+-1) *x + y in `y`. @pre `x` must be of size `numRows()` @pre `y` must be of size `numRows()+numColumns()` */ void transposeTimesNonBasic(double scalar, const double * x, double * y) const; /** Return y - A * x in `y`. @pre `x` must be of size `numRows()` @pre `y` must be of size `numRows()+numColumns()` */ void transposeTimesAll(const double * x, double * y) const; /** Return y + A * scalar(+-1) *x in `y`. @pre `x` must be of size `numRows()` @pre `y` must be of size `numRows()` */ void transposeTimesBasic(double scalar, const double * x, double * y) const; /** Return ```x * scalar * A/code> in z. 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 x *A in z but just for indices Already in z. Note - z always packed mode */ void subsetTransposeTimes(const CoinIndexedVector & x, CoinIndexedVector & z) const; /// Return -x *A in z 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 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 x * -1 * A in z. 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 x * -1 * A in z. 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 ```