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 int messageLevel ( ) const {
00239 return messageLevel_ ;
00240 }
00241 void messageLevel ( int value );
00243 inline int maximumPivots ( ) const {
00244 return maximumPivots_ ;
00245 }
00246 void maximumPivots ( int value );
00247
00249 inline int denseThreshold() const
00250 { return denseThreshold_;}
00252 inline void setDenseThreshold(int value)
00253 { denseThreshold_ = value;}
00255 inline double pivotTolerance ( ) const {
00256 return pivotTolerance_ ;
00257 }
00258 void pivotTolerance ( double value );
00260 inline double zeroTolerance ( ) const {
00261 return zeroTolerance_ ;
00262 }
00263 void zeroTolerance ( double value );
00264 #ifndef COIN_FAST_CODE
00266 inline double slackValue ( ) const {
00267 return slackValue_ ;
00268 }
00269 void slackValue ( double value );
00270 #endif
00272 double maximumCoefficient() const;
00274 inline bool forrestTomlin() const
00275 { return doForrestTomlin_;}
00276 inline void setForrestTomlin(bool value)
00277 { doForrestTomlin_=value;}
00279 inline bool spaceForForrestTomlin() const
00280 {
00281 CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
00282 CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
00283 return (space>=0)&&doForrestTomlin_;
00284 }
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;
00392 int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
00393 CoinIndexedVector * regionSparse2,
00394 CoinIndexedVector * regionSparse3,
00395 bool noPermuteRegion3=false) ;
00400 int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00401 CoinIndexedVector * regionSparse2) const;
00403 void goSparse();
00405 inline int sparseThreshold ( ) const
00406 { return sparseThreshold_;}
00408 void sparseThreshold ( int value );
00410
00411
00415
00416 inline void clearArrays()
00417 { gutsOfDestructor();}
00419
00422
00425 int add ( CoinBigIndex numberElements,
00426 int indicesRow[],
00427 int indicesColumn[], double elements[] );
00428
00431 int addColumn ( CoinBigIndex numberElements,
00432 int indicesRow[], double elements[] );
00433
00436 int addRow ( CoinBigIndex numberElements,
00437 int indicesColumn[], double elements[] );
00438
00440 int deleteColumn ( int Row );
00442 int deleteRow ( int Row );
00443
00447 int replaceRow ( int whichRow, int numberElements,
00448 const int indicesColumn[], const double elements[] );
00450 void emptyRows(int numberToEmpty, const int which[]);
00452
00453
00454 void checkSparse();
00456 inline bool collectStatistics() const
00457 { return collectStatistics_;}
00459 inline void setCollectStatistics(bool onOff) const
00460 { collectStatistics_ = onOff;}
00462 void gutsOfDestructor(int type=1);
00464 void gutsOfInitialize(int type);
00465 void gutsOfCopy(const CoinFactorization &other);
00466
00468 void resetStatistics();
00469
00470
00472
00474
00475 void getAreas ( int numberRows,
00476 int numberColumns,
00477 CoinBigIndex maximumL,
00478 CoinBigIndex maximumU );
00479
00482 void preProcess ( int state,
00483 int possibleDuplicates = -1 );
00485 int factor ( );
00486 protected:
00489 int factorSparse ( );
00492 int factorDense ( );
00493
00495 bool pivotOneOtherRow ( int pivotRow,
00496 int pivotColumn );
00498 bool pivotRowSingleton ( int pivotRow,
00499 int pivotColumn );
00501 bool pivotColumnSingleton ( int pivotRow,
00502 int pivotColumn );
00503
00508 bool getColumnSpace ( int iColumn,
00509 int extraNeeded );
00510
00514 bool getColumnSpaceIterateR ( int iColumn, double value,
00515 int iRow);
00521 CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00522 int iRow);
00526 bool getRowSpace ( int iRow, int extraNeeded );
00527
00531 bool getRowSpaceIterate ( int iRow,
00532 int extraNeeded );
00534 void checkConsistency ( );
00536 inline void addLink ( int index, int count ) {
00537 int *nextCount = nextCount_.array();
00538 int *firstCount = firstCount_.array();
00539 int *lastCount = lastCount_.array();
00540 int next = firstCount[count];
00541 lastCount[index] = -2 - count;
00542 if ( next < 0 ) {
00543
00544 firstCount[count] = index;
00545 nextCount[index] = -1;
00546 } else {
00547 firstCount[count] = index;
00548 nextCount[index] = next;
00549 lastCount[next] = index;
00550 }}
00552 inline void deleteLink ( int index ) {
00553 int *nextCount = nextCount_.array();
00554 int *firstCount = firstCount_.array();
00555 int *lastCount = lastCount_.array();
00556 int next = nextCount[index];
00557 int last = lastCount[index];
00558 if ( last >= 0 ) {
00559 nextCount[last] = next;
00560 } else {
00561 int count = -last - 2;
00562
00563 firstCount[count] = next;
00564 }
00565 if ( next >= 0 ) {
00566 lastCount[next] = last;
00567 }
00568 nextCount[index] = -2;
00569 lastCount[index] = -2;
00570 return;
00571 }
00573 void separateLinks(int count,bool rowsFirst);
00575 void cleanup ( );
00576
00578 void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00580 void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00582 void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00584 void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00585
00587 void updateColumnR ( CoinIndexedVector * region ) const;
00590 void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00591
00593 void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00594
00596 void updateColumnUSparse ( CoinIndexedVector * regionSparse,
00597 int * indexIn) const;
00599 void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
00600 int * indexIn) const;
00602 int updateColumnUDensish ( double * COIN_RESTRICT region,
00603 int * COIN_RESTRICT regionIndex) const;
00605 void updateTwoColumnsUDensish (
00606 int & numberNonZero1,
00607 double * COIN_RESTRICT region1,
00608 int * COIN_RESTRICT index1,
00609 int & numberNonZero2,
00610 double * COIN_RESTRICT region2,
00611 int * COIN_RESTRICT index2) const;
00613 void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
00615 void permuteBack ( CoinIndexedVector * regionSparse,
00616 CoinIndexedVector * outVector) const;
00617
00619 void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00622 void updateColumnTransposeU ( CoinIndexedVector * region,
00623 int smallestIndex) const;
00626 void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00627 int smallestIndex) const;
00630 void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00631 int smallestIndex) const;
00634 void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00635
00637 void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00639 void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00641 void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00642
00644 void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00646 void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00648 void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00650 void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00652 void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00653 public:
00658 int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00659 int pivotRow, double alpha);
00660 protected:
00663 int checkPivot(double saveFromU, double oldPivot) const;
00664
00665 #ifdef INT_IS_8
00666 #define COINFACTORIZATION_BITS_PER_INT 64
00667 #define COINFACTORIZATION_SHIFT_PER_INT 6
00668 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00669 #else
00670 #define COINFACTORIZATION_BITS_PER_INT 32
00671 #define COINFACTORIZATION_SHIFT_PER_INT 5
00672 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00673 #endif
00674 template <class T> inline bool
00675 pivot ( int pivotRow,
00676 int pivotColumn,
00677 CoinBigIndex pivotRowPosition,
00678 CoinBigIndex pivotColumnPosition,
00679 double work[],
00680 unsigned int workArea2[],
00681 int increment,
00682 int increment2,
00683 T markRow[] ,
00684 int largeInteger)
00685 {
00686 int *indexColumnU = indexColumnU_.array();
00687 CoinBigIndex *startColumnU = startColumnU_.array();
00688 int *numberInColumn = numberInColumn_.array();
00689 double *elementU = elementU_.array();
00690 int *indexRowU = indexRowU_.array();
00691 CoinBigIndex *startRowU = startRowU_.array();
00692 int *numberInRow = numberInRow_.array();
00693 double *elementL = elementL_.array();
00694 int *indexRowL = indexRowL_.array();
00695 int *saveColumn = saveColumn_.array();
00696 int *nextRow = nextRow_.array();
00697 int *lastRow = lastRow_.array() ;
00698
00699
00700 int numberInPivotRow = numberInRow[pivotRow] - 1;
00701 CoinBigIndex startColumn = startColumnU[pivotColumn];
00702 int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00703 CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00704 int put = 0;
00705 CoinBigIndex startRow = startRowU[pivotRow];
00706 CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00707
00708 if ( pivotColumnPosition < 0 ) {
00709 for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00710 int iColumn = indexColumnU[pivotColumnPosition];
00711 if ( iColumn != pivotColumn ) {
00712 saveColumn[put++] = iColumn;
00713 } else {
00714 break;
00715 }
00716 }
00717 } else {
00718 for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00719 saveColumn[put++] = indexColumnU[i];
00720 }
00721 }
00722 assert (pivotColumnPosition<endRow);
00723 assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00724 pivotColumnPosition++;
00725 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00726 saveColumn[put++] = indexColumnU[pivotColumnPosition];
00727 }
00728
00729 int next = nextRow[pivotRow];
00730 int last = lastRow[pivotRow];
00731
00732 nextRow[last] = next;
00733 lastRow[next] = last;
00734 nextRow[pivotRow] = numberGoodU_;
00735 lastRow[pivotRow] = -2;
00736 numberInRow[pivotRow] = 0;
00737
00738 CoinBigIndex l = lengthL_;
00739
00740 if ( l + numberInPivotColumn > lengthAreaL_ ) {
00741
00742 printf("more memory needed in middle of invert\n");
00743 return false;
00744 }
00745
00746 CoinBigIndex lSave = l;
00747
00748 pivotRowL_.array()[numberGoodL_] = pivotRow;
00749 CoinBigIndex * startColumnL = startColumnL_.array();
00750 startColumnL[numberGoodL_] = l;
00751 numberGoodL_++;
00752 startColumnL[numberGoodL_] = l + numberInPivotColumn;
00753 lengthL_ += numberInPivotColumn;
00754 if ( pivotRowPosition < 0 ) {
00755 for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00756 int iRow = indexRowU[pivotRowPosition];
00757 if ( iRow != pivotRow ) {
00758 indexRowL[l] = iRow;
00759 elementL[l] = elementU[pivotRowPosition];
00760 markRow[iRow] = l - lSave;
00761 l++;
00762
00763 CoinBigIndex start = startRowU[iRow];
00764 CoinBigIndex end = start + numberInRow[iRow];
00765 CoinBigIndex where = start;
00766
00767 while ( indexColumnU[where] != pivotColumn ) {
00768 where++;
00769 }
00770 #if DEBUG_COIN
00771 if ( where >= end ) {
00772 abort ( );
00773 }
00774 #endif
00775 indexColumnU[where] = indexColumnU[end - 1];
00776 numberInRow[iRow]--;
00777 } else {
00778 break;
00779 }
00780 }
00781 } else {
00782 CoinBigIndex i;
00783
00784 for ( i = startColumn; i < pivotRowPosition; i++ ) {
00785 int iRow = indexRowU[i];
00786
00787 markRow[iRow] = l - lSave;
00788 indexRowL[l] = iRow;
00789 elementL[l] = elementU[i];
00790 l++;
00791
00792 CoinBigIndex start = startRowU[iRow];
00793 CoinBigIndex end = start + numberInRow[iRow];
00794 CoinBigIndex where = start;
00795
00796 while ( indexColumnU[where] != pivotColumn ) {
00797 where++;
00798 }
00799 #if DEBUG_COIN
00800 if ( where >= end ) {
00801 abort ( );
00802 }
00803 #endif
00804 indexColumnU[where] = indexColumnU[end - 1];
00805 numberInRow[iRow]--;
00806 assert (numberInRow[iRow]>=0);
00807 }
00808 }
00809 assert (pivotRowPosition<endColumn);
00810 assert (indexRowU[pivotRowPosition]==pivotRow);
00811 double pivotElement = elementU[pivotRowPosition];
00812 double pivotMultiplier = 1.0 / pivotElement;
00813
00814 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00815 pivotRowPosition++;
00816 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00817 int iRow = indexRowU[pivotRowPosition];
00818
00819 markRow[iRow] = l - lSave;
00820 indexRowL[l] = iRow;
00821 elementL[l] = elementU[pivotRowPosition];
00822 l++;
00823
00824 CoinBigIndex start = startRowU[iRow];
00825 CoinBigIndex end = start + numberInRow[iRow];
00826 CoinBigIndex where = start;
00827
00828 while ( indexColumnU[where] != pivotColumn ) {
00829 where++;
00830 }
00831 #if DEBUG_COIN
00832 if ( where >= end ) {
00833 abort ( );
00834 }
00835 #endif
00836 indexColumnU[where] = indexColumnU[end - 1];
00837 numberInRow[iRow]--;
00838 assert (numberInRow[iRow]>=0);
00839 }
00840 markRow[pivotRow] = largeInteger;
00841
00842 numberInColumn[pivotColumn] = 0;
00843
00844 int *indexL = &indexRowL[lSave];
00845 double *multipliersL = &elementL[lSave];
00846
00847
00848 int j;
00849
00850 for ( j = 0; j < numberInPivotColumn; j++ ) {
00851 multipliersL[j] *= pivotMultiplier;
00852 }
00853
00854 CoinBigIndex iErase;
00855 for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00856 iErase++ ) {
00857 workArea2[iErase] = 0;
00858 }
00859 CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00860 unsigned int *temp2 = workArea2;
00861 int * nextColumn = nextColumn_.array();
00862
00863
00864 int jColumn;
00865 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
00866 int iColumn = saveColumn[jColumn];
00867 CoinBigIndex startColumn = startColumnU[iColumn];
00868 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
00869 int iRow = indexRowU[startColumn];
00870 double value = elementU[startColumn];
00871 double largest;
00872 CoinBigIndex put = startColumn;
00873 CoinBigIndex positionLargest = -1;
00874 double thisPivotValue = 0.0;
00875
00876
00877 bool checkLargest;
00878 int mark = markRow[iRow];
00879
00880 if ( mark < 0 ) {
00881 largest = fabs ( value );
00882 positionLargest = put;
00883 put++;
00884 checkLargest = false;
00885 } else {
00886
00887 largest = 0.0;
00888 checkLargest = true;
00889 if ( mark != largeInteger ) {
00890
00891 work[mark] = value;
00892 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00893 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00894
00895 temp2[word] = temp2[word] | ( 1 << bit );
00896 added--;
00897 } else {
00898 thisPivotValue = value;
00899 }
00900 }
00901 CoinBigIndex i;
00902 for ( i = startColumn + 1; i < endColumn; i++ ) {
00903 iRow = indexRowU[i];
00904 value = elementU[i];
00905 int mark = markRow[iRow];
00906
00907 if ( mark < 0 ) {
00908
00909 indexRowU[put] = iRow;
00910 elementU[put] = value;
00911 if ( checkLargest ) {
00912 double absValue = fabs ( value );
00913
00914 if ( absValue > largest ) {
00915 largest = absValue;
00916 positionLargest = put;
00917 }
00918 }
00919 put++;
00920 } else if ( mark != largeInteger ) {
00921
00922 work[mark] = value;
00923 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00924 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00925
00926 temp2[word] = temp2[word] | ( 1 << bit );
00927 added--;
00928 } else {
00929 thisPivotValue = value;
00930 }
00931 }
00932
00933 elementU[put] = elementU[startColumn];
00934 indexRowU[put] = indexRowU[startColumn];
00935 if ( positionLargest == startColumn ) {
00936 positionLargest = put;
00937 }
00938 put++;
00939 elementU[startColumn] = thisPivotValue;
00940 indexRowU[startColumn] = pivotRow;
00941
00942 startColumn++;
00943 numberInColumn[iColumn] = put - startColumn;
00944 int * numberInColumnPlus = numberInColumnPlus_.array();
00945 numberInColumnPlus[iColumn]++;
00946 startColumnU[iColumn]++;
00947
00948 int next = nextColumn[iColumn];
00949 CoinBigIndex space;
00950
00951 space = startColumnU[next] - put - numberInColumnPlus[next];
00952
00953 if ( numberInPivotColumn > space ) {
00954
00955 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
00956 return false;
00957 }
00958
00959 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
00960 startColumn = startColumnU[iColumn];
00961 put = startColumn + numberInColumn[iColumn];
00962 }
00963 double tolerance = zeroTolerance_;
00964
00965 int *nextCount = nextCount_.array();
00966 for ( j = 0; j < numberInPivotColumn; j++ ) {
00967 value = work[j] - thisPivotValue * multipliersL[j];
00968 double absValue = fabs ( value );
00969
00970 if ( absValue > tolerance ) {
00971 work[j] = 0.0;
00972 elementU[put] = value;
00973 indexRowU[put] = indexL[j];
00974 if ( absValue > largest ) {
00975 largest = absValue;
00976 positionLargest = put;
00977 }
00978 put++;
00979 } else {
00980 work[j] = 0.0;
00981 added--;
00982 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
00983 int bit = j & COINFACTORIZATION_MASK_PER_INT;
00984
00985 if ( temp2[word] & ( 1 << bit ) ) {
00986
00987 iRow = indexL[j];
00988 CoinBigIndex start = startRowU[iRow];
00989 CoinBigIndex end = start + numberInRow[iRow];
00990 CoinBigIndex where = start;
00991
00992 while ( indexColumnU[where] != iColumn ) {
00993 where++;
00994 }
00995 #if DEBUG_COIN
00996 if ( where >= end ) {
00997 abort ( );
00998 }
00999 #endif
01000 indexColumnU[where] = indexColumnU[end - 1];
01001 numberInRow[iRow]--;
01002 } else {
01003
01004 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01005 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01006
01007 temp2[word] = temp2[word] | ( 1 << bit );
01008 }
01009 }
01010 }
01011 numberInColumn[iColumn] = put - startColumn;
01012
01013 if ( positionLargest >= 0 ) {
01014 value = elementU[positionLargest];
01015 iRow = indexRowU[positionLargest];
01016 elementU[positionLargest] = elementU[startColumn];
01017 indexRowU[positionLargest] = indexRowU[startColumn];
01018 elementU[startColumn] = value;
01019 indexRowU[startColumn] = iRow;
01020 }
01021
01022 if ( nextCount[iColumn + numberRows_] != -2 ) {
01023
01024 deleteLink ( iColumn + numberRows_ );
01025 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01026 }
01027 temp2 += increment2;
01028 }
01029
01030 unsigned int *putBase = workArea2;
01031 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01032 int i = 0;
01033
01034
01035 while ( bigLoops ) {
01036 bigLoops--;
01037 int bit;
01038 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01039 unsigned int *putThis = putBase;
01040 int iRow = indexL[i];
01041
01042
01043 int number = 0;
01044 int jColumn;
01045
01046 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01047 unsigned int test = *putThis;
01048
01049 putThis += increment2;
01050 test = 1 - ( ( test >> bit ) & 1 );
01051 number += test;
01052 }
01053 int next = nextRow[iRow];
01054 CoinBigIndex space;
01055
01056 space = startRowU[next] - startRowU[iRow];
01057 number += numberInRow[iRow];
01058 if ( space < number ) {
01059 if ( !getRowSpace ( iRow, number ) ) {
01060 return false;
01061 }
01062 }
01063
01064 putThis = putBase;
01065 next = nextRow[iRow];
01066 number = numberInRow[iRow];
01067 CoinBigIndex end = startRowU[iRow] + number;
01068 int saveIndex = indexColumnU[startRowU[next]];
01069
01070
01071 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01072 unsigned int test = *putThis;
01073
01074 putThis += increment2;
01075 test = 1 - ( ( test >> bit ) & 1 );
01076 indexColumnU[end] = saveColumn[jColumn];
01077 end += test;
01078 }
01079
01080 indexColumnU[startRowU[next]] = saveIndex;
01081 markRow[iRow] = -1;
01082 number = end - startRowU[iRow];
01083 numberInRow[iRow] = number;
01084 deleteLink ( iRow );
01085 addLink ( iRow, number );
01086 }
01087 putBase++;
01088 }
01089 int bit;
01090
01091 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01092 unsigned int *putThis = putBase;
01093 int iRow = indexL[i];
01094
01095
01096 int number = 0;
01097 int jColumn;
01098
01099 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01100 unsigned int test = *putThis;
01101
01102 putThis += increment2;
01103 test = 1 - ( ( test >> bit ) & 1 );
01104 number += test;
01105 }
01106 int next = nextRow[iRow];
01107 CoinBigIndex space;
01108
01109 space = startRowU[next] - startRowU[iRow];
01110 number += numberInRow[iRow];
01111 if ( space < number ) {
01112 if ( !getRowSpace ( iRow, number ) ) {
01113 return false;
01114 }
01115 }
01116
01117 putThis = putBase;
01118 next = nextRow[iRow];
01119 number = numberInRow[iRow];
01120 CoinBigIndex end = startRowU[iRow] + number;
01121 int saveIndex;
01122
01123 saveIndex = indexColumnU[startRowU[next]];
01124
01125
01126 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01127 unsigned int test = *putThis;
01128
01129 putThis += increment2;
01130 test = 1 - ( ( test >> bit ) & 1 );
01131
01132 indexColumnU[end] = saveColumn[jColumn];
01133 end += test;
01134 }
01135 indexColumnU[startRowU[next]] = saveIndex;
01136 markRow[iRow] = -1;
01137 number = end - startRowU[iRow];
01138 numberInRow[iRow] = number;
01139 deleteLink ( iRow );
01140 addLink ( iRow, number );
01141 }
01142 markRow[pivotRow] = -1;
01143
01144 deleteLink ( pivotRow );
01145 deleteLink ( pivotColumn + numberRows_ );
01146 totalElements_ += added;
01147 return true;
01148 }
01149
01150
01152
01153 protected:
01154
01157
01158 double pivotTolerance_;
01160 double zeroTolerance_;
01161 #ifndef COIN_FAST_CODE
01163 double slackValue_;
01164 #else
01165 #ifndef slackValue_
01166 #define slackValue_ -1.0
01167 #endif
01168 #endif
01170 double areaFactor_;
01172 double relaxCheck_;
01174 int numberRows_;
01176 int numberRowsExtra_;
01178 int maximumRowsExtra_;
01180 int numberColumns_;
01182 int numberColumnsExtra_;
01184 int maximumColumnsExtra_;
01186 int numberGoodU_;
01188 int numberGoodL_;
01190 int maximumPivots_;
01192 int numberPivots_;
01195 CoinBigIndex totalElements_;
01197 CoinBigIndex factorElements_;
01199 CoinIntArrayWithLength pivotColumn_;
01201 CoinIntArrayWithLength permute_;
01203 CoinIntArrayWithLength permuteBack_;
01205 CoinIntArrayWithLength pivotColumnBack_;
01207 int status_;
01208
01213
01214
01216 int numberTrials_;
01218 CoinBigIndexArrayWithLength startRowU_;
01219
01221 CoinIntArrayWithLength numberInRow_;
01222
01224 CoinIntArrayWithLength numberInColumn_;
01225
01227 CoinIntArrayWithLength numberInColumnPlus_;
01228
01231 CoinIntArrayWithLength firstCount_;
01232
01234 CoinIntArrayWithLength nextCount_;
01235
01237 CoinIntArrayWithLength lastCount_;
01238
01240 CoinIntArrayWithLength nextColumn_;
01241
01243 CoinIntArrayWithLength lastColumn_;
01244
01246 CoinIntArrayWithLength nextRow_;
01247
01249 CoinIntArrayWithLength lastRow_;
01250
01252 CoinIntArrayWithLength saveColumn_;
01253
01255 CoinIntArrayWithLength markRow_;
01256
01258 int messageLevel_;
01259
01261 int biggerDimension_;
01262
01264 CoinIntArrayWithLength indexColumnU_;
01265
01267 CoinIntArrayWithLength pivotRowL_;
01268
01270 CoinDoubleArrayWithLength pivotRegion_;
01271
01273 int numberSlacks_;
01274
01276 int numberU_;
01277
01279 CoinBigIndex maximumU_;
01280
01282
01283
01285 CoinBigIndex lengthU_;
01286
01288 CoinBigIndex lengthAreaU_;
01289
01291 CoinDoubleArrayWithLength elementU_;
01292
01294 CoinIntArrayWithLength indexRowU_;
01295
01297 CoinBigIndexArrayWithLength startColumnU_;
01298
01300 CoinBigIndexArrayWithLength convertRowToColumnU_;
01301
01303 CoinBigIndex numberL_;
01304
01306 CoinBigIndex baseL_;
01307
01309 CoinBigIndex lengthL_;
01310
01312 CoinBigIndex lengthAreaL_;
01313
01315 CoinDoubleArrayWithLength elementL_;
01316
01318 CoinIntArrayWithLength indexRowL_;
01319
01321 CoinBigIndexArrayWithLength startColumnL_;
01322
01324 bool doForrestTomlin_;
01325
01327 int numberR_;
01328
01330 CoinBigIndex lengthR_;
01331
01333 CoinBigIndex lengthAreaR_;
01334
01336 double *elementR_;
01337
01339 int *indexRowR_;
01340
01342 CoinBigIndexArrayWithLength startColumnR_;
01343
01345 double * denseArea_;
01346
01348 int * densePermute_;
01349
01351 int numberDense_;
01352
01354 int denseThreshold_;
01355
01357 CoinDoubleArrayWithLength workArea_;
01358
01360 CoinUnsignedIntArrayWithLength workArea2_;
01361
01363 CoinBigIndex numberCompressions_;
01364
01366 mutable double ftranCountInput_;
01367 mutable double ftranCountAfterL_;
01368 mutable double ftranCountAfterR_;
01369 mutable double ftranCountAfterU_;
01370 mutable double btranCountInput_;
01371 mutable double btranCountAfterU_;
01372 mutable double btranCountAfterR_;
01373 mutable double btranCountAfterL_;
01374
01376 mutable int numberFtranCounts_;
01377 mutable int numberBtranCounts_;
01378
01380 double ftranAverageAfterL_;
01381 double ftranAverageAfterR_;
01382 double ftranAverageAfterU_;
01383 double btranAverageAfterU_;
01384 double btranAverageAfterR_;
01385 double btranAverageAfterL_;
01386
01388 mutable bool collectStatistics_;
01389
01391 int sparseThreshold_;
01392
01394 int sparseThreshold2_;
01395
01397 CoinBigIndexArrayWithLength startRowL_;
01398
01400 CoinIntArrayWithLength indexColumnL_;
01401
01403 CoinDoubleArrayWithLength elementByRowL_;
01404
01406 mutable CoinIntArrayWithLength sparse_;
01410 int biasLU_;
01416 int persistenceFlag_;
01418 };
01419
01420 #ifdef COIN_HAS_LAPACK
01421 #define DENSE_CODE 1
01422
01423 #ifndef ipfint
01424
01425 typedef int ipfint;
01426 typedef const int cipfint;
01427 #endif
01428 #endif
01429 #endif