00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef CoinFactorization_H
00013 #define CoinFactorization_H
00014
00015
00016 #include <iostream>
00017 #include <string>
00018 #include <cassert>
00019 #include <cstdio>
00020 #include <cmath>
00021 #include "CoinTypes.hpp"
00022 #include "CoinIndexedVector.hpp"
00023
00024 class CoinPackedMatrix;
00050 class CoinFactorization {
00051 friend void CoinFactorizationUnitTest( const std::string & mpsDir );
00052
00053 public:
00054
00057
00058 CoinFactorization ( );
00060 CoinFactorization ( const CoinFactorization &other);
00061
00063 ~CoinFactorization ( );
00065 void almostDestructor();
00067 void show_self ( ) const;
00069 int saveFactorization (const char * file ) const;
00073 int restoreFactorization (const char * file , bool factor=false) ;
00075 void sort ( ) const;
00077 CoinFactorization & operator = ( const CoinFactorization & other );
00079
00089 int factorize ( const CoinPackedMatrix & matrix,
00090 int rowIsBasic[], int columnIsBasic[] ,
00091 double areaFactor = 0.0 );
00102 int factorize ( int numberRows,
00103 int numberColumns,
00104 CoinBigIndex numberElements,
00105 CoinBigIndex maximumL,
00106 CoinBigIndex maximumU,
00107 const int indicesRow[],
00108 const int indicesColumn[], const double elements[] ,
00109 int permutation[],
00110 double areaFactor = 0.0);
00115 int factorizePart1 ( int numberRows,
00116 int numberColumns,
00117 CoinBigIndex estimateNumberElements,
00118 int * indicesRow[],
00119 int * indicesColumn[],
00120 CoinFactorizationDouble * elements[],
00121 double areaFactor = 0.0);
00128 int factorizePart2 (int permutation[],int exactNumberElements);
00130 double conditionNumber() const;
00131
00133
00136
00137 inline int status ( ) const {
00138 return status_;
00139 }
00141 inline void setStatus ( int value)
00142 { status_=value; }
00144 inline int pivots ( ) const {
00145 return numberPivots_;
00146 }
00148 inline void setPivots ( int value )
00149 { numberPivots_=value; }
00151 inline int *permute ( ) const {
00152 return permute_.array();
00153 }
00155 inline int *pivotColumn ( ) const {
00156 return pivotColumn_.array();
00157 }
00159 inline CoinFactorizationDouble *pivotRegion ( ) const {
00160 return pivotRegion_.array();
00161 }
00163 inline int *permuteBack ( ) const {
00164 return permuteBack_.array();
00165 }
00167 inline int *lastRow ( ) const {
00168 return lastRow_.array();
00169 }
00172 inline int *pivotColumnBack ( ) const {
00173
00174 return pivotColumnBack_.array();
00175 }
00177 inline CoinBigIndex * startRowL() const
00178 { return startRowL_.array();}
00179
00181 inline CoinBigIndex * startColumnL() const
00182 { return startColumnL_.array();}
00183
00185 inline int * indexColumnL() const
00186 { return indexColumnL_.array();}
00187
00189 inline int * indexRowL() const
00190 { return indexRowL_.array();}
00191
00193 inline CoinFactorizationDouble * elementByRowL() const
00194 { return elementByRowL_.array();}
00195
00197 inline int numberRowsExtra ( ) const {
00198 return numberRowsExtra_;
00199 }
00201 inline void setNumberRows(int value)
00202 { numberRows_ = value; }
00204 inline int numberRows ( ) const {
00205 return numberRows_;
00206 }
00208 inline CoinBigIndex numberL() const
00209 { return numberL_;}
00210
00212 inline CoinBigIndex baseL() const
00213 { return baseL_;}
00215 inline int maximumRowsExtra ( ) const {
00216 return maximumRowsExtra_;
00217 }
00219 inline int numberColumns ( ) const {
00220 return numberColumns_;
00221 }
00223 inline int numberElements ( ) const {
00224 return totalElements_;
00225 }
00227 inline int numberForrestTomlin ( ) const {
00228 return numberInColumn_.array()[numberColumnsExtra_];
00229 }
00231 inline int numberGoodColumns ( ) const {
00232 return numberGoodU_;
00233 }
00235 inline double areaFactor ( ) const {
00236 return areaFactor_;
00237 }
00238 inline void areaFactor ( double value ) {
00239 areaFactor_=value;
00240 }
00242 double adjustedAreaFactor() const;
00244 inline void relaxAccuracyCheck(double value)
00245 { relaxCheck_ = value;}
00246 inline double getAccuracyCheck() const
00247 { return relaxCheck_;}
00249 inline int messageLevel ( ) const {
00250 return messageLevel_ ;
00251 }
00252 void messageLevel ( int value );
00254 inline int maximumPivots ( ) const {
00255 return maximumPivots_ ;
00256 }
00257 void maximumPivots ( int value );
00258
00260 inline int denseThreshold() const
00261 { return denseThreshold_;}
00263 inline void setDenseThreshold(int value)
00264 { denseThreshold_ = value;}
00266 inline double pivotTolerance ( ) const {
00267 return pivotTolerance_ ;
00268 }
00269 void pivotTolerance ( double value );
00271 inline double zeroTolerance ( ) const {
00272 return zeroTolerance_ ;
00273 }
00274 void zeroTolerance ( double value );
00275 #ifndef COIN_FAST_CODE
00277 inline double slackValue ( ) const {
00278 return slackValue_ ;
00279 }
00280 void slackValue ( double value );
00281 #endif
00283 double maximumCoefficient() const;
00285 inline bool forrestTomlin() const
00286 { return doForrestTomlin_;}
00287 inline void setForrestTomlin(bool value)
00288 { doForrestTomlin_=value;}
00290 inline bool spaceForForrestTomlin() const
00291 {
00292 CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
00293 CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
00294 return (space>=0)&&doForrestTomlin_;
00295 }
00297
00300
00302 inline int numberDense() const
00303 { return numberDense_;}
00304
00306 inline CoinBigIndex numberElementsU ( ) const {
00307 return lengthU_;
00308 }
00310 inline void setNumberElementsU(CoinBigIndex value)
00311 { lengthU_ = value; }
00313 inline CoinBigIndex lengthAreaU ( ) const {
00314 return lengthAreaU_;
00315 }
00317 inline CoinBigIndex numberElementsL ( ) const {
00318 return lengthL_;
00319 }
00321 inline CoinBigIndex lengthAreaL ( ) const {
00322 return lengthAreaL_;
00323 }
00325 inline CoinBigIndex numberElementsR ( ) const {
00326 return lengthR_;
00327 }
00329 inline CoinBigIndex numberCompressions() const
00330 { return numberCompressions_;}
00332 inline int * numberInRow() const
00333 { return numberInRow_.array();}
00335 inline int * numberInColumn() const
00336 { return numberInColumn_.array();}
00338 inline CoinFactorizationDouble * elementU() const
00339 { return elementU_.array();}
00341 inline int * indexRowU() const
00342 { return indexRowU_.array();}
00344 inline CoinBigIndex * startColumnU() const
00345 { return startColumnU_.array();}
00347 inline int maximumColumnsExtra()
00348 { return maximumColumnsExtra_;}
00352 inline int biasLU() const
00353 { return biasLU_;}
00354 inline void setBiasLU(int value)
00355 { biasLU_=value;}
00361 inline int persistenceFlag() const
00362 { return persistenceFlag_;}
00363 void setPersistenceFlag(int value);
00365
00368
00376 int replaceColumn ( CoinIndexedVector * regionSparse,
00377 int pivotRow,
00378 double pivotCheck ,
00379 bool checkBeforeModifying=false,
00380 double acceptablePivot=1.0e-8);
00385 void replaceColumnU ( CoinIndexedVector * regionSparse,
00386 CoinBigIndex * deleted,
00387 int internalPivotRow);
00388 #ifdef ABC_USE_COIN_FACTORIZATION
00389
00391 CoinIndexedVector * fakeVector(CoinIndexedVector * vector,
00392 int already=0) const;
00393 void deleteFakeVector(CoinIndexedVector * vector,
00394 CoinIndexedVector * fakeVector) const;
00399 double checkReplacePart1 ( CoinIndexedVector * regionSparse,
00400 int pivotRow);
00405 double checkReplacePart1 ( CoinIndexedVector * regionSparse,
00406 CoinIndexedVector * partialUpdate,
00407 int pivotRow);
00410 int checkReplacePart2 ( int pivotRow,
00411 double btranAlpha,
00412 double ftranAlpha,
00413 double ftAlpha,
00414 double acceptablePivot = 1.0e-8);
00417 void replaceColumnPart3 ( CoinIndexedVector * regionSparse,
00418 int pivotRow,
00419 double alpha );
00422 void replaceColumnPart3 ( CoinIndexedVector * regionSparse,
00423 CoinIndexedVector * partialUpdate,
00424 int pivotRow,
00425 double alpha );
00433 int updateColumnFT ( CoinIndexedVector & regionSparse);
00434 int updateColumnFTPart1 ( CoinIndexedVector & regionSparse) ;
00435 void updateColumnFTPart2 ( CoinIndexedVector & regionSparse) ;
00439 void updateColumnFT ( CoinIndexedVector & regionSparseFT,
00440 CoinIndexedVector & partialUpdate,
00441 int which);
00443 int updateColumn ( CoinIndexedVector & regionSparse) const;
00448 int updateTwoColumnsFT ( CoinIndexedVector & regionSparseFT,
00449 CoinIndexedVector & regionSparseOther);
00451 int updateColumnTranspose ( CoinIndexedVector & regionSparse) const;
00453 void updateColumnCpu ( CoinIndexedVector & regionSparse,int whichCpu) const;
00455 void updateColumnTransposeCpu ( CoinIndexedVector & regionSparse,int whichCpu) const;
00457 void updateFullColumn ( CoinIndexedVector & regionSparse) const;
00459 void updateFullColumnTranspose ( CoinIndexedVector & regionSparse) const;
00461 void updateWeights ( CoinIndexedVector & regionSparse) const;
00463 inline bool wantsTableauColumn() const
00464 {return false;}
00466 inline double minimumPivotTolerance ( ) const {
00467 return pivotTolerance_ ;
00468 }
00469 inline void minimumPivotTolerance ( double value )
00470 { pivotTolerance(value);}
00472 inline void setParallelMode(int value)
00473 { parallelMode_=value;}
00475 inline void setSolveMode(int value)
00476 { parallelMode_ &= 3;parallelMode_ |= (value<<2);}
00478 inline int solveMode() const
00479 { return parallelMode_ >> 2;}
00481 void updatePartialUpdate(CoinIndexedVector & partialUpdate);
00483 void makeNonSingular(int * COIN_RESTRICT sequence);
00484 #endif
00485
00486
00496 int updateColumnFT ( CoinIndexedVector * regionSparse,
00497 CoinIndexedVector * regionSparse2);
00500 int updateColumn ( CoinIndexedVector * regionSparse,
00501 CoinIndexedVector * regionSparse2,
00502 bool noPermute=false) const;
00508 int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
00509 CoinIndexedVector * regionSparse2,
00510 CoinIndexedVector * regionSparse3,
00511 bool noPermuteRegion3=false) ;
00516 int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00517 CoinIndexedVector * regionSparse2) const;
00519 void goSparse();
00521 inline int sparseThreshold ( ) const
00522 { return sparseThreshold_;}
00524 void sparseThreshold ( int value );
00526
00527
00531
00532 inline void clearArrays()
00533 { gutsOfDestructor();}
00535
00538
00541 int add ( CoinBigIndex numberElements,
00542 int indicesRow[],
00543 int indicesColumn[], double elements[] );
00544
00547 int addColumn ( CoinBigIndex numberElements,
00548 int indicesRow[], double elements[] );
00549
00552 int addRow ( CoinBigIndex numberElements,
00553 int indicesColumn[], double elements[] );
00554
00556 int deleteColumn ( int Row );
00558 int deleteRow ( int Row );
00559
00563 int replaceRow ( int whichRow, int numberElements,
00564 const int indicesColumn[], const double elements[] );
00566 void emptyRows(int numberToEmpty, const int which[]);
00568
00569
00570 void checkSparse();
00572 #if 0 //def CLP_FACTORIZATION_INSTRUMENT
00573 inline bool collectStatistics() const
00574 { return collectStatistics_;}
00576 inline void setCollectStatistics(bool onOff) const
00577 { collectStatistics_ = onOff;}
00578 #else
00579 inline bool collectStatistics() const
00580 { return true;}
00582 inline void setCollectStatistics(bool onOff) const
00583 { }
00584 #endif
00586 void gutsOfDestructor(int type=1);
00588 void gutsOfInitialize(int type);
00589 void gutsOfCopy(const CoinFactorization &other);
00590
00592 void resetStatistics();
00593
00594
00596
00598
00599 void getAreas ( int numberRows,
00600 int numberColumns,
00601 CoinBigIndex maximumL,
00602 CoinBigIndex maximumU );
00603
00606 void preProcess ( int state,
00607 int possibleDuplicates = -1 );
00609 int factor ( );
00610 protected:
00613 int factorSparse ( );
00616 int factorSparseSmall ( );
00619 int factorSparseLarge ( );
00622 int factorDense ( );
00623
00625 bool pivotOneOtherRow ( int pivotRow,
00626 int pivotColumn );
00628 bool pivotRowSingleton ( int pivotRow,
00629 int pivotColumn );
00631 bool pivotColumnSingleton ( int pivotRow,
00632 int pivotColumn );
00633
00638 bool getColumnSpace ( int iColumn,
00639 int extraNeeded );
00640
00643 bool reorderU();
00647 bool getColumnSpaceIterateR ( int iColumn, double value,
00648 int iRow);
00654 CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00655 int iRow);
00659 bool getRowSpace ( int iRow, int extraNeeded );
00660
00664 bool getRowSpaceIterate ( int iRow,
00665 int extraNeeded );
00667 void checkConsistency ( );
00669 inline void addLink ( int index, int count ) {
00670 int *nextCount = nextCount_.array();
00671 int *firstCount = firstCount_.array();
00672 int *lastCount = lastCount_.array();
00673 int next = firstCount[count];
00674 lastCount[index] = -2 - count;
00675 if ( next < 0 ) {
00676
00677 firstCount[count] = index;
00678 nextCount[index] = -1;
00679 } else {
00680 firstCount[count] = index;
00681 nextCount[index] = next;
00682 lastCount[next] = index;
00683 }}
00685 inline void deleteLink ( int index ) {
00686 int *nextCount = nextCount_.array();
00687 int *firstCount = firstCount_.array();
00688 int *lastCount = lastCount_.array();
00689 int next = nextCount[index];
00690 int last = lastCount[index];
00691 if ( last >= 0 ) {
00692 nextCount[last] = next;
00693 } else {
00694 int count = -last - 2;
00695
00696 firstCount[count] = next;
00697 }
00698 if ( next >= 0 ) {
00699 lastCount[next] = last;
00700 }
00701 nextCount[index] = -2;
00702 lastCount[index] = -2;
00703 return;
00704 }
00706 void separateLinks(int count,bool rowsFirst);
00708 void cleanup ( );
00709
00711 void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00713 void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00715 void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00717 void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00718
00720 void updateColumnR ( CoinIndexedVector * region ) const;
00723 void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00724
00726 void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00727
00729 void updateColumnUSparse ( CoinIndexedVector * regionSparse,
00730 int * indexIn) const;
00732 void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
00733 int * indexIn) const;
00735 int updateColumnUDensish ( double * COIN_RESTRICT region,
00736 int * COIN_RESTRICT regionIndex) const;
00738 void updateTwoColumnsUDensish (
00739 int & numberNonZero1,
00740 double * COIN_RESTRICT region1,
00741 int * COIN_RESTRICT index1,
00742 int & numberNonZero2,
00743 double * COIN_RESTRICT region2,
00744 int * COIN_RESTRICT index2) const;
00746 void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
00748 void permuteBack ( CoinIndexedVector * regionSparse,
00749 CoinIndexedVector * outVector) const;
00750
00752 void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00755 void updateColumnTransposeU ( CoinIndexedVector * region,
00756 int smallestIndex) const;
00759 void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00760 int smallestIndex) const;
00763 void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00764 int smallestIndex) const;
00767 void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00770 void updateColumnTransposeUByColumn ( CoinIndexedVector * region,
00771 int smallestIndex) const;
00772
00774 void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00776 void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00778 void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00779
00781 void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00783 void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00785 void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00787 void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00789 void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00790 public:
00795 int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00796 int pivotRow, double alpha);
00797 protected:
00800 int checkPivot(double saveFromU, double oldPivot) const;
00801
00802 #ifdef INT_IS_8
00803 #define COINFACTORIZATION_BITS_PER_INT 64
00804 #define COINFACTORIZATION_SHIFT_PER_INT 6
00805 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00806 #else
00807 #define COINFACTORIZATION_BITS_PER_INT 32
00808 #define COINFACTORIZATION_SHIFT_PER_INT 5
00809 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00810 #endif
00811 template <class T> inline bool
00812 pivot ( int pivotRow,
00813 int pivotColumn,
00814 CoinBigIndex pivotRowPosition,
00815 CoinBigIndex pivotColumnPosition,
00816 CoinFactorizationDouble work[],
00817 unsigned int workArea2[],
00818 int increment2,
00819 T markRow[] ,
00820 int largeInteger)
00821 {
00822 int *indexColumnU = indexColumnU_.array();
00823 CoinBigIndex *startColumnU = startColumnU_.array();
00824 int *numberInColumn = numberInColumn_.array();
00825 CoinFactorizationDouble *elementU = elementU_.array();
00826 int *indexRowU = indexRowU_.array();
00827 CoinBigIndex *startRowU = startRowU_.array();
00828 int *numberInRow = numberInRow_.array();
00829 CoinFactorizationDouble *elementL = elementL_.array();
00830 int *indexRowL = indexRowL_.array();
00831 int *saveColumn = saveColumn_.array();
00832 int *nextRow = nextRow_.array();
00833 int *lastRow = lastRow_.array() ;
00834
00835
00836 int numberInPivotRow = numberInRow[pivotRow] - 1;
00837 CoinBigIndex startColumn = startColumnU[pivotColumn];
00838 int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00839 CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00840 int put = 0;
00841 CoinBigIndex startRow = startRowU[pivotRow];
00842 CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00843
00844 if ( pivotColumnPosition < 0 ) {
00845 for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00846 int iColumn = indexColumnU[pivotColumnPosition];
00847 if ( iColumn != pivotColumn ) {
00848 saveColumn[put++] = iColumn;
00849 } else {
00850 break;
00851 }
00852 }
00853 } else {
00854 for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00855 saveColumn[put++] = indexColumnU[i];
00856 }
00857 }
00858 assert (pivotColumnPosition<endRow);
00859 assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00860 pivotColumnPosition++;
00861 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00862 saveColumn[put++] = indexColumnU[pivotColumnPosition];
00863 }
00864
00865 int next = nextRow[pivotRow];
00866 int last = lastRow[pivotRow];
00867
00868 nextRow[last] = next;
00869 lastRow[next] = last;
00870 nextRow[pivotRow] = numberGoodU_;
00871 lastRow[pivotRow] = -2;
00872 numberInRow[pivotRow] = 0;
00873
00874 CoinBigIndex l = lengthL_;
00875
00876 if ( l + numberInPivotColumn > lengthAreaL_ ) {
00877
00878 if ((messageLevel_&4)!=0)
00879 printf("more memory needed in middle of invert\n");
00880 return false;
00881 }
00882
00883 CoinBigIndex lSave = l;
00884
00885 CoinBigIndex * startColumnL = startColumnL_.array();
00886 startColumnL[numberGoodL_] = l;
00887 numberGoodL_++;
00888 startColumnL[numberGoodL_] = l + numberInPivotColumn;
00889 lengthL_ += numberInPivotColumn;
00890 if ( pivotRowPosition < 0 ) {
00891 for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00892 int iRow = indexRowU[pivotRowPosition];
00893 if ( iRow != pivotRow ) {
00894 indexRowL[l] = iRow;
00895 elementL[l] = elementU[pivotRowPosition];
00896 markRow[iRow] = static_cast<T>(l - lSave);
00897 l++;
00898
00899 CoinBigIndex start = startRowU[iRow];
00900 CoinBigIndex end = start + numberInRow[iRow];
00901 CoinBigIndex where = start;
00902
00903 while ( indexColumnU[where] != pivotColumn ) {
00904 where++;
00905 }
00906 #if DEBUG_COIN
00907 if ( where >= end ) {
00908 abort ( );
00909 }
00910 #endif
00911 indexColumnU[where] = indexColumnU[end - 1];
00912 numberInRow[iRow]--;
00913 } else {
00914 break;
00915 }
00916 }
00917 } else {
00918 CoinBigIndex i;
00919
00920 for ( i = startColumn; i < pivotRowPosition; i++ ) {
00921 int iRow = indexRowU[i];
00922
00923 markRow[iRow] = static_cast<T>(l - lSave);
00924 indexRowL[l] = iRow;
00925 elementL[l] = elementU[i];
00926 l++;
00927
00928 CoinBigIndex start = startRowU[iRow];
00929 CoinBigIndex end = start + numberInRow[iRow];
00930 CoinBigIndex where = start;
00931
00932 while ( indexColumnU[where] != pivotColumn ) {
00933 where++;
00934 }
00935 #if DEBUG_COIN
00936 if ( where >= end ) {
00937 abort ( );
00938 }
00939 #endif
00940 indexColumnU[where] = indexColumnU[end - 1];
00941 numberInRow[iRow]--;
00942 assert (numberInRow[iRow]>=0);
00943 }
00944 }
00945 assert (pivotRowPosition<endColumn);
00946 assert (indexRowU[pivotRowPosition]==pivotRow);
00947 CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
00948 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
00949
00950 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00951 pivotRowPosition++;
00952 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00953 int iRow = indexRowU[pivotRowPosition];
00954
00955 markRow[iRow] = static_cast<T>(l - lSave);
00956 indexRowL[l] = iRow;
00957 elementL[l] = elementU[pivotRowPosition];
00958 l++;
00959
00960 CoinBigIndex start = startRowU[iRow];
00961 CoinBigIndex end = start + numberInRow[iRow];
00962 CoinBigIndex where = start;
00963
00964 while ( indexColumnU[where] != pivotColumn ) {
00965 where++;
00966 }
00967 #if DEBUG_COIN
00968 if ( where >= end ) {
00969 abort ( );
00970 }
00971 #endif
00972 indexColumnU[where] = indexColumnU[end - 1];
00973 numberInRow[iRow]--;
00974 assert (numberInRow[iRow]>=0);
00975 }
00976 markRow[pivotRow] = static_cast<T>(largeInteger);
00977
00978 numberInColumn[pivotColumn] = 0;
00979
00980 int *indexL = &indexRowL[lSave];
00981 CoinFactorizationDouble *multipliersL = &elementL[lSave];
00982
00983
00984 int j;
00985
00986 for ( j = 0; j < numberInPivotColumn; j++ ) {
00987 multipliersL[j] *= pivotMultiplier;
00988 }
00989
00990 CoinBigIndex iErase;
00991 for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00992 iErase++ ) {
00993 workArea2[iErase] = 0;
00994 }
00995 CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00996 unsigned int *temp2 = workArea2;
00997 int * nextColumn = nextColumn_.array();
00998
00999
01000 int jColumn;
01001 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01002 int iColumn = saveColumn[jColumn];
01003 CoinBigIndex startColumn = startColumnU[iColumn];
01004 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
01005 int iRow = indexRowU[startColumn];
01006 CoinFactorizationDouble value = elementU[startColumn];
01007 double largest;
01008 CoinBigIndex put = startColumn;
01009 CoinBigIndex positionLargest = -1;
01010 CoinFactorizationDouble thisPivotValue = 0.0;
01011
01012
01013 bool checkLargest;
01014 int mark = markRow[iRow];
01015
01016 if ( mark == largeInteger+1 ) {
01017 largest = fabs ( value );
01018 positionLargest = put;
01019 put++;
01020 checkLargest = false;
01021 } else {
01022
01023 largest = 0.0;
01024 checkLargest = true;
01025 if ( mark != largeInteger ) {
01026
01027 work[mark] = value;
01028 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01029 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01030
01031 temp2[word] = temp2[word] | ( 1 << bit );
01032 added--;
01033 } else {
01034 thisPivotValue = value;
01035 }
01036 }
01037 CoinBigIndex i;
01038 for ( i = startColumn + 1; i < endColumn; i++ ) {
01039 iRow = indexRowU[i];
01040 value = elementU[i];
01041 int mark = markRow[iRow];
01042
01043 if ( mark == largeInteger+1 ) {
01044
01045 indexRowU[put] = iRow;
01046 elementU[put] = value;
01047 if ( checkLargest ) {
01048 double absValue = fabs ( value );
01049
01050 if ( absValue > largest ) {
01051 largest = absValue;
01052 positionLargest = put;
01053 }
01054 }
01055 put++;
01056 } else if ( mark != largeInteger ) {
01057
01058 work[mark] = value;
01059 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01060 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01061
01062 temp2[word] = temp2[word] | ( 1 << bit );
01063 added--;
01064 } else {
01065 thisPivotValue = value;
01066 }
01067 }
01068
01069 elementU[put] = elementU[startColumn];
01070 indexRowU[put] = indexRowU[startColumn];
01071 if ( positionLargest == startColumn ) {
01072 positionLargest = put;
01073 }
01074 put++;
01075 elementU[startColumn] = thisPivotValue;
01076 indexRowU[startColumn] = pivotRow;
01077
01078 startColumn++;
01079 numberInColumn[iColumn] = put - startColumn;
01080 int * numberInColumnPlus = numberInColumnPlus_.array();
01081 numberInColumnPlus[iColumn]++;
01082 startColumnU[iColumn]++;
01083
01084 int next = nextColumn[iColumn];
01085 CoinBigIndex space;
01086
01087 space = startColumnU[next] - put - numberInColumnPlus[next];
01088
01089 if ( numberInPivotColumn > space ) {
01090
01091 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
01092 return false;
01093 }
01094
01095 if (positionLargest >= 0)
01096 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
01097 startColumn = startColumnU[iColumn];
01098 put = startColumn + numberInColumn[iColumn];
01099 }
01100 double tolerance = zeroTolerance_;
01101
01102 int *nextCount = nextCount_.array();
01103 for ( j = 0; j < numberInPivotColumn; j++ ) {
01104 value = work[j] - thisPivotValue * multipliersL[j];
01105 double absValue = fabs ( value );
01106
01107 if ( absValue > tolerance ) {
01108 work[j] = 0.0;
01109 assert (put<lengthAreaU_);
01110 elementU[put] = value;
01111 indexRowU[put] = indexL[j];
01112 if ( absValue > largest ) {
01113 largest = absValue;
01114 positionLargest = put;
01115 }
01116 put++;
01117 } else {
01118 work[j] = 0.0;
01119 added--;
01120 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01121 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01122
01123 if ( temp2[word] & ( 1 << bit ) ) {
01124
01125 iRow = indexL[j];
01126 CoinBigIndex start = startRowU[iRow];
01127 CoinBigIndex end = start + numberInRow[iRow];
01128 CoinBigIndex where = start;
01129
01130 while ( indexColumnU[where] != iColumn ) {
01131 where++;
01132 }
01133 #if DEBUG_COIN
01134 if ( where >= end ) {
01135 abort ( );
01136 }
01137 #endif
01138 indexColumnU[where] = indexColumnU[end - 1];
01139 numberInRow[iRow]--;
01140 } else {
01141
01142 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01143 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01144
01145 temp2[word] = temp2[word] | ( 1 << bit );
01146 }
01147 }
01148 }
01149 numberInColumn[iColumn] = put - startColumn;
01150
01151 if ( positionLargest >= 0 ) {
01152 value = elementU[positionLargest];
01153 iRow = indexRowU[positionLargest];
01154 elementU[positionLargest] = elementU[startColumn];
01155 indexRowU[positionLargest] = indexRowU[startColumn];
01156 elementU[startColumn] = value;
01157 indexRowU[startColumn] = iRow;
01158 }
01159
01160 if ( nextCount[iColumn + numberRows_] != -2 ) {
01161
01162 deleteLink ( iColumn + numberRows_ );
01163 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01164 }
01165 temp2 += increment2;
01166 }
01167
01168 unsigned int *putBase = workArea2;
01169 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01170 int i = 0;
01171
01172
01173 while ( bigLoops ) {
01174 bigLoops--;
01175 int bit;
01176 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01177 unsigned int *putThis = putBase;
01178 int iRow = indexL[i];
01179
01180
01181 int number = 0;
01182 int jColumn;
01183
01184 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01185 unsigned int test = *putThis;
01186
01187 putThis += increment2;
01188 test = 1 - ( ( test >> bit ) & 1 );
01189 number += test;
01190 }
01191 int next = nextRow[iRow];
01192 CoinBigIndex space;
01193
01194 space = startRowU[next] - startRowU[iRow];
01195 number += numberInRow[iRow];
01196 if ( space < number ) {
01197 if ( !getRowSpace ( iRow, number ) ) {
01198 return false;
01199 }
01200 }
01201
01202 putThis = putBase;
01203 next = nextRow[iRow];
01204 number = numberInRow[iRow];
01205 CoinBigIndex end = startRowU[iRow] + number;
01206 int saveIndex = indexColumnU[startRowU[next]];
01207
01208
01209 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01210 unsigned int test = *putThis;
01211
01212 putThis += increment2;
01213 test = 1 - ( ( test >> bit ) & 1 );
01214 indexColumnU[end] = saveColumn[jColumn];
01215 end += test;
01216 }
01217
01218 indexColumnU[startRowU[next]] = saveIndex;
01219 markRow[iRow] = static_cast<T>(largeInteger+1);
01220 number = end - startRowU[iRow];
01221 numberInRow[iRow] = number;
01222 deleteLink ( iRow );
01223 addLink ( iRow, number );
01224 }
01225 putBase++;
01226 }
01227 int bit;
01228
01229 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01230 unsigned int *putThis = putBase;
01231 int iRow = indexL[i];
01232
01233
01234 int number = 0;
01235 int jColumn;
01236
01237 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01238 unsigned int test = *putThis;
01239
01240 putThis += increment2;
01241 test = 1 - ( ( test >> bit ) & 1 );
01242 number += test;
01243 }
01244 int next = nextRow[iRow];
01245 CoinBigIndex space;
01246
01247 space = startRowU[next] - startRowU[iRow];
01248 number += numberInRow[iRow];
01249 if ( space < number ) {
01250 if ( !getRowSpace ( iRow, number ) ) {
01251 return false;
01252 }
01253 }
01254
01255 putThis = putBase;
01256 next = nextRow[iRow];
01257 number = numberInRow[iRow];
01258 CoinBigIndex end = startRowU[iRow] + number;
01259 int saveIndex;
01260
01261 saveIndex = indexColumnU[startRowU[next]];
01262
01263
01264 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01265 unsigned int test = *putThis;
01266
01267 putThis += increment2;
01268 test = 1 - ( ( test >> bit ) & 1 );
01269
01270 indexColumnU[end] = saveColumn[jColumn];
01271 end += test;
01272 }
01273 indexColumnU[startRowU[next]] = saveIndex;
01274 markRow[iRow] = static_cast<T>(largeInteger+1);
01275 number = end - startRowU[iRow];
01276 numberInRow[iRow] = number;
01277 deleteLink ( iRow );
01278 addLink ( iRow, number );
01279 }
01280 markRow[pivotRow] = static_cast<T>(largeInteger+1);
01281
01282 deleteLink ( pivotRow );
01283 deleteLink ( pivotColumn + numberRows_ );
01284 totalElements_ += added;
01285 return true;
01286 }
01287
01288
01290
01291 protected:
01292
01295
01296 double pivotTolerance_;
01298 double zeroTolerance_;
01299 #ifndef COIN_FAST_CODE
01301 double slackValue_;
01302 #else
01303 #ifndef slackValue_
01304 #define slackValue_ -1.0
01305 #endif
01306 #endif
01308 double areaFactor_;
01310 double relaxCheck_;
01312 int numberRows_;
01314 int numberRowsExtra_;
01316 int maximumRowsExtra_;
01318 int numberColumns_;
01320 int numberColumnsExtra_;
01322 int maximumColumnsExtra_;
01324 int numberGoodU_;
01326 int numberGoodL_;
01328 int maximumPivots_;
01330 int numberPivots_;
01333 CoinBigIndex totalElements_;
01335 CoinBigIndex factorElements_;
01337 CoinIntArrayWithLength pivotColumn_;
01339 CoinIntArrayWithLength permute_;
01341 CoinIntArrayWithLength permuteBack_;
01343 CoinIntArrayWithLength pivotColumnBack_;
01345 int status_;
01346
01351
01352
01354 int numberTrials_;
01356 CoinBigIndexArrayWithLength startRowU_;
01357
01359 CoinIntArrayWithLength numberInRow_;
01360
01362 CoinIntArrayWithLength numberInColumn_;
01363
01365 CoinIntArrayWithLength numberInColumnPlus_;
01366
01369 CoinIntArrayWithLength firstCount_;
01370
01372 CoinIntArrayWithLength nextCount_;
01373
01375 CoinIntArrayWithLength lastCount_;
01376
01378 CoinIntArrayWithLength nextColumn_;
01379
01381 CoinIntArrayWithLength lastColumn_;
01382
01384 CoinIntArrayWithLength nextRow_;
01385
01387 CoinIntArrayWithLength lastRow_;
01388
01390 CoinIntArrayWithLength saveColumn_;
01391
01393 CoinIntArrayWithLength markRow_;
01394
01396 int messageLevel_;
01397
01399 int biggerDimension_;
01400
01402 CoinIntArrayWithLength indexColumnU_;
01403
01405 CoinIntArrayWithLength pivotRowL_;
01406
01408 CoinFactorizationDoubleArrayWithLength pivotRegion_;
01409
01411 int numberSlacks_;
01412
01414 int numberU_;
01415
01417 CoinBigIndex maximumU_;
01418
01420
01421
01423 CoinBigIndex lengthU_;
01424
01426 CoinBigIndex lengthAreaU_;
01427
01429 CoinFactorizationDoubleArrayWithLength elementU_;
01430
01432 CoinIntArrayWithLength indexRowU_;
01433
01435 CoinBigIndexArrayWithLength startColumnU_;
01436
01438 CoinBigIndexArrayWithLength convertRowToColumnU_;
01439
01441 CoinBigIndex numberL_;
01442
01444 CoinBigIndex baseL_;
01445
01447 CoinBigIndex lengthL_;
01448
01450 CoinBigIndex lengthAreaL_;
01451
01453 CoinFactorizationDoubleArrayWithLength elementL_;
01454
01456 CoinIntArrayWithLength indexRowL_;
01457
01459 CoinBigIndexArrayWithLength startColumnL_;
01460
01462 bool doForrestTomlin_;
01463
01465 int numberR_;
01466
01468 CoinBigIndex lengthR_;
01469
01471 CoinBigIndex lengthAreaR_;
01472
01474 CoinFactorizationDouble *elementR_;
01475
01477 int *indexRowR_;
01478
01480 CoinBigIndexArrayWithLength startColumnR_;
01481
01483 double * denseArea_;
01484
01486 double * denseAreaAddress_;
01487
01489 int * densePermute_;
01490
01492 int numberDense_;
01493
01495 int denseThreshold_;
01496
01498 CoinFactorizationDoubleArrayWithLength workArea_;
01499
01501 CoinUnsignedIntArrayWithLength workArea2_;
01502
01504 CoinBigIndex numberCompressions_;
01505
01506 public:
01508 mutable double ftranCountInput_;
01509 mutable double ftranCountAfterL_;
01510 mutable double ftranCountAfterR_;
01511 mutable double ftranCountAfterU_;
01512 mutable double btranCountInput_;
01513 mutable double btranCountAfterU_;
01514 mutable double btranCountAfterR_;
01515 mutable double btranCountAfterL_;
01516
01518 mutable int numberFtranCounts_;
01519 mutable int numberBtranCounts_;
01520
01522 double ftranAverageAfterL_;
01523 double ftranAverageAfterR_;
01524 double ftranAverageAfterU_;
01525 double btranAverageAfterU_;
01526 double btranAverageAfterR_;
01527 double btranAverageAfterL_;
01528 protected:
01529
01531 #if 0
01532 mutable bool collectStatistics_;
01533 #else
01534 #define collectStatistics_ 1
01535 #endif
01536
01538 int sparseThreshold_;
01539
01541 int sparseThreshold2_;
01542
01544 CoinBigIndexArrayWithLength startRowL_;
01545
01547 CoinIntArrayWithLength indexColumnL_;
01548
01550 CoinFactorizationDoubleArrayWithLength elementByRowL_;
01551
01553 mutable CoinIntArrayWithLength sparse_;
01557 int biasLU_;
01563 int persistenceFlag_;
01564 #ifdef ABC_USE_COIN_FACTORIZATION
01566 int parallelMode_;
01567 #endif
01568
01569 };
01570
01571 #ifdef COIN_HAS_LAPACK
01572 #ifndef COIN_FACTORIZATION_DENSE_CODE
01573 #define COIN_FACTORIZATION_DENSE_CODE 1
01574 #endif
01575 #endif
01576 #ifdef COIN_FACTORIZATION_DENSE_CODE
01577
01578 #ifndef ipfint
01579
01580 typedef int ipfint;
01581 typedef const int cipfint;
01582 #endif
01583 #endif
01584 #endif
01585
01586 #ifdef UGLY_COIN_FACTOR_CODING
01587 #define FAC_UNSET (FAC_SET+1)
01588 {
01589 goodPivot=false;
01590
01591 CoinBigIndex startColumnThis = startColumn[iPivotColumn];
01592 CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
01593 int put = 0;
01594 CoinBigIndex startRowThis = startRow[iPivotRow];
01595 CoinBigIndex endRow = startRowThis + numberDoRow + 1;
01596 if ( pivotColumnPosition < 0 ) {
01597 for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01598 int iColumn = indexColumn[pivotColumnPosition];
01599 if ( iColumn != iPivotColumn ) {
01600 saveColumn[put++] = iColumn;
01601 } else {
01602 break;
01603 }
01604 }
01605 } else {
01606 for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
01607 saveColumn[put++] = indexColumn[i];
01608 }
01609 }
01610 assert (pivotColumnPosition<endRow);
01611 assert (indexColumn[pivotColumnPosition]==iPivotColumn);
01612 pivotColumnPosition++;
01613 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01614 saveColumn[put++] = indexColumn[pivotColumnPosition];
01615 }
01616
01617 int next = nextRow[iPivotRow];
01618 int last = lastRow[iPivotRow];
01619
01620 nextRow[last] = next;
01621 lastRow[next] = last;
01622 nextRow[iPivotRow] = numberGoodU_;
01623 lastRow[iPivotRow] = -2;
01624 numberInRow[iPivotRow] = 0;
01625
01626 CoinBigIndex l = lengthL_;
01627
01628 {
01629 if ( l + numberDoColumn > lengthAreaL_ ) {
01630
01631 if ((messageLevel_&4)!=0)
01632 printf("more memory needed in middle of invert\n");
01633 goto BAD_PIVOT;
01634 }
01635
01636 CoinBigIndex lSave = l;
01637
01638 CoinBigIndex * startColumnL = startColumnL_.array();
01639 startColumnL[numberGoodL_] = l;
01640 numberGoodL_++;
01641 startColumnL[numberGoodL_] = l + numberDoColumn;
01642 lengthL_ += numberDoColumn;
01643 if ( pivotRowPosition < 0 ) {
01644 for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01645 int iRow = indexRow[pivotRowPosition];
01646 if ( iRow != iPivotRow ) {
01647 indexRowL[l] = iRow;
01648 elementL[l] = element[pivotRowPosition];
01649 markRow[iRow] = l - lSave;
01650 l++;
01651
01652 CoinBigIndex start = startRow[iRow];
01653 CoinBigIndex end = start + numberInRow[iRow];
01654 CoinBigIndex where = start;
01655
01656 while ( indexColumn[where] != iPivotColumn ) {
01657 where++;
01658 }
01659 #if DEBUG_COIN
01660 if ( where >= end ) {
01661 abort ( );
01662 }
01663 #endif
01664 indexColumn[where] = indexColumn[end - 1];
01665 numberInRow[iRow]--;
01666 } else {
01667 break;
01668 }
01669 }
01670 } else {
01671 CoinBigIndex i;
01672
01673 for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
01674 int iRow = indexRow[i];
01675
01676 markRow[iRow] = l - lSave;
01677 indexRowL[l] = iRow;
01678 elementL[l] = element[i];
01679 l++;
01680
01681 CoinBigIndex start = startRow[iRow];
01682 CoinBigIndex end = start + numberInRow[iRow];
01683 CoinBigIndex where = start;
01684
01685 while ( indexColumn[where] != iPivotColumn ) {
01686 where++;
01687 }
01688 #if DEBUG_COIN
01689 if ( where >= end ) {
01690 abort ( );
01691 }
01692 #endif
01693 indexColumn[where] = indexColumn[end - 1];
01694 numberInRow[iRow]--;
01695 assert (numberInRow[iRow]>=0);
01696 }
01697 }
01698 assert (pivotRowPosition<endColumn);
01699 assert (indexRow[pivotRowPosition]==iPivotRow);
01700 CoinFactorizationDouble pivotElement = element[pivotRowPosition];
01701 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
01702
01703 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
01704 pivotRowPosition++;
01705 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01706 int iRow = indexRow[pivotRowPosition];
01707
01708 markRow[iRow] = l - lSave;
01709 indexRowL[l] = iRow;
01710 elementL[l] = element[pivotRowPosition];
01711 l++;
01712
01713 CoinBigIndex start = startRow[iRow];
01714 CoinBigIndex end = start + numberInRow[iRow];
01715 CoinBigIndex where = start;
01716
01717 while ( indexColumn[where] != iPivotColumn ) {
01718 where++;
01719 }
01720 #if DEBUG_COIN
01721 if ( where >= end ) {
01722 abort ( );
01723 }
01724 #endif
01725 indexColumn[where] = indexColumn[end - 1];
01726 numberInRow[iRow]--;
01727 assert (numberInRow[iRow]>=0);
01728 }
01729 markRow[iPivotRow] = FAC_SET;
01730
01731 numberInColumn[iPivotColumn] = 0;
01732
01733 int *indexL = &indexRowL[lSave];
01734 CoinFactorizationDouble *multipliersL = &elementL[lSave];
01735
01736
01737 int j;
01738
01739 for ( j = 0; j < numberDoColumn; j++ ) {
01740 multipliersL[j] *= pivotMultiplier;
01741 }
01742
01743 CoinBigIndex iErase;
01744 for ( iErase = 0; iErase < increment2 * numberDoRow;
01745 iErase++ ) {
01746 workArea2[iErase] = 0;
01747 }
01748 CoinBigIndex added = numberDoRow * numberDoColumn;
01749 unsigned int *temp2 = workArea2;
01750 int * nextColumn = nextColumn_.array();
01751
01752
01753 int jColumn;
01754 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01755 int iColumn = saveColumn[jColumn];
01756 CoinBigIndex startColumnThis = startColumn[iColumn];
01757 CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
01758 int iRow = indexRow[startColumnThis];
01759 CoinFactorizationDouble value = element[startColumnThis];
01760 double largest;
01761 CoinBigIndex put = startColumnThis;
01762 CoinBigIndex positionLargest = -1;
01763 CoinFactorizationDouble thisPivotValue = 0.0;
01764
01765
01766 bool checkLargest;
01767 int mark = markRow[iRow];
01768
01769 if ( mark == FAC_UNSET ) {
01770 largest = fabs ( value );
01771 positionLargest = put;
01772 put++;
01773 checkLargest = false;
01774 } else {
01775
01776 largest = 0.0;
01777 checkLargest = true;
01778 if ( mark != FAC_SET ) {
01779
01780 workArea[mark] = value;
01781 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01782 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01783
01784 temp2[word] = temp2[word] | ( 1 << bit );
01785 added--;
01786 } else {
01787 thisPivotValue = value;
01788 }
01789 }
01790 CoinBigIndex i;
01791 for ( i = startColumnThis + 1; i < endColumn; i++ ) {
01792 iRow = indexRow[i];
01793 value = element[i];
01794 int mark = markRow[iRow];
01795
01796 if ( mark == FAC_UNSET ) {
01797
01798 indexRow[put] = iRow;
01799 element[put] = value;
01800 if ( checkLargest ) {
01801 double absValue = fabs ( value );
01802
01803 if ( absValue > largest ) {
01804 largest = absValue;
01805 positionLargest = put;
01806 }
01807 }
01808 put++;
01809 } else if ( mark != FAC_SET ) {
01810
01811 workArea[mark] = value;
01812 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01813 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01814
01815 temp2[word] = temp2[word] | ( 1 << bit );
01816 added--;
01817 } else {
01818 thisPivotValue = value;
01819 }
01820 }
01821
01822 element[put] = element[startColumnThis];
01823 indexRow[put] = indexRow[startColumnThis];
01824 if ( positionLargest == startColumnThis ) {
01825 positionLargest = put;
01826 }
01827 put++;
01828 element[startColumnThis] = thisPivotValue;
01829 indexRow[startColumnThis] = iPivotRow;
01830
01831 startColumnThis++;
01832 numberInColumn[iColumn] = put - startColumnThis;
01833 int * numberInColumnPlus = numberInColumnPlus_.array();
01834 numberInColumnPlus[iColumn]++;
01835 startColumn[iColumn]++;
01836
01837 int next = nextColumn[iColumn];
01838 CoinBigIndex space;
01839
01840 space = startColumn[next] - put - numberInColumnPlus[next];
01841
01842 if ( numberDoColumn > space ) {
01843
01844 if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
01845 goto BAD_PIVOT;
01846 }
01847
01848 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
01849 startColumnThis = startColumn[iColumn];
01850 put = startColumnThis + numberInColumn[iColumn];
01851 }
01852 double tolerance = zeroTolerance_;
01853
01854 int *nextCount = nextCount_.array();
01855 for ( j = 0; j < numberDoColumn; j++ ) {
01856 value = workArea[j] - thisPivotValue * multipliersL[j];
01857 double absValue = fabs ( value );
01858
01859 if ( absValue > tolerance ) {
01860 workArea[j] = 0.0;
01861 element[put] = value;
01862 indexRow[put] = indexL[j];
01863 if ( absValue > largest ) {
01864 largest = absValue;
01865 positionLargest = put;
01866 }
01867 put++;
01868 } else {
01869 workArea[j] = 0.0;
01870 added--;
01871 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01872 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01873
01874 if ( temp2[word] & ( 1 << bit ) ) {
01875
01876 iRow = indexL[j];
01877 CoinBigIndex start = startRow[iRow];
01878 CoinBigIndex end = start + numberInRow[iRow];
01879 CoinBigIndex where = start;
01880
01881 while ( indexColumn[where] != iColumn ) {
01882 where++;
01883 }
01884 #if DEBUG_COIN
01885 if ( where >= end ) {
01886 abort ( );
01887 }
01888 #endif
01889 indexColumn[where] = indexColumn[end - 1];
01890 numberInRow[iRow]--;
01891 } else {
01892
01893 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01894 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01895
01896 temp2[word] = temp2[word] | ( 1 << bit );
01897 }
01898 }
01899 }
01900 numberInColumn[iColumn] = put - startColumnThis;
01901
01902 if ( positionLargest >= 0 ) {
01903 value = element[positionLargest];
01904 iRow = indexRow[positionLargest];
01905 element[positionLargest] = element[startColumnThis];
01906 indexRow[positionLargest] = indexRow[startColumnThis];
01907 element[startColumnThis] = value;
01908 indexRow[startColumnThis] = iRow;
01909 }
01910
01911 if ( nextCount[iColumn + numberRows_] != -2 ) {
01912
01913 deleteLink ( iColumn + numberRows_ );
01914 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01915 }
01916 temp2 += increment2;
01917 }
01918
01919 unsigned int *putBase = workArea2;
01920 int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01921 int i = 0;
01922
01923
01924 while ( bigLoops ) {
01925 bigLoops--;
01926 int bit;
01927 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01928 unsigned int *putThis = putBase;
01929 int iRow = indexL[i];
01930
01931
01932 int number = 0;
01933 int jColumn;
01934
01935 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01936 unsigned int test = *putThis;
01937
01938 putThis += increment2;
01939 test = 1 - ( ( test >> bit ) & 1 );
01940 number += test;
01941 }
01942 int next = nextRow[iRow];
01943 CoinBigIndex space;
01944
01945 space = startRow[next] - startRow[iRow];
01946 number += numberInRow[iRow];
01947 if ( space < number ) {
01948 if ( !getRowSpace ( iRow, number ) ) {
01949 goto BAD_PIVOT;
01950 }
01951 }
01952
01953 putThis = putBase;
01954 next = nextRow[iRow];
01955 number = numberInRow[iRow];
01956 CoinBigIndex end = startRow[iRow] + number;
01957 int saveIndex = indexColumn[startRow[next]];
01958
01959
01960 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01961 unsigned int test = *putThis;
01962
01963 putThis += increment2;
01964 test = 1 - ( ( test >> bit ) & 1 );
01965 indexColumn[end] = saveColumn[jColumn];
01966 end += test;
01967 }
01968
01969 indexColumn[startRow[next]] = saveIndex;
01970 markRow[iRow] = FAC_UNSET;
01971 number = end - startRow[iRow];
01972 numberInRow[iRow] = number;
01973 deleteLink ( iRow );
01974 addLink ( iRow, number );
01975 }
01976 putBase++;
01977 }
01978 int bit;
01979
01980 for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
01981 unsigned int *putThis = putBase;
01982 int iRow = indexL[i];
01983
01984
01985 int number = 0;
01986 int jColumn;
01987
01988 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01989 unsigned int test = *putThis;
01990
01991 putThis += increment2;
01992 test = 1 - ( ( test >> bit ) & 1 );
01993 number += test;
01994 }
01995 int next = nextRow[iRow];
01996 CoinBigIndex space;
01997
01998 space = startRow[next] - startRow[iRow];
01999 number += numberInRow[iRow];
02000 if ( space < number ) {
02001 if ( !getRowSpace ( iRow, number ) ) {
02002 goto BAD_PIVOT;
02003 }
02004 }
02005
02006 putThis = putBase;
02007 next = nextRow[iRow];
02008 number = numberInRow[iRow];
02009 CoinBigIndex end = startRow[iRow] + number;
02010 int saveIndex;
02011
02012 saveIndex = indexColumn[startRow[next]];
02013
02014
02015 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
02016 unsigned int test = *putThis;
02017
02018 putThis += increment2;
02019 test = 1 - ( ( test >> bit ) & 1 );
02020
02021 indexColumn[end] = saveColumn[jColumn];
02022 end += test;
02023 }
02024 indexColumn[startRow[next]] = saveIndex;
02025 markRow[iRow] = FAC_UNSET;
02026 number = end - startRow[iRow];
02027 numberInRow[iRow] = number;
02028 deleteLink ( iRow );
02029 addLink ( iRow, number );
02030 }
02031 markRow[iPivotRow] = FAC_UNSET;
02032
02033 deleteLink ( iPivotRow );
02034 deleteLink ( iPivotColumn + numberRows_ );
02035 totalElements_ += added;
02036 goodPivot= true;
02037
02038 }
02039 BAD_PIVOT:
02040
02041 ;
02042 }
02043 #undef FAC_UNSET
02044 #endif