00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #ifndef CoinFactorization_H
00011 #define CoinFactorization_H
00012
00013 #include <iostream>
00014 #include <string>
00015 #include <cassert>
00016 #include "CoinFinite.hpp"
00017 #include "CoinIndexedVector.hpp"
00018 class CoinPackedMatrix;
00045 class CoinFactorization {
00046 friend void CoinFactorizationUnitTest( const std::string & mpsDir );
00047
00048 public:
00049
00052
00053 CoinFactorization ( );
00055 CoinFactorization ( const CoinFactorization &other);
00056
00058 ~CoinFactorization ( );
00060 void almostDestructor();
00062 void show_self ( ) const;
00064 int saveFactorization (const char * file ) const;
00068 int restoreFactorization (const char * file , bool factor=false) ;
00070 void sort ( ) const;
00072 CoinFactorization & operator = ( const CoinFactorization & other );
00074
00084 int factorize ( const CoinPackedMatrix & matrix,
00085 int rowIsBasic[], int columnIsBasic[] ,
00086 double areaFactor = 0.0 );
00097 int factorize ( int numberRows,
00098 int numberColumns,
00099 CoinBigIndex numberElements,
00100 CoinBigIndex maximumL,
00101 CoinBigIndex maximumU,
00102 const int indicesRow[],
00103 const int indicesColumn[], const double elements[] ,
00104 int permutation[],
00105 double areaFactor = 0.0);
00110 int factorizePart1 ( int numberRows,
00111 int numberColumns,
00112 CoinBigIndex estimateNumberElements,
00113 int * indicesRow[],
00114 int * indicesColumn[],
00115 double * elements[],
00116 double areaFactor = 0.0);
00123 int factorizePart2 (int permutation[],int exactNumberElements);
00125 double conditionNumber() const;
00126
00128
00131
00132 inline int status ( ) const {
00133 return status_;
00134 }
00136 inline void setStatus ( int value)
00137 { status_=value; }
00139 inline int pivots ( ) const {
00140 return numberPivots_;
00141 }
00143 inline void setPivots ( int value )
00144 { numberPivots_=value; }
00146 inline int *permute ( ) const {
00147 return permute_.array();
00148 }
00150 inline int *pivotColumn ( ) const {
00151 return pivotColumn_.array();
00152 }
00154 inline double *pivotRegion ( ) const {
00155 return pivotRegion_.array();
00156 }
00158 inline int *permuteBack ( ) const {
00159 return permuteBack_.array();
00160 }
00162 inline int *pivotColumnBack ( ) const {
00163 return pivotColumnBack_.array();
00164 }
00166 inline CoinBigIndex * startRowL() const
00167 { return startRowL_.array();}
00168
00170 inline CoinBigIndex * startColumnL() const
00171 { return startColumnL_.array();}
00172
00174 inline int * indexColumnL() const
00175 { return indexColumnL_.array();}
00176
00178 inline int * indexRowL() const
00179 { return indexRowL_.array();}
00180
00182 inline double * elementByRowL() const
00183 { return elementByRowL_.array();}
00184
00186 inline int numberRowsExtra ( ) const {
00187 return numberRowsExtra_;
00188 }
00190 inline void setNumberRows(int value)
00191 { numberRows_ = value; }
00193 inline int numberRows ( ) const {
00194 return numberRows_;
00195 }
00197 inline CoinBigIndex numberL() const
00198 { return numberL_;}
00199
00201 inline CoinBigIndex baseL() const
00202 { return baseL_;}
00204 inline int maximumRowsExtra ( ) const {
00205 return maximumRowsExtra_;
00206 }
00208 inline int numberColumns ( ) const {
00209 return numberColumns_;
00210 }
00212 inline int numberElements ( ) const {
00213 return totalElements_;
00214 }
00216 inline int numberForrestTomlin ( ) const {
00217 return numberInColumn_.array()[numberColumnsExtra_];
00218 }
00220 inline int numberGoodColumns ( ) const {
00221 return numberGoodU_;
00222 }
00224 inline double areaFactor ( ) const {
00225 return areaFactor_;
00226 }
00227 inline void areaFactor ( double value ) {
00228 areaFactor_=value;
00229 }
00231 double adjustedAreaFactor() const;
00233 inline void relaxAccuracyCheck(double value)
00234 { relaxCheck_ = value;}
00235 inline double getAccuracyCheck() const
00236 { return relaxCheck_;}
00238 inline bool increasingRows ( ) const
00239 { return true; }
00245 inline void increasingRows ( int value ) {}
00247 inline int messageLevel ( ) const {
00248 return messageLevel_ ;
00249 }
00250 void messageLevel ( int value );
00252 inline int maximumPivots ( ) const {
00253 return maximumPivots_ ;
00254 }
00255 void maximumPivots ( int value );
00256
00258 inline int denseThreshold() const
00259 { return denseThreshold_;}
00261 inline void setDenseThreshold(int value)
00262 { denseThreshold_ = value;}
00264 inline double pivotTolerance ( ) const {
00265 return pivotTolerance_ ;
00266 }
00267 void pivotTolerance ( double value );
00269 inline double zeroTolerance ( ) const {
00270 return zeroTolerance_ ;
00271 }
00272 void zeroTolerance ( double value );
00274 inline double slackValue ( ) const {
00275 return slackValue_ ;
00276 }
00277 void slackValue ( double value );
00279 double maximumCoefficient() const;
00281 inline bool forrestTomlin() const
00282 { return doForrestTomlin_;}
00283 inline void setForrestTomlin(bool value)
00284 { doForrestTomlin_=value;}
00286
00289
00291 inline int numberDense() const
00292 { return numberDense_;}
00293
00295 inline CoinBigIndex numberElementsU ( ) const {
00296 return lengthU_;
00297 }
00299 inline void setNumberElementsU(CoinBigIndex value)
00300 { lengthU_ = value; }
00302 inline CoinBigIndex lengthAreaU ( ) const {
00303 return lengthAreaU_;
00304 }
00306 inline CoinBigIndex numberElementsL ( ) const {
00307 return lengthL_;
00308 }
00310 inline CoinBigIndex lengthAreaL ( ) const {
00311 return lengthAreaL_;
00312 }
00314 inline CoinBigIndex numberElementsR ( ) const {
00315 return lengthR_;
00316 }
00318 inline CoinBigIndex numberCompressions() const
00319 { return numberCompressions_;}
00321 inline int * numberInRow() const
00322 { return numberInRow_.array();}
00324 inline int * numberInColumn() const
00325 { return numberInColumn_.array();}
00327 inline double * elementU() const
00328 { return elementU_.array();}
00330 inline int * indexRowU() const
00331 { return indexRowU_.array();}
00333 inline CoinBigIndex * startColumnU() const
00334 { return startColumnU_.array();}
00336 inline int maximumColumnsExtra()
00337 { return maximumColumnsExtra_;}
00341 inline int biasLU() const
00342 { return biasLU_;}
00343 inline void setBiasLU(int value)
00344 { biasLU_=value;}
00350 inline int persistenceFlag() const
00351 { return persistenceFlag_;}
00352 void setPersistenceFlag(int value);
00354
00357
00365 int replaceColumn ( CoinIndexedVector * regionSparse,
00366 int pivotRow,
00367 double pivotCheck ,
00368 bool checkBeforeModifying=false);
00370
00380 int updateColumnFT ( CoinIndexedVector * regionSparse,
00381 CoinIndexedVector * regionSparse2);
00384 int updateColumn ( CoinIndexedVector * regionSparse,
00385 CoinIndexedVector * regionSparse2,
00386 bool noPermute=false) const;
00391 int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00392 CoinIndexedVector * regionSparse2) const;
00394 void goSparse();
00396 inline int sparseThreshold ( ) const
00397 { return sparseThreshold_;}
00399 void sparseThreshold ( int value );
00401
00402
00406
00407 inline void clearArrays()
00408 { gutsOfDestructor();}
00410
00413
00416 int add ( CoinBigIndex numberElements,
00417 int indicesRow[],
00418 int indicesColumn[], double elements[] );
00419
00422 int addColumn ( CoinBigIndex numberElements,
00423 int indicesRow[], double elements[] );
00424
00427 int addRow ( CoinBigIndex numberElements,
00428 int indicesColumn[], double elements[] );
00429
00431 int deleteColumn ( int Row );
00433 int deleteRow ( int Row );
00434
00438 int replaceRow ( int whichRow, int numberElements,
00439 const int indicesColumn[], const double elements[] );
00441 void emptyRows(int numberToEmpty, const int which[]);
00443
00444
00445 void checkSparse();
00447 inline bool collectStatistics() const
00448 { return collectStatistics_;}
00450 inline void setCollectStatistics(bool onOff) const
00451 { collectStatistics_ = onOff;}
00453 void gutsOfDestructor(int type=1);
00455 void gutsOfInitialize(int type);
00456 void gutsOfCopy(const CoinFactorization &other);
00457
00459 void resetStatistics();
00460
00461
00463
00465
00466 void getAreas ( int numberRows,
00467 int numberColumns,
00468 CoinBigIndex maximumL,
00469 CoinBigIndex maximumU );
00470
00473 void preProcess ( int state,
00474 int possibleDuplicates = -1 );
00476 int factor ( );
00477 protected:
00480 int factorSparse ( );
00483 int factorDense ( );
00484
00486 bool pivotOneOtherRow ( int pivotRow,
00487 int pivotColumn );
00489 bool pivotRowSingleton ( int pivotRow,
00490 int pivotColumn );
00492 bool pivotColumnSingleton ( int pivotRow,
00493 int pivotColumn );
00494
00499 bool getColumnSpace ( int iColumn,
00500 int extraNeeded );
00501
00505 bool getColumnSpaceIterateR ( int iColumn, double value,
00506 int iRow);
00512 CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00513 int iRow);
00517 bool getRowSpace ( int iRow, int extraNeeded );
00518
00522 bool getRowSpaceIterate ( int iRow,
00523 int extraNeeded );
00525 void checkConsistency ( );
00527 inline void addLink ( int index, int count ) {
00528 int *nextCount = nextCount_.array();
00529 int *firstCount = firstCount_.array();
00530 int *lastCount = lastCount_.array();
00531 int next = firstCount[count];
00532 lastCount[index] = -2 - count;
00533 if ( next < 0 ) {
00534
00535 firstCount[count] = index;
00536 nextCount[index] = -1;
00537 } else {
00538 firstCount[count] = index;
00539 nextCount[index] = next;
00540 lastCount[next] = index;
00541 }}
00543 inline void deleteLink ( int index ) {
00544 int *nextCount = nextCount_.array();
00545 int *firstCount = firstCount_.array();
00546 int *lastCount = lastCount_.array();
00547 int next = nextCount[index];
00548 int last = lastCount[index];
00549 if ( last >= 0 ) {
00550 nextCount[last] = next;
00551 } else {
00552 int count = -last - 2;
00553
00554 firstCount[count] = next;
00555 }
00556 if ( next >= 0 ) {
00557 lastCount[next] = last;
00558 }
00559 nextCount[index] = -2;
00560 lastCount[index] = -2;
00561 return;
00562 }
00564 void separateLinks(int count,bool rowsFirst);
00566 void cleanup ( );
00567
00569 void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00571 void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00573 void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00575 void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00576
00578 void updateColumnR ( CoinIndexedVector * region ) const;
00581 void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00582
00584 void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00585
00587 void updateColumnUSparse ( CoinIndexedVector * regionSparse,
00588 int * indexIn) const;
00590 void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
00591 int * indexIn) const;
00593 void updateColumnUDensish ( CoinIndexedVector * regionSparse,
00594 int * indexIn) const;
00596 void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
00598 void permuteBack ( CoinIndexedVector * regionSparse,
00599 CoinIndexedVector * outVector) const;
00600
00602 void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00605 void updateColumnTransposeU ( CoinIndexedVector * region,
00606 int smallestIndex) const;
00609 void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00610 int smallestIndex) const;
00613 void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00614 int smallestIndex) const;
00617 void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00618
00620 void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00622 void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00624 void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00625
00627 void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00629 void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00631 void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00633 void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00635 void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00636 public:
00641 int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00642 int pivotRow, double alpha);
00643 protected:
00646 int checkPivot(double saveFromU, double oldPivot) const;
00647
00648 #ifdef INT_IS_8
00649 #define COINFACTORIZATION_BITS_PER_INT 64
00650 #define COINFACTORIZATION_SHIFT_PER_INT 6
00651 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00652 #else
00653 #define COINFACTORIZATION_BITS_PER_INT 32
00654 #define COINFACTORIZATION_SHIFT_PER_INT 5
00655 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00656 #endif
00657 template <class T> inline bool
00658 pivot ( int pivotRow,
00659 int pivotColumn,
00660 CoinBigIndex pivotRowPosition,
00661 CoinBigIndex pivotColumnPosition,
00662 double work[],
00663 unsigned int workArea2[],
00664 int increment,
00665 int increment2,
00666 T markRow[] ,
00667 int largeInteger)
00668 {
00669 int *indexColumnU = indexColumnU_.array();
00670 CoinBigIndex *startColumnU = startColumnU_.array();
00671 int *numberInColumn = numberInColumn_.array();
00672 double *elementU = elementU_.array();
00673 int *indexRowU = indexRowU_.array();
00674 CoinBigIndex *startRowU = startRowU_.array();
00675 int *numberInRow = numberInRow_.array();
00676 double *elementL = elementL_.array();
00677 int *indexRowL = indexRowL_.array();
00678 int *saveColumn = saveColumn_.array();
00679 int *nextRow = nextRow_.array();
00680 int *lastRow = lastRow_.array() ;
00681
00682
00683 int numberInPivotRow = numberInRow[pivotRow] - 1;
00684 CoinBigIndex startColumn = startColumnU[pivotColumn];
00685 int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00686 CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00687 int put = 0;
00688 CoinBigIndex startRow = startRowU[pivotRow];
00689 CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00690
00691 if ( pivotColumnPosition < 0 ) {
00692 for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00693 int iColumn = indexColumnU[pivotColumnPosition];
00694 if ( iColumn != pivotColumn ) {
00695 saveColumn[put++] = iColumn;
00696 } else {
00697 break;
00698 }
00699 }
00700 } else {
00701 for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00702 saveColumn[put++] = indexColumnU[i];
00703 }
00704 }
00705 assert (pivotColumnPosition<endRow);
00706 assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00707 pivotColumnPosition++;
00708 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00709 saveColumn[put++] = indexColumnU[pivotColumnPosition];
00710 }
00711
00712 int next = nextRow[pivotRow];
00713 int last = lastRow[pivotRow];
00714
00715 nextRow[last] = next;
00716 lastRow[next] = last;
00717 nextRow[pivotRow] = numberGoodU_;
00718 lastRow[pivotRow] = -2;
00719 numberInRow[pivotRow] = 0;
00720
00721 CoinBigIndex l = lengthL_;
00722
00723 if ( l + numberInPivotColumn > lengthAreaL_ ) {
00724
00725 printf("more memory needed in middle of invert\n");
00726 return false;
00727 }
00728
00729 CoinBigIndex lSave = l;
00730
00731 pivotRowL_.array()[numberGoodL_] = pivotRow;
00732 CoinBigIndex * startColumnL = startColumnL_.array();
00733 startColumnL[numberGoodL_] = l;
00734 numberGoodL_++;
00735 startColumnL[numberGoodL_] = l + numberInPivotColumn;
00736 lengthL_ += numberInPivotColumn;
00737 if ( pivotRowPosition < 0 ) {
00738 for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00739 int iRow = indexRowU[pivotRowPosition];
00740 if ( iRow != pivotRow ) {
00741 indexRowL[l] = iRow;
00742 elementL[l] = elementU[pivotRowPosition];
00743 markRow[iRow] = l - lSave;
00744 l++;
00745
00746 CoinBigIndex start = startRowU[iRow];
00747 CoinBigIndex end = start + numberInRow[iRow];
00748 CoinBigIndex where = start;
00749
00750 while ( indexColumnU[where] != pivotColumn ) {
00751 where++;
00752 }
00753 #if DEBUG_COIN
00754 if ( where >= end ) {
00755 abort ( );
00756 }
00757 #endif
00758 indexColumnU[where] = indexColumnU[end - 1];
00759 numberInRow[iRow]--;
00760 } else {
00761 break;
00762 }
00763 }
00764 } else {
00765 CoinBigIndex i;
00766
00767 for ( i = startColumn; i < pivotRowPosition; i++ ) {
00768 int iRow = indexRowU[i];
00769
00770 markRow[iRow] = l - lSave;
00771 indexRowL[l] = iRow;
00772 elementL[l] = elementU[i];
00773 l++;
00774
00775 CoinBigIndex start = startRowU[iRow];
00776 CoinBigIndex end = start + numberInRow[iRow];
00777 CoinBigIndex where = start;
00778
00779 while ( indexColumnU[where] != pivotColumn ) {
00780 where++;
00781 }
00782 #if DEBUG_COIN
00783 if ( where >= end ) {
00784 abort ( );
00785 }
00786 #endif
00787 indexColumnU[where] = indexColumnU[end - 1];
00788 numberInRow[iRow]--;
00789 assert (numberInRow[iRow]>=0);
00790 }
00791 }
00792 assert (pivotRowPosition<endColumn);
00793 assert (indexRowU[pivotRowPosition]==pivotRow);
00794 double pivotElement = elementU[pivotRowPosition];
00795 double pivotMultiplier = 1.0 / pivotElement;
00796
00797 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00798 pivotRowPosition++;
00799 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00800 int iRow = indexRowU[pivotRowPosition];
00801
00802 markRow[iRow] = l - lSave;
00803 indexRowL[l] = iRow;
00804 elementL[l] = elementU[pivotRowPosition];
00805 l++;
00806
00807 CoinBigIndex start = startRowU[iRow];
00808 CoinBigIndex end = start + numberInRow[iRow];
00809 CoinBigIndex where = start;
00810
00811 while ( indexColumnU[where] != pivotColumn ) {
00812 where++;
00813 }
00814 #if DEBUG_COIN
00815 if ( where >= end ) {
00816 abort ( );
00817 }
00818 #endif
00819 indexColumnU[where] = indexColumnU[end - 1];
00820 numberInRow[iRow]--;
00821 assert (numberInRow[iRow]>=0);
00822 }
00823 markRow[pivotRow] = largeInteger;
00824
00825 numberInColumn[pivotColumn] = 0;
00826
00827 int *indexL = &indexRowL[lSave];
00828 double *multipliersL = &elementL[lSave];
00829
00830
00831 int j;
00832
00833 for ( j = 0; j < numberInPivotColumn; j++ ) {
00834 multipliersL[j] *= pivotMultiplier;
00835 }
00836
00837 CoinBigIndex iErase;
00838 for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00839 iErase++ ) {
00840 workArea2[iErase] = 0;
00841 }
00842 CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00843 unsigned int *temp2 = workArea2;
00844 int * nextColumn = nextColumn_.array();
00845
00846
00847 int jColumn;
00848 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
00849 int iColumn = saveColumn[jColumn];
00850 CoinBigIndex startColumn = startColumnU[iColumn];
00851 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
00852 int iRow = indexRowU[startColumn];
00853 double value = elementU[startColumn];
00854 double largest;
00855 CoinBigIndex put = startColumn;
00856 CoinBigIndex positionLargest = -1;
00857 double thisPivotValue = 0.0;
00858
00859
00860 bool checkLargest;
00861 int mark = markRow[iRow];
00862
00863 if ( mark < 0 ) {
00864 largest = fabs ( value );
00865 positionLargest = put;
00866 put++;
00867 checkLargest = false;
00868 } else {
00869
00870 largest = 0.0;
00871 checkLargest = true;
00872 if ( mark != largeInteger ) {
00873
00874 work[mark] = value;
00875 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00876 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00877
00878 temp2[word] = temp2[word] | ( 1 << bit );
00879 added--;
00880 } else {
00881 thisPivotValue = value;
00882 }
00883 }
00884 CoinBigIndex i;
00885 for ( i = startColumn + 1; i < endColumn; i++ ) {
00886 iRow = indexRowU[i];
00887 value = elementU[i];
00888 int mark = markRow[iRow];
00889
00890 if ( mark < 0 ) {
00891
00892 indexRowU[put] = iRow;
00893 elementU[put] = value;
00894 if ( checkLargest ) {
00895 double absValue = fabs ( value );
00896
00897 if ( absValue > largest ) {
00898 largest = absValue;
00899 positionLargest = put;
00900 }
00901 }
00902 put++;
00903 } else if ( mark != largeInteger ) {
00904
00905 work[mark] = value;
00906 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00907 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00908
00909 temp2[word] = temp2[word] | ( 1 << bit );
00910 added--;
00911 } else {
00912 thisPivotValue = value;
00913 }
00914 }
00915
00916 elementU[put] = elementU[startColumn];
00917 indexRowU[put] = indexRowU[startColumn];
00918 if ( positionLargest == startColumn ) {
00919 positionLargest = put;
00920 }
00921 put++;
00922 elementU[startColumn] = thisPivotValue;
00923 indexRowU[startColumn] = pivotRow;
00924
00925 startColumn++;
00926 numberInColumn[iColumn] = put - startColumn;
00927 int * numberInColumnPlus = numberInColumnPlus_.array();
00928 numberInColumnPlus[iColumn]++;
00929 startColumnU[iColumn]++;
00930
00931 int next = nextColumn[iColumn];
00932 CoinBigIndex space;
00933
00934 space = startColumnU[next] - put - numberInColumnPlus[next];
00935
00936 if ( numberInPivotColumn > space ) {
00937
00938 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
00939 return false;
00940 }
00941
00942 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
00943 startColumn = startColumnU[iColumn];
00944 put = startColumn + numberInColumn[iColumn];
00945 }
00946 double tolerance = zeroTolerance_;
00947
00948 int *nextCount = nextCount_.array();
00949 for ( j = 0; j < numberInPivotColumn; j++ ) {
00950 value = work[j] - thisPivotValue * multipliersL[j];
00951 double absValue = fabs ( value );
00952
00953 if ( absValue > tolerance ) {
00954 work[j] = 0.0;
00955 elementU[put] = value;
00956 indexRowU[put] = indexL[j];
00957 if ( absValue > largest ) {
00958 largest = absValue;
00959 positionLargest = put;
00960 }
00961 put++;
00962 } else {
00963 work[j] = 0.0;
00964 added--;
00965 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
00966 int bit = j & COINFACTORIZATION_MASK_PER_INT;
00967
00968 if ( temp2[word] & ( 1 << bit ) ) {
00969
00970 iRow = indexL[j];
00971 CoinBigIndex start = startRowU[iRow];
00972 CoinBigIndex end = start + numberInRow[iRow];
00973 CoinBigIndex where = start;
00974
00975 while ( indexColumnU[where] != iColumn ) {
00976 where++;
00977 }
00978 #if DEBUG_COIN
00979 if ( where >= end ) {
00980 abort ( );
00981 }
00982 #endif
00983 indexColumnU[where] = indexColumnU[end - 1];
00984 numberInRow[iRow]--;
00985 } else {
00986
00987 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
00988 int bit = j & COINFACTORIZATION_MASK_PER_INT;
00989
00990 temp2[word] = temp2[word] | ( 1 << bit );
00991 }
00992 }
00993 }
00994 numberInColumn[iColumn] = put - startColumn;
00995
00996 if ( positionLargest >= 0 ) {
00997 value = elementU[positionLargest];
00998 iRow = indexRowU[positionLargest];
00999 elementU[positionLargest] = elementU[startColumn];
01000 indexRowU[positionLargest] = indexRowU[startColumn];
01001 elementU[startColumn] = value;
01002 indexRowU[startColumn] = iRow;
01003 }
01004
01005 if ( nextCount[iColumn + numberRows_] != -2 ) {
01006
01007 deleteLink ( iColumn + numberRows_ );
01008 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01009 }
01010 temp2 += increment2;
01011 }
01012
01013 unsigned int *putBase = workArea2;
01014 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01015 int i = 0;
01016
01017
01018 while ( bigLoops ) {
01019 bigLoops--;
01020 int bit;
01021 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01022 unsigned int *putThis = putBase;
01023 int iRow = indexL[i];
01024
01025
01026 int number = 0;
01027 int jColumn;
01028
01029 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01030 unsigned int test = *putThis;
01031
01032 putThis += increment2;
01033 test = 1 - ( ( test >> bit ) & 1 );
01034 number += test;
01035 }
01036 int next = nextRow[iRow];
01037 CoinBigIndex space;
01038
01039 space = startRowU[next] - startRowU[iRow];
01040 number += numberInRow[iRow];
01041 if ( space < number ) {
01042 if ( !getRowSpace ( iRow, number ) ) {
01043 return false;
01044 }
01045 }
01046
01047 putThis = putBase;
01048 next = nextRow[iRow];
01049 number = numberInRow[iRow];
01050 CoinBigIndex end = startRowU[iRow] + number;
01051 int saveIndex = indexColumnU[startRowU[next]];
01052
01053
01054 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01055 unsigned int test = *putThis;
01056
01057 putThis += increment2;
01058 test = 1 - ( ( test >> bit ) & 1 );
01059 indexColumnU[end] = saveColumn[jColumn];
01060 end += test;
01061 }
01062
01063 indexColumnU[startRowU[next]] = saveIndex;
01064 markRow[iRow] = -1;
01065 number = end - startRowU[iRow];
01066 numberInRow[iRow] = number;
01067 deleteLink ( iRow );
01068 addLink ( iRow, number );
01069 }
01070 putBase++;
01071 }
01072 int bit;
01073
01074 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01075 unsigned int *putThis = putBase;
01076 int iRow = indexL[i];
01077
01078
01079 int number = 0;
01080 int jColumn;
01081
01082 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01083 unsigned int test = *putThis;
01084
01085 putThis += increment2;
01086 test = 1 - ( ( test >> bit ) & 1 );
01087 number += test;
01088 }
01089 int next = nextRow[iRow];
01090 CoinBigIndex space;
01091
01092 space = startRowU[next] - startRowU[iRow];
01093 number += numberInRow[iRow];
01094 if ( space < number ) {
01095 if ( !getRowSpace ( iRow, number ) ) {
01096 return false;
01097 }
01098 }
01099
01100 putThis = putBase;
01101 next = nextRow[iRow];
01102 number = numberInRow[iRow];
01103 CoinBigIndex end = startRowU[iRow] + number;
01104 int saveIndex;
01105
01106 saveIndex = indexColumnU[startRowU[next]];
01107
01108
01109 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01110 unsigned int test = *putThis;
01111
01112 putThis += increment2;
01113 test = 1 - ( ( test >> bit ) & 1 );
01114
01115 indexColumnU[end] = saveColumn[jColumn];
01116 end += test;
01117 }
01118 indexColumnU[startRowU[next]] = saveIndex;
01119 markRow[iRow] = -1;
01120 number = end - startRowU[iRow];
01121 numberInRow[iRow] = number;
01122 deleteLink ( iRow );
01123 addLink ( iRow, number );
01124 }
01125 markRow[pivotRow] = -1;
01126
01127 deleteLink ( pivotRow );
01128 deleteLink ( pivotColumn + numberRows_ );
01129 totalElements_ += added;
01130 return true;
01131 }
01132
01133
01135
01136 protected:
01137
01140
01141 double pivotTolerance_;
01143 double zeroTolerance_;
01145 double slackValue_;
01147 double areaFactor_;
01149 double relaxCheck_;
01151 int numberRows_;
01153 int numberRowsExtra_;
01155 int maximumRowsExtra_;
01157 int numberColumns_;
01159 int numberColumnsExtra_;
01161 int maximumColumnsExtra_;
01163 int numberGoodU_;
01165 int numberGoodL_;
01167 int maximumPivots_;
01169 int numberPivots_;
01172 CoinBigIndex totalElements_;
01174 CoinBigIndex factorElements_;
01176 CoinIntArrayWithLength pivotColumn_;
01178 CoinIntArrayWithLength permute_;
01180 CoinIntArrayWithLength permuteBack_;
01182 CoinIntArrayWithLength pivotColumnBack_;
01184 int status_;
01185
01190
01191
01193 int numberTrials_;
01195 CoinBigIndexArrayWithLength startRowU_;
01196
01198 CoinIntArrayWithLength numberInRow_;
01199
01201 CoinIntArrayWithLength numberInColumn_;
01202
01204 CoinIntArrayWithLength numberInColumnPlus_;
01205
01208 CoinIntArrayWithLength firstCount_;
01209
01211 CoinIntArrayWithLength nextCount_;
01212
01214 CoinIntArrayWithLength lastCount_;
01215
01217 CoinIntArrayWithLength nextColumn_;
01218
01220 CoinIntArrayWithLength lastColumn_;
01221
01223 CoinIntArrayWithLength nextRow_;
01224
01226 CoinIntArrayWithLength lastRow_;
01227
01229 CoinIntArrayWithLength saveColumn_;
01230
01232 CoinIntArrayWithLength markRow_;
01233
01235 int messageLevel_;
01236
01238 int biggerDimension_;
01239
01241 CoinIntArrayWithLength indexColumnU_;
01242
01244 CoinIntArrayWithLength pivotRowL_;
01245
01247 CoinDoubleArrayWithLength pivotRegion_;
01248
01250 int numberSlacks_;
01251
01253 int numberU_;
01254
01256 CoinBigIndex maximumU_;
01257
01259
01260
01262 CoinBigIndex lengthU_;
01263
01265 CoinBigIndex lengthAreaU_;
01266
01268 CoinDoubleArrayWithLength elementU_;
01269
01271 CoinIntArrayWithLength indexRowU_;
01272
01274 CoinBigIndexArrayWithLength startColumnU_;
01275
01277 CoinBigIndexArrayWithLength convertRowToColumnU_;
01278
01280 CoinBigIndex numberL_;
01281
01283 CoinBigIndex baseL_;
01284
01286 CoinBigIndex lengthL_;
01287
01289 CoinBigIndex lengthAreaL_;
01290
01292 CoinDoubleArrayWithLength elementL_;
01293
01295 CoinIntArrayWithLength indexRowL_;
01296
01298 CoinBigIndexArrayWithLength startColumnL_;
01299
01301 bool doForrestTomlin_;
01302
01304 int numberR_;
01305
01307 CoinBigIndex lengthR_;
01308
01310 CoinBigIndex lengthAreaR_;
01311
01313 double *elementR_;
01314
01316 int *indexRowR_;
01317
01319 CoinBigIndexArrayWithLength startColumnR_;
01320
01322 double * denseArea_;
01323
01325 int * densePermute_;
01326
01328 int numberDense_;
01329
01331 int denseThreshold_;
01332
01334 CoinDoubleArrayWithLength workArea_;
01335
01337 CoinUnsignedIntArrayWithLength workArea2_;
01338
01340 CoinBigIndex numberCompressions_;
01341
01343 mutable double ftranCountInput_;
01344 mutable double ftranCountAfterL_;
01345 mutable double ftranCountAfterR_;
01346 mutable double ftranCountAfterU_;
01347 mutable double btranCountInput_;
01348 mutable double btranCountAfterU_;
01349 mutable double btranCountAfterR_;
01350 mutable double btranCountAfterL_;
01351
01353 mutable int numberFtranCounts_;
01354 mutable int numberBtranCounts_;
01355
01357 double ftranAverageAfterL_;
01358 double ftranAverageAfterR_;
01359 double ftranAverageAfterU_;
01360 double btranAverageAfterU_;
01361 double btranAverageAfterR_;
01362 double btranAverageAfterL_;
01363
01365 mutable bool collectStatistics_;
01366
01368 int sparseThreshold_;
01369
01371 int sparseThreshold2_;
01372
01374 CoinBigIndexArrayWithLength startRowL_;
01375
01377 CoinIntArrayWithLength indexColumnL_;
01378
01380 CoinDoubleArrayWithLength elementByRowL_;
01381
01383 mutable CoinIntArrayWithLength sparse_;
01387 int biasLU_;
01393 int persistenceFlag_;
01395 };
01396
01397 #ifdef COIN_HAS_LAPACK
01398 #define DENSE_CODE 1
01399
01400 #ifndef ipfint
01401
01402 typedef int ipfint;
01403 typedef const int cipfint;
01404 #endif
01405 #endif
01406 #endif