/* $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). /* Authors John Forrest */ /** This deals with Factorization and Updates I am assuming that 32 bits is enough for number of rows or columns, but CoinBigIndex may be redefined to get 64 bits. */ #include "AbcCommon.hpp" #include "CoinAbcHelperFunctions.hpp" #if ABC_PARALLEL #define FACTOR_CPU 4 #else #define FACTOR_CPU 1 #endif #define LARGE_SET COIN_INT_MAX - 10 #define LARGE_UNSET (LARGE_SET + 1) class CoinAbcTypeFactorization : public CoinAbcAnyFactorization { friend void CoinAbcFactorizationUnitTest(const std::string &mpsDir); public: /**@name Constructors and destructor and copy */ //@{ /// Default constructor CoinAbcTypeFactorization(); /// Copy constructor CoinAbcTypeFactorization(const CoinAbcTypeFactorization &other); /// Copy constructor CoinAbcTypeFactorization(const CoinFactorization &other); /// Destructor virtual ~CoinAbcTypeFactorization(); /// Clone virtual CoinAbcAnyFactorization *clone() const; /// Delete all stuff (leaves as after CoinAbcFactorization()) void almostDestructor(); /// Debug show object (shows one representation) void show_self() const; /// Debug - sort so can compare void sort() const; /// = copy CoinAbcTypeFactorization &operator=(const CoinAbcTypeFactorization &other); //@} /**@name Do factorization */ //@{ /// Condition number - product of pivots after factorization CoinSimplexDouble conditionNumber() const; //@} /**@name general stuff such as permutation or status */ //@{ /// Returns address of permute region inline CoinSimplexInt *permute() const { return NULL; //permute_.array(); } /// Returns array to put basis indices in virtual inline CoinSimplexInt *indices() const { return indexRowU_.array(); } /// Returns address of pivotColumn region (also used for permuting) virtual inline CoinSimplexInt *pivotColumn() const { return pivotColumn_.array(); } /// Returns address of pivot region virtual inline CoinFactorizationDouble *pivotRegion() const { return pivotRegionAddress_; } #if ABC_SMALL < 2 /// Start of each row in L inline CoinBigIndex *startRowL() const { return startRowL_.array(); } #endif /// Start of each column in L inline CoinBigIndex *startColumnL() const { return startColumnL_.array(); } #if ABC_SMALL < 2 /// Index of column in row for L inline CoinSimplexInt *indexColumnL() const { return indexColumnL_.array(); } #endif /// Row indices of L inline CoinSimplexInt *indexRowL() const { return indexRowL_.array(); } #if ABC_SMALL < 2 /// Elements in L (row copy) inline CoinFactorizationDouble *elementByRowL() const { return elementByRowL_.array(); } #endif /** Forward and backward linked lists (numberRows_+2) **/ inline CoinSimplexInt *pivotLinkedBackwards() const { return firstCount_.array() + numberRows_ + 1; } inline CoinSimplexInt *pivotLinkedForwards() const { return firstCount_.array() + 2 * numberRows_ + 3; } inline CoinSimplexInt *pivotLOrder() const { return firstCount_.array(); } #if ABC_SMALL < 0 #define ABC_USE_FUNCTION_POINTERS 0 #define SMALL_PERMUTE #endif #ifdef ABC_USE_FUNCTION_POINTERS typedef void (*scatterUpdate)(int, CoinFactorizationDouble, const CoinFactorizationDouble *, CoinFactorizationDouble *); #if ABC_USE_FUNCTION_POINTERS typedef struct { scatterUpdate functionPointer; CoinBigIndex offset; int number; } scatterStruct; #else typedef struct { CoinBigIndex offset; int number; } scatterStruct; #endif /// Array of function pointers PLUS for U Column inline scatterStruct *scatterUColumn() const { return scatterPointersUColumnAddress_; } #endif /// For equal counts in factorization /** First Row/Column with count of k, can tell which by offset - Rows then Columns actually comes before nextCount*/ inline CoinSimplexInt *firstCount() const { return firstCount_.array(); } /// Next Row/Column with count inline CoinSimplexInt *nextCount() const { return firstCount_.array() + numberRows_ + 2; } /// Previous Row/Column with count inline CoinSimplexInt *lastCount() const { return firstCount_.array() + 3 * numberRows_ + 2; } /// Number of Rows after iterating inline CoinSimplexInt numberRowsExtra() const { return numberRowsExtra_; } /// Number in L inline CoinBigIndex numberL() const { return numberL_; } /// Base of L inline CoinBigIndex baseL() const { return baseL_; } /// Maximum of Rows after iterating inline CoinSimplexInt maximumRowsExtra() const { return maximumRowsExtra_; } /// Total number of elements in factorization virtual inline CoinBigIndex numberElements() const { return totalElements_; } /// Length of FT vector inline CoinSimplexInt numberForrestTomlin() const { return numberInColumn_.array()[numberRowsExtra_]; } /// Returns areaFactor but adjusted for dense CoinSimplexDouble adjustedAreaFactor() const; /// Level of detail of messages inline CoinSimplexInt messageLevel() const { return messageLevel_; } void messageLevel(CoinSimplexInt value); /// Set maximum pivots virtual void maximumPivots(CoinSimplexInt value); #if ABC_SMALL < 4 /// Gets dense threshold inline CoinSimplexInt denseThreshold() const { return denseThreshold_; } /// Sets dense threshold inline void setDenseThreshold(CoinSimplexInt value) { denseThreshold_ = value; } #endif /// Returns maximum absolute value in factorization CoinSimplexDouble maximumCoefficient() const; #if 0 /// true if Forrest Tomlin update, false if PFI inline bool forrestTomlin() const { return doForrestTomlin_;} inline void setForrestTomlin(bool value) { doForrestTomlin_=value;} #endif /// True if FT update and space inline bool spaceForForrestTomlin() const { CoinBigIndex start = lastEntryByColumnU_; CoinBigIndex space = lengthAreaU_ - (start + numberRowsExtra_); return (space >= 0); //&&doForrestTomlin_; } //@} /**@name some simple stuff */ //@{ /// Returns number in U area inline CoinBigIndex numberElementsU() const { return lengthU_; } /// Setss number in U area inline void setNumberElementsU(CoinBigIndex value) { lengthU_ = value; } /// Returns length of U area inline CoinBigIndex lengthAreaU() const { return lengthAreaU_; } /// Returns number in L area inline CoinBigIndex numberElementsL() const { return lengthL_; } /// Returns length of L area inline CoinBigIndex lengthAreaL() const { return lengthAreaL_; } /// Returns number in R area inline CoinBigIndex numberElementsR() const { return lengthR_; } /// Number of compressions done inline CoinBigIndex numberCompressions() const { return numberCompressions_; } /// Returns pivot row //virtual CoinSimplexInt * pivotRow() const; /// Returns work area //virtual CoinFactorizationDouble * workArea() const; /// Returns CoinSimplexInt work area //virtual CoinSimplexInt * intWorkArea() const; /// Returns array to put basis starts in virtual inline CoinBigIndex *starts() const { return startColumnU_.array(); } /// Number of entries in each row virtual inline CoinSimplexInt *numberInRow() const { return numberInRow_.array(); } /// Number of entries in each column virtual inline CoinSimplexInt *numberInColumn() const { return numberInColumn_.array(); } /// Returns array to put basis elements in virtual inline CoinFactorizationDouble *elements() const { return elementU_.array(); } /// Start of columns for R inline CoinBigIndex *startColumnR() const { return reinterpret_cast< CoinBigIndex * >(firstCount_.array() + 3 * numberRows_ + 4); } /// Elements of U inline CoinFactorizationDouble *elementU() const { return elementU_.array(); } /// Row indices of U inline CoinSimplexInt *indexRowU() const { return indexRowU_.array(); } /// Start of each column in U inline CoinBigIndex *startColumnU() const { return startColumnU_.array(); } #if COIN_BIG_DOUBLE == 1 /// To a work array and associate vector void toLongArray(CoinIndexedVector *vector, int which) const; /// From a work array and dis-associate vector void fromLongArray(CoinIndexedVector *vector) const; /// From a work array and dis-associate vector void fromLongArray(int which) const; /// Returns long double * associated with vector long double *denseVector(CoinIndexedVector *vector) const; /// Returns long double * associated with vector long double *denseVector(CoinIndexedVector &vector) const; /// Returns long double * associated with vector const long double *denseVector(const CoinIndexedVector *vector) const; /// Returns long double * associated with vector const long double *denseVector(const CoinIndexedVector &vector) const; /// Scans region to find nonzeros void scan(CoinIndexedVector *vector) const; /// Clear all hidden arrays void clearHiddenArrays(); #else /// Returns double * associated with vector inline double *denseVector(CoinIndexedVector *vector) const { return vector->denseVector(); } inline double *denseVector(CoinIndexedVector &vector) const { return vector.denseVector(); } /// Returns double * associated with vector inline const double *denseVector(const CoinIndexedVector *vector) const { return vector->denseVector(); } inline const double *denseVector(const CoinIndexedVector &vector) const { return vector.denseVector(); } /// To a work array and associate vector inline void toLongArray(CoinIndexedVector *vector, int which) const {} /// From a work array and dis-associate vector inline void fromLongArray(CoinIndexedVector *vector) const {} /// From a work array and dis-associate vector inline void fromLongArray(int which) const {} /// Scans region to find nonzeros inline void scan(CoinIndexedVector *vector) const { vector->scan(0, numberRows_, zeroTolerance_); } #endif #ifdef ABC_ORDERED_FACTORIZATION /// Permute in for Ftran void permuteInForFtran(CoinIndexedVector ®ionSparse, bool full = false) const; /// Permute in for Btran and multiply by pivot Region void permuteInForBtranAndMultiply(CoinIndexedVector ®ionSparse, bool full = false) const; /// Permute out for Btran void permuteOutForBtran(CoinIndexedVector ®ionSparse) const; #endif /** Array persistence flag If 0 then as now (delete/new) 1 then only do arrays if bigger needed 2 as 1 but give a bit extra if bigger needed */ //inline CoinSimplexInt persistenceFlag() const //{ return persistenceFlag_;} //@} /**@name rank one updates which do exist */ //@{ #if 0 /** Checks if can replace one Column to basis, returns 0=OK, 1=Probably OK, 2=singular, 3=no room, 5 max pivots Fills in region for use later partial update already in U */ virtual int checkReplace ( const AbcSimplex * model, CoinIndexedVector * regionSparse, int pivotRow, CoinSimplexDouble & pivotCheck, double acceptablePivot = 1.0e-8); /** Replaces one Column to basis, returns 0=OK, 1=Probably OK, 2=singular, 3=no room If skipBtranU is false will do btran part partial update already in U */ virtual CoinSimplexInt replaceColumn ( CoinIndexedVector * regionSparse, CoinSimplexInt pivotRow, CoinSimplexDouble pivotCheck , bool skipBtranU=false, CoinSimplexDouble acceptablePivot=1.0e-8); #endif /** Checks if can replace one Column to basis, returns update alpha Fills in region for use later partial update already in U */ virtual #ifdef ABC_LONG_FACTORIZATION long #endif double checkReplacePart1(CoinIndexedVector *regionSparse, int pivotRow); /** Checks if can replace one Column to basis, returns update alpha Fills in region for use later partial update in vector */ virtual #ifdef ABC_LONG_FACTORIZATION long #endif double checkReplacePart1(CoinIndexedVector *regionSparse, CoinIndexedVector *partialUpdate, int pivotRow); #ifdef MOVE_REPLACE_PART1A /** Checks if can replace one Column to basis, returns update alpha Fills in region for use later partial update already in U */ virtual void checkReplacePart1a(CoinIndexedVector *regionSparse, int pivotRow); virtual #ifdef ABC_LONG_FACTORIZATION long #endif double checkReplacePart1b(CoinIndexedVector *regionSparse, int pivotRow); #endif /** Checks if can replace one Column to basis, returns 0=OK, 1=Probably OK, 2=singular, 3=no room, 5 max pivots */ virtual int checkReplacePart2(int pivotRow, CoinSimplexDouble btranAlpha, double ftranAlpha, #ifdef ABC_LONG_FACTORIZATION long #endif double ftAlpha, double acceptablePivot = 1.0e-8); /** Replaces one Column to basis, partial update already in U */ virtual void replaceColumnPart3(const AbcSimplex *model, CoinIndexedVector *regionSparse, CoinIndexedVector *tableauColumn, int pivotRow, #ifdef ABC_LONG_FACTORIZATION long #endif double alpha); /** Replaces one Column to basis, partial update in vector */ virtual void replaceColumnPart3(const AbcSimplex *model, CoinIndexedVector *regionSparse, CoinIndexedVector *tableauColumn, CoinIndexedVector *partialUpdate, int pivotRow, #ifdef ABC_LONG_FACTORIZATION long #endif double alpha); #ifdef EARLY_FACTORIZE /// 0 success, -1 can't +1 accuracy problems virtual int replaceColumns(const AbcSimplex *model, CoinIndexedVector &stuff, int firstPivot, int lastPivot, bool cleanUp); #endif /// Update partial Ftran by R update void updatePartialUpdate(CoinIndexedVector &partialUpdate); /// Returns true if wants tableauColumn in replaceColumn inline virtual bool wantsTableauColumn() const { return false; } /** Combines BtranU and store which elements are to be deleted returns number to be deleted */ int replaceColumnU(CoinIndexedVector *regionSparse, CoinBigIndex *deletedPosition, CoinSimplexInt *deletedColumns, CoinSimplexInt pivotRow); //@} /**@name various uses of factorization (return code number elements) which user may want to know about */ /// Later take out return codes (apart from +- 1 on FT) //@{ /** Updates one column (FTRAN) from regionSparse2 Tries to do FT update number returned is negative if no room regionSparse starts as zero and is zero at end. Note - if regionSparse2 packed on input - will be packed on output */ virtual CoinSimplexInt updateColumnFT(CoinIndexedVector ®ionSparse); virtual int updateColumnFTPart1(CoinIndexedVector ®ionSparse); virtual void updateColumnFTPart2(CoinIndexedVector ®ionSparse); /** Updates one column (FTRAN) Tries to do FT update puts partial update in vector */ virtual void updateColumnFT(CoinIndexedVector ®ionSparseFT, CoinIndexedVector &partialUpdate, int which); /** This version has same effect as above with FTUpdate==false so number returned is always >=0 */ virtual CoinSimplexInt updateColumn(CoinIndexedVector ®ionSparse) const; /** Updates one column (FTRAN) from region2 Tries to do FT update number returned is negative if no room. Also updates region3 region1 starts as zero and is zero at end */ virtual CoinSimplexInt updateTwoColumnsFT(CoinIndexedVector ®ionFT, CoinIndexedVector ®ionOther); /** Updates one column (BTRAN) from regionSparse2 regionSparse starts as zero and is zero at end Note - if regionSparse2 packed on input - will be packed on output */ virtual CoinSimplexInt updateColumnTranspose(CoinIndexedVector ®ionSparse) const; /** Updates one full column (FTRAN) */ virtual void updateFullColumn(CoinIndexedVector ®ionSparse) const; /** Updates one full column (BTRAN) */ virtual void updateFullColumnTranspose(CoinIndexedVector ®ionSparse) const; /** Updates one column for dual steepest edge weights (FTRAN) */ virtual void updateWeights(CoinIndexedVector ®ionSparse) const; /** Updates one column (FTRAN) */ virtual void updateColumnCpu(CoinIndexedVector ®ionSparse, int whichCpu) const; /** Updates one column (BTRAN) */ virtual void updateColumnTransposeCpu(CoinIndexedVector ®ionSparse, int whichCpu) const; void unpack(CoinIndexedVector *regionFrom, CoinIndexedVector *regionTo) const; void pack(CoinIndexedVector *regionFrom, CoinIndexedVector *regionTo) const; /** makes a row copy of L for speed and to allow very sparse problems */ inline void goSparse() {} void goSparse2(); #ifndef NDEBUG virtual void checkMarkArrays() const; #endif #if ABC_SMALL < 2 /** get sparse threshold */ inline CoinSimplexInt sparseThreshold() const { return sparseThreshold_; } #endif /** set sparse threshold */ void sparseThreshold(CoinSimplexInt value); //@} /// *** Below this user may not want to know about /**@name various uses of factorization (return code number elements) which user may not want to know about (left over from my LP code) */ //@{ /// Get rid of all memory inline void clearArrays() { gutsOfDestructor(); } //@} /**@name used by ClpFactorization */ /// See if worth going sparse void checkSparse(); /// The real work of constructors etc 0 just scalars, 1 bit normal void gutsOfDestructor(CoinSimplexInt type = 1); /// 1 bit - tolerances etc, 2 more, 4 dummy arrays void gutsOfInitialize(CoinSimplexInt type); void gutsOfCopy(const CoinAbcTypeFactorization &other); /// Reset all sparsity etc statistics void resetStatistics(); void printRegion(const CoinIndexedVector &vector, const char *where) const; //@} /**@name used by factorization */ /// Gets space for a factorization, called by constructors virtual void getAreas(CoinSimplexInt numberRows, CoinSimplexInt numberColumns, CoinBigIndex maximumL, CoinBigIndex maximumU); /// PreProcesses column ordered copy of basis virtual void preProcess(); void preProcess(CoinSimplexInt); /// Return largest element double preProcess3(); void preProcess4(); /// Does most of factorization virtual CoinSimplexInt factor(AbcSimplex *model); #ifdef EARLY_FACTORIZE /// Returns -2 if can't, -1 if singular, -99 memory, 0 OK virtual int factorize(AbcSimplex *model, CoinIndexedVector &stuff); #endif /// Does post processing on valid factorization - putting variables on correct rows virtual void postProcess(const CoinSimplexInt *sequence, CoinSimplexInt *pivotVariable); /// Makes a non-singular basis by replacing variables virtual void makeNonSingular(CoinSimplexInt *sequence); protected: /** Does sparse phase of factorization return code is <0 error, 0= finished */ CoinSimplexInt factorSparse(); /** Does dense phase of factorization return code is <0 error, 0= finished */ CoinSimplexInt factorDense(); /// Pivots when just one other row so faster? bool pivotOneOtherRow(CoinSimplexInt pivotRow, CoinSimplexInt pivotColumn); /// Does one pivot on Row Singleton in factorization bool pivotRowSingleton(CoinSimplexInt pivotRow, CoinSimplexInt pivotColumn); /// Does one pivot on Column Singleton in factorization (can't return false) void pivotColumnSingleton(CoinSimplexInt pivotRow, CoinSimplexInt pivotColumn); /// After pivoting void afterPivot(CoinSimplexInt pivotRow, CoinSimplexInt pivotColumn); /// After pivoting - returns true if need to go dense int wantToGoDense(); /** Gets space for one Column with given length, may have to do compression (returns True if successful), also moves existing vector, extraNeeded is over and above present */ bool getColumnSpace(CoinSimplexInt iColumn, CoinSimplexInt extraNeeded); /** Reorders U so contiguous and in order (if there is space) Returns true if it could */ bool reorderU(); /** getColumnSpaceIterateR. Gets space for one extra R element in Column may have to do compression (returns true) also moves existing vector */ bool getColumnSpaceIterateR(CoinSimplexInt iColumn, CoinFactorizationDouble value, CoinSimplexInt iRow); /** getColumnSpaceIterate. Gets space for one extra U element in Column may have to do compression (returns true) also moves existing vector. Returns -1 if no memory or where element was put Used by replaceRow (turns off R version) */ CoinBigIndex getColumnSpaceIterate(CoinSimplexInt iColumn, CoinFactorizationDouble value, CoinSimplexInt iRow); /** Gets space for one Row with given length, may have to do compression (returns True if successful), also moves existing vector */ bool getRowSpace(CoinSimplexInt iRow, CoinSimplexInt extraNeeded); /** Gets space for one Row with given length while iterating, may have to do compression (returns True if successful), also moves existing vector */ bool getRowSpaceIterate(CoinSimplexInt iRow, CoinSimplexInt extraNeeded); /// Checks that row and column copies look OK void checkConsistency(); //#define CHECK_LINKS #ifdef CHECK_LINKS void checkLinks(int x = 0); #else #define checkLinks(x) #endif /// Adds a link in chain of equal counts inline void addLink(CoinSimplexInt index, CoinSimplexInt count) { CoinSimplexInt *COIN_RESTRICT nextCount = nextCountAddress_; CoinSimplexInt *COIN_RESTRICT firstCount = this->firstCount(); CoinSimplexInt *COIN_RESTRICT lastCount = lastCountAddress_; CoinSimplexInt next = firstCount[count]; firstCount[count] = index; nextCount[index] = next; lastCount[index] = count - numberRows_ - 2; // points to firstCount[count] if (next >= 0) lastCount[next] = index; } /// Deletes a link in chain of equal counts inline void deleteLink(CoinSimplexInt index) { CoinSimplexInt *COIN_RESTRICT nextCount = nextCountAddress_; CoinSimplexInt *COIN_RESTRICT lastCount = lastCountAddress_; CoinSimplexInt next = nextCount[index]; CoinSimplexInt last = lastCount[index]; assert(next != index); assert(last != index); if (next >= 0) lastCount[next] = last; if (last >= 0) { nextCount[last] = next; } else { int count = last + numberRows_ + 2; CoinSimplexInt *COIN_RESTRICT firstCount = this->firstCount(); firstCount[count] = next; } } /// Modifies links in chain of equal counts inline void modifyLink(CoinSimplexInt index, CoinSimplexInt count) { CoinSimplexInt *COIN_RESTRICT nextCount = nextCountAddress_; CoinSimplexInt *COIN_RESTRICT lastCount = lastCountAddress_; CoinSimplexInt *COIN_RESTRICT firstCount = this->firstCount(); CoinSimplexInt next2 = firstCount[count]; if (next2 == index) return; firstCount[count] = index; CoinSimplexInt next = nextCount[index]; CoinSimplexInt last = lastCount[index]; assert(next != index); assert(last != index); nextCount[index] = next2; lastCount[index] = count - numberRows_ - 2; // points to firstCount[count] if (next >= 0) lastCount[next] = last; if (next2 >= 0) lastCount[next2] = index; if (last >= 0) { nextCount[last] = next; } else { int count = last + numberRows_ + 2; firstCount[count] = next; } } /// Separate out links with same row/column count void separateLinks(); void separateLinks(CoinSimplexInt, CoinSimplexInt); /// Cleans up at end of factorization void cleanup(); /// Set up addresses from arrays void doAddresses(); /// Updates part of column (FTRANL) void updateColumnL(CoinIndexedVector *region #if ABC_SMALL < 2 , CoinAbcStatistics &statistics #endif #if ABC_PARALLEL , int whichSparse = 0 #endif ) const; /// Updates part of column (FTRANL) when densish void updateColumnLDensish(CoinIndexedVector *region) const; /// Updates part of column (FTRANL) when dense (i.e. do as inner products) void updateColumnLDense(CoinIndexedVector *region) const; /// Updates part of column (FTRANL) when sparse void updateColumnLSparse(CoinIndexedVector *region #if ABC_PARALLEL , int whichSparse #endif ) const; /// Updates part of column (FTRANR) without FT update void updateColumnR(CoinIndexedVector *region #if ABC_SMALL < 2 , CoinAbcStatistics &statistics #endif #if ABC_PARALLEL , int whichSparse = 0 #endif ) const; /// Store update after doing L and R - retuns false if no room bool storeFT( #if ABC_SMALL < 3 const #endif CoinIndexedVector *regionFT); /// Updates part of column (FTRANU) void updateColumnU(CoinIndexedVector *region #if ABC_SMALL < 2 , CoinAbcStatistics &statistics #endif #if ABC_PARALLEL , int whichSparse = 0 #endif ) const; /// Updates part of column (FTRANU) when sparse void updateColumnUSparse(CoinIndexedVector *regionSparse #if ABC_PARALLEL , int whichSparse #endif ) const; /// Updates part of column (FTRANU) void updateColumnUDensish(CoinIndexedVector *regionSparse) const; /// Updates part of column (FTRANU) when dense (i.e. do as inner products) void updateColumnUDense(CoinIndexedVector *regionSparse) const; /// Updates part of 2 columns (FTRANU) real work void updateTwoColumnsUDensish( CoinSimplexInt &numberNonZero1, CoinFactorizationDouble *COIN_RESTRICT region1, CoinSimplexInt *COIN_RESTRICT index1, CoinSimplexInt &numberNonZero2, CoinFactorizationDouble *COIN_RESTRICT region2, CoinSimplexInt *COIN_RESTRICT index2) const; /// Updates part of column PFI (FTRAN) (after rest) void updateColumnPFI(CoinIndexedVector *regionSparse) const; /// Updates part of column transpose PFI (BTRAN) (before rest) void updateColumnTransposePFI(CoinIndexedVector *region) const; /** Updates part of column transpose (BTRANU), assumes index is sorted i.e. region is correct */ void updateColumnTransposeU(CoinIndexedVector *region, CoinSimplexInt smallestIndex #if ABC_SMALL < 2 , CoinAbcStatistics &statistics #endif #if ABC_PARALLEL , int whichCpu #endif ) const; /** Updates part of column transpose (BTRANU) when densish, assumes index is sorted i.e. region is correct */ void updateColumnTransposeUDensish(CoinIndexedVector *region, CoinSimplexInt smallestIndex) const; /** Updates part of column transpose (BTRANU) when sparse, assumes index is sorted i.e. region is correct */ void updateColumnTransposeUSparse(CoinIndexedVector *region #if ABC_PARALLEL , int whichSparse #endif ) const; /** Updates part of column transpose (BTRANU) by column assumes index is sorted i.e. region is correct */ void updateColumnTransposeUByColumn(CoinIndexedVector *region, CoinSimplexInt smallestIndex) const; /// Updates part of column transpose (BTRANR) void updateColumnTransposeR(CoinIndexedVector *region #if ABC_SMALL < 2 , CoinAbcStatistics &statistics #endif ) const; /// Updates part of column transpose (BTRANR) when dense void updateColumnTransposeRDensish(CoinIndexedVector *region) const; /// Updates part of column transpose (BTRANR) when sparse void updateColumnTransposeRSparse(CoinIndexedVector *region) const; /// Updates part of column transpose (BTRANL) void updateColumnTransposeL(CoinIndexedVector *region #if ABC_SMALL < 2 , CoinAbcStatistics &statistics #endif #if ABC_PARALLEL , int whichSparse #endif ) const; /// Updates part of column transpose (BTRANL) when densish by column void updateColumnTransposeLDensish(CoinIndexedVector *region) const; /// Updates part of column transpose (BTRANL) when densish by row void updateColumnTransposeLByRow(CoinIndexedVector *region) const; /// Updates part of column transpose (BTRANL) when sparse (by Row) void updateColumnTransposeLSparse(CoinIndexedVector *region #if ABC_PARALLEL , int whichSparse #endif ) const; public: /** Replaces one Column to basis for PFI returns 0=OK, 1=Probably OK, 2=singular, 3=no room. In this case region is not empty - it is incoming variable (updated) */ CoinSimplexInt replaceColumnPFI(CoinIndexedVector *regionSparse, CoinSimplexInt pivotRow, CoinSimplexDouble alpha); protected: /** Returns accuracy status of replaceColumn returns 0=OK, 1=Probably OK, 2=singular */ CoinSimplexInt checkPivot(CoinSimplexDouble saveFromU, CoinSimplexDouble oldPivot) const; /// 0 fine, -99 singular, 2 dense int pivot(CoinSimplexInt pivotRow, CoinSimplexInt pivotColumn, CoinBigIndex pivotRowPosition, CoinBigIndex pivotColumnPosition, CoinFactorizationDouble *COIN_RESTRICT work, CoinSimplexUnsignedInt *COIN_RESTRICT workArea2, CoinSimplexInt increment2, int *COIN_RESTRICT markRow); int pivot(CoinSimplexInt &pivotRow, CoinSimplexInt &pivotColumn, CoinBigIndex pivotRowPosition, CoinBigIndex pivotColumnPosition, int *COIN_RESTRICT markRow); #if ABC_SMALL < 2 #define CONVERTROW 2 #elif ABC_SMALL < 4 #else #undef ABC_DENSE_CODE #define ABC_DENSE_CODE 0 #endif //@} ////////////////// data ////////////////// protected: /**@name data */ //@{ CoinSimplexInt *pivotColumnAddress_; CoinSimplexInt *permuteAddress_; CoinFactorizationDouble *pivotRegionAddress_; CoinFactorizationDouble *elementUAddress_; CoinSimplexInt *indexRowUAddress_; CoinSimplexInt *numberInColumnAddress_; CoinSimplexInt *numberInColumnPlusAddress_; #ifdef ABC_USE_FUNCTION_POINTERS /// Array of function pointers scatterStruct *scatterPointersUColumnAddress_; CoinFactorizationDouble *elementUColumnPlusAddress_; #endif CoinBigIndex *startColumnUAddress_; #if CONVERTROW CoinBigIndex *convertRowToColumnUAddress_; #if CONVERTROW > 1 CoinBigIndex *convertColumnToRowUAddress_; #endif #endif #if ABC_SMALL < 2 CoinFactorizationDouble *elementRowUAddress_; #endif CoinBigIndex *startRowUAddress_; CoinSimplexInt *numberInRowAddress_; CoinSimplexInt *indexColumnUAddress_; CoinSimplexInt *firstCountAddress_; /// Next Row/Column with count CoinSimplexInt *nextCountAddress_; /// Previous Row/Column with count CoinSimplexInt *lastCountAddress_; CoinSimplexInt *nextColumnAddress_; CoinSimplexInt *lastColumnAddress_; CoinSimplexInt *nextRowAddress_; CoinSimplexInt *lastRowAddress_; CoinSimplexInt *saveColumnAddress_; //CoinSimplexInt * saveColumnAddress2_; CoinCheckZero *markRowAddress_; CoinSimplexInt *listAddress_; CoinFactorizationDouble *elementLAddress_; CoinSimplexInt *indexRowLAddress_; CoinBigIndex *startColumnLAddress_; #if ABC_SMALL < 2 CoinBigIndex *startRowLAddress_; #endif CoinSimplexInt *pivotLinkedBackwardsAddress_; CoinSimplexInt *pivotLinkedForwardsAddress_; CoinSimplexInt *pivotLOrderAddress_; CoinBigIndex *startColumnRAddress_; /// Elements of R CoinFactorizationDouble *elementRAddress_; /// Row indices for R CoinSimplexInt *indexRowRAddress_; CoinSimplexInt *indexColumnLAddress_; CoinFactorizationDouble *elementByRowLAddress_; #if ABC_SMALL < 4 CoinFactorizationDouble *denseAreaAddress_; #endif CoinFactorizationDouble *workAreaAddress_; CoinSimplexUnsignedInt *workArea2Address_; mutable CoinSimplexInt *sparseAddress_; #ifdef SMALL_PERMUTE CoinSimplexInt *fromSmallToBigRow_; CoinSimplexInt *fromSmallToBigColumn_; #endif /// Number of Rows after iterating CoinSimplexInt numberRowsExtra_; /// Maximum number of Rows after iterating CoinSimplexInt maximumRowsExtra_; /// Size of small inverse CoinSimplexInt numberRowsSmall_; /// Number factorized in L CoinSimplexInt numberGoodL_; /// Number Rows left (numberRows-numberGood) CoinSimplexInt numberRowsLeft_; /// Number of elements in U (to go) /// or while iterating total overall CoinBigIndex totalElements_; /// First place in funny copy zeroed out CoinBigIndex firstZeroed_; #if ABC_SMALL < 2 /// Below this use sparse technology - if 0 then no L row copy CoinSimplexInt sparseThreshold_; #endif /// Number in R CoinSimplexInt numberR_; /// Length of R stuff CoinBigIndex lengthR_; /// length of area reserved for R CoinBigIndex lengthAreaR_; /// Number in L CoinBigIndex numberL_; /// Base of L CoinBigIndex baseL_; /// Length of L CoinBigIndex lengthL_; /// Length of area reserved for L CoinBigIndex lengthAreaL_; /// Number in U CoinSimplexInt numberU_; /// Maximum space used in U CoinBigIndex maximumU_; /// Length of U CoinBigIndex lengthU_; /// Length of area reserved for U CoinBigIndex lengthAreaU_; /// Last entry by column for U CoinBigIndex lastEntryByColumnU_; #ifdef ABC_USE_FUNCTION_POINTERS /// Last entry by column for U CoinBigIndex lastEntryByColumnUPlus_; /// Length of U CoinBigIndex lengthAreaUPlus_; #endif /// Last entry by row for U CoinBigIndex lastEntryByRowU_; /// Number of trials before rejection CoinSimplexInt numberTrials_; #if ABC_SMALL < 4 /// Leading dimension for dense CoinSimplexInt leadingDimension_; #endif #if COIN_BIG_DOUBLE == 1 /// Work arrays mutable CoinFactorizationLongDoubleArrayWithLength longArray_[FACTOR_CPU]; /// Associated CoinIndexedVector mutable CoinIndexedVector *associatedVector_[FACTOR_CPU]; #endif /// Pivot order for each Column CoinIntArrayWithLength pivotColumn_; /// Permutation vector for pivot row order CoinIntArrayWithLength permute_; /// Start of each Row as pointer CoinBigIndexArrayWithLength startRowU_; /// Number in each Row CoinIntArrayWithLength numberInRow_; /// Number in each Column CoinIntArrayWithLength numberInColumn_; /// Number in each Column including pivoted CoinIntArrayWithLength numberInColumnPlus_; /** First Row/Column with count of k, can tell which by offset - Rows then Columns */ CoinIntArrayWithLength firstCount_; /// Next Column in memory order CoinIntArrayWithLength nextColumn_; /// Previous Column in memory order CoinIntArrayWithLength lastColumn_; /// Next Row in memory order CoinIntArrayWithLength nextRow_; /// Previous Row in memory order CoinIntArrayWithLength lastRow_; /// Columns left to do in a single pivot CoinIntArrayWithLength saveColumn_; /// Marks rows to be updated CoinIntArrayWithLength markRow_; /// Base address for U (may change) CoinIntArrayWithLength indexColumnU_; /// Inverses of pivot values CoinFactorizationDoubleArrayWithLength pivotRegion_; /// Elements of U CoinFactorizationDoubleArrayWithLength elementU_; /// Row indices of U CoinIntArrayWithLength indexRowU_; /// Start of each column in U CoinBigIndexArrayWithLength startColumnU_; #ifdef ABC_USE_FUNCTION_POINTERS /// Array of structs for U Column CoinArbitraryArrayWithLength scatterUColumn_; #endif #if CONVERTROW /// Converts rows to columns in U CoinBigIndexArrayWithLength convertRowToColumnU_; #if CONVERTROW > 1 /// Converts columns to rows in U CoinBigIndexArrayWithLength convertColumnToRowU_; #endif #endif #if ABC_SMALL < 2 /// Elements of U by row CoinFactorizationDoubleArrayWithLength elementRowU_; #endif /// Elements of L CoinFactorizationDoubleArrayWithLength elementL_; /// Row indices of L CoinIntArrayWithLength indexRowL_; /// Start of each column in L CoinBigIndexArrayWithLength startColumnL_; #if ABC_SMALL < 4 /// Dense area CoinFactorizationDoubleArrayWithLength denseArea_; #endif /// First work area CoinFactorizationDoubleArrayWithLength workArea_; /// Second work area CoinUnsignedIntArrayWithLength workArea2_; #if ABC_SMALL < 2 /// Start of each row in L CoinBigIndexArrayWithLength startRowL_; /// Index of column in row for L CoinIntArrayWithLength indexColumnL_; /// Elements in L (row copy) CoinFactorizationDoubleArrayWithLength elementByRowL_; /// Sparse regions mutable CoinIntArrayWithLength sparse_; #endif /// Detail in messages CoinSimplexInt messageLevel_; /// Number of compressions done CoinBigIndex numberCompressions_; // last slack pivot row CoinSimplexInt lastSlack_; #if ABC_SMALL < 2 /// To decide how to solve mutable double ftranCountInput_; mutable double ftranCountAfterL_; mutable double ftranCountAfterR_; mutable double ftranCountAfterU_; double ftranAverageAfterL_; double ftranAverageAfterR_; double ftranAverageAfterU_; #if FACTORIZATION_STATISTICS double ftranTwiddleFactor1_; double ftranTwiddleFactor2_; #endif mutable CoinSimplexInt numberFtranCounts_; #endif /// Maximum rows (ever) (here to use double alignment) CoinSimplexInt maximumRows_; #if ABC_SMALL < 2 mutable double ftranFTCountInput_; mutable double ftranFTCountAfterL_; mutable double ftranFTCountAfterR_; mutable double ftranFTCountAfterU_; double ftranFTAverageAfterL_; double ftranFTAverageAfterR_; double ftranFTAverageAfterU_; #if FACTORIZATION_STATISTICS double ftranFTTwiddleFactor1_; double ftranFTTwiddleFactor2_; #endif mutable CoinSimplexInt numberFtranFTCounts_; #endif #if ABC_SMALL < 4 /// Dense threshold (here to use double alignment) CoinSimplexInt denseThreshold_; #endif #if ABC_SMALL < 2 mutable double btranCountInput_; mutable double btranCountAfterU_; mutable double btranCountAfterR_; mutable double btranCountAfterL_; double btranAverageAfterU_; double btranAverageAfterR_; double btranAverageAfterL_; #if FACTORIZATION_STATISTICS double btranTwiddleFactor1_; double btranTwiddleFactor2_; #endif mutable CoinSimplexInt numberBtranCounts_; #endif /// Maximum maximum pivots CoinSimplexInt maximumMaximumPivots_; #if ABC_SMALL < 2 /// To decide how to solve mutable double ftranFullCountInput_; mutable double ftranFullCountAfterL_; mutable double ftranFullCountAfterR_; mutable double ftranFullCountAfterU_; double ftranFullAverageAfterL_; double ftranFullAverageAfterR_; double ftranFullAverageAfterU_; #if FACTORIZATION_STATISTICS double ftranFullTwiddleFactor1_; double ftranFullTwiddleFactor2_; #endif mutable CoinSimplexInt numberFtranFullCounts_; #endif /// Rows first time nonzero CoinSimplexInt initialNumberRows_; #if ABC_SMALL < 2 /// To decide how to solve mutable double btranFullCountInput_; mutable double btranFullCountAfterL_; mutable double btranFullCountAfterR_; mutable double btranFullCountAfterU_; double btranFullAverageAfterL_; double btranFullAverageAfterR_; double btranFullAverageAfterU_; #if FACTORIZATION_STATISTICS double btranFullTwiddleFactor1_; double btranFullTwiddleFactor2_; #endif mutable CoinSimplexInt numberBtranFullCounts_; #endif /** State of saved version and what can be done 0 - nothing saved 1 - saved and can go back to previous save by unwinding 2 - saved - getting on for a full copy higher bits - see ABC_FAC.... */ CoinSimplexInt state_; /// Size in bytes of a sparseArray CoinBigIndex sizeSparseArray_; public: #if ABC_SMALL < 2 #if ABC_SMALL >= 0 inline bool gotLCopy() const { return ((state_ & ABC_FAC_GOT_LCOPY) != 0); } inline void setNoGotLCopy() { state_ &= ~ABC_FAC_GOT_LCOPY; } inline void setYesGotLCopy() { state_ |= ABC_FAC_GOT_LCOPY; } inline bool gotRCopy() const { return ((state_ & ABC_FAC_GOT_RCOPY) != 0); } inline void setNoGotRCopy() { state_ &= ~ABC_FAC_GOT_RCOPY; } inline void setYesGotRCopy() { state_ |= ABC_FAC_GOT_RCOPY; } inline bool gotUCopy() const { return ((state_ & ABC_FAC_GOT_UCOPY) != 0); } inline void setNoGotUCopy() { state_ &= ~ABC_FAC_GOT_UCOPY; } inline void setYesGotUCopy() { state_ |= ABC_FAC_GOT_UCOPY; } inline bool gotSparse() const { return ((state_ & ABC_FAC_GOT_SPARSE) != 0); } inline void setNoGotSparse() { state_ &= ~ABC_FAC_GOT_SPARSE; } inline void setYesGotSparse() { state_ |= ABC_FAC_GOT_SPARSE; } #else // force use of copies inline bool gotLCopy() const { return true; } inline void setNoGotLCopy() {} inline void setYesGotLCopy() {} inline bool gotRCopy() const { return true; } inline void setNoGotRCopy() {} inline void setYesGotRCopy() {} inline bool gotUCopy() const { return true; } inline void setNoGotUCopy() {} inline void setYesGotUCopy() {} inline bool gotSparse() const { return true; } inline void setNoGotSparse() {} inline void setYesGotSparse() {} #endif #else // force no use of copies inline bool gotLCopy() const { return false; } inline void setNoGotLCopy() {} inline void setYesGotLCopy() {} inline bool gotRCopy() const { return false; } inline void setNoGotRCopy() {} inline void setYesGotRCopy() {} inline bool gotUCopy() const { return false; } inline void setNoGotUCopy() {} inline void setYesGotUCopy() {} inline bool gotSparse() const { return false; } inline void setNoGotSparse() {} inline void setYesGotSparse() {} #endif /** Array persistence flag If 0 then as now (delete/new) 1 then only do arrays if bigger needed 2 as 1 but give a bit extra if bigger needed */ //CoinSimplexInt persistenceFlag_; //@} }; /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 */