/* $Id$ */
// 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
/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/