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 <cstdio>
00017 #include "CoinFinite.hpp"
00018 #include "CoinIndexedVector.hpp"
00019 class CoinPackedMatrix;
00046 class CoinFactorization {
00047 friend void CoinFactorizationUnitTest( const std::string & mpsDir );
00048
00049 public:
00050
00053
00054 CoinFactorization ( );
00056 CoinFactorization ( const CoinFactorization &other);
00057
00059 ~CoinFactorization ( );
00061 void almostDestructor();
00063 void show_self ( ) const;
00065 int saveFactorization (const char * file ) const;
00069 int restoreFactorization (const char * file , bool factor=false) ;
00071 void sort ( ) const;
00073 CoinFactorization & operator = ( const CoinFactorization & other );
00075
00085 int factorize ( const CoinPackedMatrix & matrix,
00086 int rowIsBasic[], int columnIsBasic[] ,
00087 double areaFactor = 0.0 );
00098 int factorize ( int numberRows,
00099 int numberColumns,
00100 CoinBigIndex numberElements,
00101 CoinBigIndex maximumL,
00102 CoinBigIndex maximumU,
00103 const int indicesRow[],
00104 const int indicesColumn[], const double elements[] ,
00105 int permutation[],
00106 double areaFactor = 0.0);
00111 int factorizePart1 ( int numberRows,
00112 int numberColumns,
00113 CoinBigIndex estimateNumberElements,
00114 int * indicesRow[],
00115 int * indicesColumn[],
00116 double * elements[],
00117 double areaFactor = 0.0);
00124 int factorizePart2 (int permutation[],int exactNumberElements);
00126 double conditionNumber() const;
00127
00129
00132
00133 inline int status ( ) const {
00134 return status_;
00135 }
00137 inline void setStatus ( int value)
00138 { status_=value; }
00140 inline int pivots ( ) const {
00141 return numberPivots_;
00142 }
00144 inline void setPivots ( int value )
00145 { numberPivots_=value; }
00147 inline int *permute ( ) const {
00148 return permute_.array();
00149 }
00151 inline int *pivotColumn ( ) const {
00152 return pivotColumn_.array();
00153 }
00155 inline double *pivotRegion ( ) const {
00156 return pivotRegion_.array();
00157 }
00159 inline int *permuteBack ( ) const {
00160 return permuteBack_.array();
00161 }
00163 inline int *pivotColumnBack ( ) const {
00164 return pivotColumnBack_.array();
00165 }
00167 inline CoinBigIndex * startRowL() const
00168 { return startRowL_.array();}
00169
00171 inline CoinBigIndex * startColumnL() const
00172 { return startColumnL_.array();}
00173
00175 inline int * indexColumnL() const
00176 { return indexColumnL_.array();}
00177
00179 inline int * indexRowL() const
00180 { return indexRowL_.array();}
00181
00183 inline double * elementByRowL() const
00184 { return elementByRowL_.array();}
00185
00187 inline int numberRowsExtra ( ) const {
00188 return numberRowsExtra_;
00189 }
00191 inline void setNumberRows(int value)
00192 { numberRows_ = value; }
00194 inline int numberRows ( ) const {
00195 return numberRows_;
00196 }
00198 inline CoinBigIndex numberL() const
00199 { return numberL_;}
00200
00202 inline CoinBigIndex baseL() const
00203 { return baseL_;}
00205 inline int maximumRowsExtra ( ) const {
00206 return maximumRowsExtra_;
00207 }
00209 inline int numberColumns ( ) const {
00210 return numberColumns_;
00211 }
00213 inline int numberElements ( ) const {
00214 return totalElements_;
00215 }
00217 inline int numberForrestTomlin ( ) const {
00218 return numberInColumn_.array()[numberColumnsExtra_];
00219 }
00221 inline int numberGoodColumns ( ) const {
00222 return numberGoodU_;
00223 }
00225 inline double areaFactor ( ) const {
00226 return areaFactor_;
00227 }
00228 inline void areaFactor ( double value ) {
00229 areaFactor_=value;
00230 }
00232 double adjustedAreaFactor() const;
00234 inline void relaxAccuracyCheck(double value)
00235 { relaxCheck_ = value;}
00236 inline double getAccuracyCheck() const
00237 { return relaxCheck_;}
00239 inline int messageLevel ( ) const {
00240 return messageLevel_ ;
00241 }
00242 void messageLevel ( int value );
00244 inline int maximumPivots ( ) const {
00245 return maximumPivots_ ;
00246 }
00247 void maximumPivots ( int value );
00248
00250 inline int denseThreshold() const
00251 { return denseThreshold_;}
00253 inline void setDenseThreshold(int value)
00254 { denseThreshold_ = value;}
00256 inline double pivotTolerance ( ) const {
00257 return pivotTolerance_ ;
00258 }
00259 void pivotTolerance ( double value );
00261 inline double zeroTolerance ( ) const {
00262 return zeroTolerance_ ;
00263 }
00264 void zeroTolerance ( double value );
00265 #ifndef COIN_FAST_CODE
00267 inline double slackValue ( ) const {
00268 return slackValue_ ;
00269 }
00270 void slackValue ( double value );
00271 #endif
00273 double maximumCoefficient() const;
00275 inline bool forrestTomlin() const
00276 { return doForrestTomlin_;}
00277 inline void setForrestTomlin(bool value)
00278 { doForrestTomlin_=value;}
00280 inline bool spaceForForrestTomlin() const
00281 {
00282 CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
00283 CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
00284 return (space>=0)&&doForrestTomlin_;
00285 }
00287
00290
00292 inline int numberDense() const
00293 { return numberDense_;}
00294
00296 inline CoinBigIndex numberElementsU ( ) const {
00297 return lengthU_;
00298 }
00300 inline void setNumberElementsU(CoinBigIndex value)
00301 { lengthU_ = value; }
00303 inline CoinBigIndex lengthAreaU ( ) const {
00304 return lengthAreaU_;
00305 }
00307 inline CoinBigIndex numberElementsL ( ) const {
00308 return lengthL_;
00309 }
00311 inline CoinBigIndex lengthAreaL ( ) const {
00312 return lengthAreaL_;
00313 }
00315 inline CoinBigIndex numberElementsR ( ) const {
00316 return lengthR_;
00317 }
00319 inline CoinBigIndex numberCompressions() const
00320 { return numberCompressions_;}
00322 inline int * numberInRow() const
00323 { return numberInRow_.array();}
00325 inline int * numberInColumn() const
00326 { return numberInColumn_.array();}
00328 inline double * elementU() const
00329 { return elementU_.array();}
00331 inline int * indexRowU() const
00332 { return indexRowU_.array();}
00334 inline CoinBigIndex * startColumnU() const
00335 { return startColumnU_.array();}
00337 inline int maximumColumnsExtra()
00338 { return maximumColumnsExtra_;}
00342 inline int biasLU() const
00343 { return biasLU_;}
00344 inline void setBiasLU(int value)
00345 { biasLU_=value;}
00351 inline int persistenceFlag() const
00352 { return persistenceFlag_;}
00353 void setPersistenceFlag(int value);
00355
00358
00366 int replaceColumn ( CoinIndexedVector * regionSparse,
00367 int pivotRow,
00368 double pivotCheck ,
00369 bool checkBeforeModifying=false);
00371
00381 int updateColumnFT ( CoinIndexedVector * regionSparse,
00382 CoinIndexedVector * regionSparse2);
00385 int updateColumn ( CoinIndexedVector * regionSparse,
00386 CoinIndexedVector * regionSparse2,
00387 bool noPermute=false) const;
00393 int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
00394 CoinIndexedVector * regionSparse2,
00395 CoinIndexedVector * regionSparse3,
00396 bool noPermuteRegion3=false) ;
00401 int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00402 CoinIndexedVector * regionSparse2) const;
00404 void goSparse();
00406 inline int sparseThreshold ( ) const
00407 { return sparseThreshold_;}
00409 void sparseThreshold ( int value );
00411
00412
00416
00417 inline void clearArrays()
00418 { gutsOfDestructor();}
00420
00423
00426 int add ( CoinBigIndex numberElements,
00427 int indicesRow[],
00428 int indicesColumn[], double elements[] );
00429
00432 int addColumn ( CoinBigIndex numberElements,
00433 int indicesRow[], double elements[] );
00434
00437 int addRow ( CoinBigIndex numberElements,
00438 int indicesColumn[], double elements[] );
00439
00441 int deleteColumn ( int Row );
00443 int deleteRow ( int Row );
00444
00448 int replaceRow ( int whichRow, int numberElements,
00449 const int indicesColumn[], const double elements[] );
00451 void emptyRows(int numberToEmpty, const int which[]);
00453
00454
00455 void checkSparse();
00457 inline bool collectStatistics() const
00458 { return collectStatistics_;}
00460 inline void setCollectStatistics(bool onOff) const
00461 { collectStatistics_ = onOff;}
00463 void gutsOfDestructor(int type=1);
00465 void gutsOfInitialize(int type);
00466 void gutsOfCopy(const CoinFactorization &other);
00467
00469 void resetStatistics();
00470
00471
00473
00475
00476 void getAreas ( int numberRows,
00477 int numberColumns,
00478 CoinBigIndex maximumL,
00479 CoinBigIndex maximumU );
00480
00483 void preProcess ( int state,
00484 int possibleDuplicates = -1 );
00486 int factor ( );
00487 protected:
00490 int factorSparse ( );
00493 int factorDense ( );
00494
00496 bool pivotOneOtherRow ( int pivotRow,
00497 int pivotColumn );
00499 bool pivotRowSingleton ( int pivotRow,
00500 int pivotColumn );
00502 bool pivotColumnSingleton ( int pivotRow,
00503 int pivotColumn );
00504
00509 bool getColumnSpace ( int iColumn,
00510 int extraNeeded );
00511
00515 bool getColumnSpaceIterateR ( int iColumn, double value,
00516 int iRow);
00522 CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00523 int iRow);
00527 bool getRowSpace ( int iRow, int extraNeeded );
00528
00532 bool getRowSpaceIterate ( int iRow,
00533 int extraNeeded );
00535 void checkConsistency ( );
00537 inline void addLink ( int index, int count ) {
00538 int *nextCount = nextCount_.array();
00539 int *firstCount = firstCount_.array();
00540 int *lastCount = lastCount_.array();
00541 int next = firstCount[count];
00542 lastCount[index] = -2 - count;
00543 if ( next < 0 ) {
00544
00545 firstCount[count] = index;
00546 nextCount[index] = -1;
00547 } else {
00548 firstCount[count] = index;
00549 nextCount[index] = next;
00550 lastCount[next] = index;
00551 }}
00553 inline void deleteLink ( int index ) {
00554 int *nextCount = nextCount_.array();
00555 int *firstCount = firstCount_.array();
00556 int *lastCount = lastCount_.array();
00557 int next = nextCount[index];
00558 int last = lastCount[index];
00559 if ( last >= 0 ) {
00560 nextCount[last] = next;
00561 } else {
00562 int count = -last - 2;
00563
00564 firstCount[count] = next;
00565 }
00566 if ( next >= 0 ) {
00567 lastCount[next] = last;
00568 }
00569 nextCount[index] = -2;
00570 lastCount[index] = -2;
00571 return;
00572 }
00574 void separateLinks(int count,bool rowsFirst);
00576 void cleanup ( );
00577
00579 void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00581 void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00583 void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00585 void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00586
00588 void updateColumnR ( CoinIndexedVector * region ) const;
00591 void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00592
00594 void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00595
00597 void updateColumnUSparse ( CoinIndexedVector * regionSparse,
00598 int * indexIn) const;
00600 void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
00601 int * indexIn) const;
00603 int updateColumnUDensish ( double * COIN_RESTRICT region,
00604 int * COIN_RESTRICT regionIndex) const;
00606 void updateTwoColumnsUDensish (
00607 int & numberNonZero1,
00608 double * COIN_RESTRICT region1,
00609 int * COIN_RESTRICT index1,
00610 int & numberNonZero2,
00611 double * COIN_RESTRICT region2,
00612 int * COIN_RESTRICT index2) const;
00614 void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
00616 void permuteBack ( CoinIndexedVector * regionSparse,
00617 CoinIndexedVector * outVector) const;
00618
00620 void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00623 void updateColumnTransposeU ( CoinIndexedVector * region,
00624 int smallestIndex) const;
00627 void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00628 int smallestIndex) const;
00631 void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00632 int smallestIndex) const;
00635 void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00636
00638 void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00640 void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00642 void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00643
00645 void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00647 void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00649 void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00651 void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00653 void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00654 public:
00659 int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00660 int pivotRow, double alpha);
00661 protected:
00664 int checkPivot(double saveFromU, double oldPivot) const;
00665
00666 #ifdef INT_IS_8
00667 #define COINFACTORIZATION_BITS_PER_INT 64
00668 #define COINFACTORIZATION_SHIFT_PER_INT 6
00669 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00670 #else
00671 #define COINFACTORIZATION_BITS_PER_INT 32
00672 #define COINFACTORIZATION_SHIFT_PER_INT 5
00673 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00674 #endif
00675 template <class T> inline bool
00676 pivot ( int pivotRow,
00677 int pivotColumn,
00678 CoinBigIndex pivotRowPosition,
00679 CoinBigIndex pivotColumnPosition,
00680 double work[],
00681 unsigned int workArea2[],
00682 int increment,
00683 int increment2,
00684 T markRow[] ,
00685 int largeInteger)
00686 {
00687 int *indexColumnU = indexColumnU_.array();
00688 CoinBigIndex *startColumnU = startColumnU_.array();
00689 int *numberInColumn = numberInColumn_.array();
00690 double *elementU = elementU_.array();
00691 int *indexRowU = indexRowU_.array();
00692 CoinBigIndex *startRowU = startRowU_.array();
00693 int *numberInRow = numberInRow_.array();
00694 double *elementL = elementL_.array();
00695 int *indexRowL = indexRowL_.array();
00696 int *saveColumn = saveColumn_.array();
00697 int *nextRow = nextRow_.array();
00698 int *lastRow = lastRow_.array() ;
00699
00700
00701 int numberInPivotRow = numberInRow[pivotRow] - 1;
00702 CoinBigIndex startColumn = startColumnU[pivotColumn];
00703 int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00704 CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00705 int put = 0;
00706 CoinBigIndex startRow = startRowU[pivotRow];
00707 CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00708
00709 if ( pivotColumnPosition < 0 ) {
00710 for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00711 int iColumn = indexColumnU[pivotColumnPosition];
00712 if ( iColumn != pivotColumn ) {
00713 saveColumn[put++] = iColumn;
00714 } else {
00715 break;
00716 }
00717 }
00718 } else {
00719 for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00720 saveColumn[put++] = indexColumnU[i];
00721 }
00722 }
00723 assert (pivotColumnPosition<endRow);
00724 assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00725 pivotColumnPosition++;
00726 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00727 saveColumn[put++] = indexColumnU[pivotColumnPosition];
00728 }
00729
00730 int next = nextRow[pivotRow];
00731 int last = lastRow[pivotRow];
00732
00733 nextRow[last] = next;
00734 lastRow[next] = last;
00735 nextRow[pivotRow] = numberGoodU_;
00736 lastRow[pivotRow] = -2;
00737 numberInRow[pivotRow] = 0;
00738
00739 CoinBigIndex l = lengthL_;
00740
00741 if ( l + numberInPivotColumn > lengthAreaL_ ) {
00742
00743 printf("more memory needed in middle of invert\n");
00744 return false;
00745 }
00746
00747 CoinBigIndex lSave = l;
00748
00749 pivotRowL_.array()[numberGoodL_] = pivotRow;
00750 CoinBigIndex * startColumnL = startColumnL_.array();
00751 startColumnL[numberGoodL_] = l;
00752 numberGoodL_++;
00753 startColumnL[numberGoodL_] = l + numberInPivotColumn;
00754 lengthL_ += numberInPivotColumn;
00755 if ( pivotRowPosition < 0 ) {
00756 for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00757 int iRow = indexRowU[pivotRowPosition];
00758 if ( iRow != pivotRow ) {
00759 indexRowL[l] = iRow;
00760 elementL[l] = elementU[pivotRowPosition];
00761 markRow[iRow] = static_cast<short>(l - lSave);
00762 l++;
00763
00764 CoinBigIndex start = startRowU[iRow];
00765 CoinBigIndex end = start + numberInRow[iRow];
00766 CoinBigIndex where = start;
00767
00768 while ( indexColumnU[where] != pivotColumn ) {
00769 where++;
00770 }
00771 #if DEBUG_COIN
00772 if ( where >= end ) {
00773 abort ( );
00774 }
00775 #endif
00776 indexColumnU[where] = indexColumnU[end - 1];
00777 numberInRow[iRow]--;
00778 } else {
00779 break;
00780 }
00781 }
00782 } else {
00783 CoinBigIndex i;
00784
00785 for ( i = startColumn; i < pivotRowPosition; i++ ) {
00786 int iRow = indexRowU[i];
00787
00788 markRow[iRow] = static_cast<short>(l - lSave);
00789 indexRowL[l] = iRow;
00790 elementL[l] = elementU[i];
00791 l++;
00792
00793 CoinBigIndex start = startRowU[iRow];
00794 CoinBigIndex end = start + numberInRow[iRow];
00795 CoinBigIndex where = start;
00796
00797 while ( indexColumnU[where] != pivotColumn ) {
00798 where++;
00799 }
00800 #if DEBUG_COIN
00801 if ( where >= end ) {
00802 abort ( );
00803 }
00804 #endif
00805 indexColumnU[where] = indexColumnU[end - 1];
00806 numberInRow[iRow]--;
00807 assert (numberInRow[iRow]>=0);
00808 }
00809 }
00810 assert (pivotRowPosition<endColumn);
00811 assert (indexRowU[pivotRowPosition]==pivotRow);
00812 double pivotElement = elementU[pivotRowPosition];
00813 double pivotMultiplier = 1.0 / pivotElement;
00814
00815 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00816 pivotRowPosition++;
00817 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00818 int iRow = indexRowU[pivotRowPosition];
00819
00820 markRow[iRow] = static_cast<short>(l - lSave);
00821 indexRowL[l] = iRow;
00822 elementL[l] = elementU[pivotRowPosition];
00823 l++;
00824
00825 CoinBigIndex start = startRowU[iRow];
00826 CoinBigIndex end = start + numberInRow[iRow];
00827 CoinBigIndex where = start;
00828
00829 while ( indexColumnU[where] != pivotColumn ) {
00830 where++;
00831 }
00832 #if DEBUG_COIN
00833 if ( where >= end ) {
00834 abort ( );
00835 }
00836 #endif
00837 indexColumnU[where] = indexColumnU[end - 1];
00838 numberInRow[iRow]--;
00839 assert (numberInRow[iRow]>=0);
00840 }
00841 markRow[pivotRow] = static_cast<short>(largeInteger);
00842
00843 numberInColumn[pivotColumn] = 0;
00844
00845 int *indexL = &indexRowL[lSave];
00846 double *multipliersL = &elementL[lSave];
00847
00848
00849 int j;
00850
00851 for ( j = 0; j < numberInPivotColumn; j++ ) {
00852 multipliersL[j] *= pivotMultiplier;
00853 }
00854
00855 CoinBigIndex iErase;
00856 for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00857 iErase++ ) {
00858 workArea2[iErase] = 0;
00859 }
00860 CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00861 unsigned int *temp2 = workArea2;
00862 int * nextColumn = nextColumn_.array();
00863
00864
00865 int jColumn;
00866 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
00867 int iColumn = saveColumn[jColumn];
00868 CoinBigIndex startColumn = startColumnU[iColumn];
00869 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
00870 int iRow = indexRowU[startColumn];
00871 double value = elementU[startColumn];
00872 double largest;
00873 CoinBigIndex put = startColumn;
00874 CoinBigIndex positionLargest = -1;
00875 double thisPivotValue = 0.0;
00876
00877
00878 bool checkLargest;
00879 int mark = markRow[iRow];
00880
00881 if ( mark < 0 ) {
00882 largest = fabs ( value );
00883 positionLargest = put;
00884 put++;
00885 checkLargest = false;
00886 } else {
00887
00888 largest = 0.0;
00889 checkLargest = true;
00890 if ( mark != largeInteger ) {
00891
00892 work[mark] = value;
00893 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00894 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00895
00896 temp2[word] = temp2[word] | ( 1 << bit );
00897 added--;
00898 } else {
00899 thisPivotValue = value;
00900 }
00901 }
00902 CoinBigIndex i;
00903 for ( i = startColumn + 1; i < endColumn; i++ ) {
00904 iRow = indexRowU[i];
00905 value = elementU[i];
00906 int mark = markRow[iRow];
00907
00908 if ( mark < 0 ) {
00909
00910 indexRowU[put] = iRow;
00911 elementU[put] = value;
00912 if ( checkLargest ) {
00913 double absValue = fabs ( value );
00914
00915 if ( absValue > largest ) {
00916 largest = absValue;
00917 positionLargest = put;
00918 }
00919 }
00920 put++;
00921 } else if ( mark != largeInteger ) {
00922
00923 work[mark] = value;
00924 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00925 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00926
00927 temp2[word] = temp2[word] | ( 1 << bit );
00928 added--;
00929 } else {
00930 thisPivotValue = value;
00931 }
00932 }
00933
00934 elementU[put] = elementU[startColumn];
00935 indexRowU[put] = indexRowU[startColumn];
00936 if ( positionLargest == startColumn ) {
00937 positionLargest = put;
00938 }
00939 put++;
00940 elementU[startColumn] = thisPivotValue;
00941 indexRowU[startColumn] = pivotRow;
00942
00943 startColumn++;
00944 numberInColumn[iColumn] = put - startColumn;
00945 int * numberInColumnPlus = numberInColumnPlus_.array();
00946 numberInColumnPlus[iColumn]++;
00947 startColumnU[iColumn]++;
00948
00949 int next = nextColumn[iColumn];
00950 CoinBigIndex space;
00951
00952 space = startColumnU[next] - put - numberInColumnPlus[next];
00953
00954 if ( numberInPivotColumn > space ) {
00955
00956 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
00957 return false;
00958 }
00959
00960 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
00961 startColumn = startColumnU[iColumn];
00962 put = startColumn + numberInColumn[iColumn];
00963 }
00964 double tolerance = zeroTolerance_;
00965
00966 int *nextCount = nextCount_.array();
00967 for ( j = 0; j < numberInPivotColumn; j++ ) {
00968 value = work[j] - thisPivotValue * multipliersL[j];
00969 double absValue = fabs ( value );
00970
00971 if ( absValue > tolerance ) {
00972 work[j] = 0.0;
00973 elementU[put] = value;
00974 indexRowU[put] = indexL[j];
00975 if ( absValue > largest ) {
00976 largest = absValue;
00977 positionLargest = put;
00978 }
00979 put++;
00980 } else {
00981 work[j] = 0.0;
00982 added--;
00983 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
00984 int bit = j & COINFACTORIZATION_MASK_PER_INT;
00985
00986 if ( temp2[word] & ( 1 << bit ) ) {
00987
00988 iRow = indexL[j];
00989 CoinBigIndex start = startRowU[iRow];
00990 CoinBigIndex end = start + numberInRow[iRow];
00991 CoinBigIndex where = start;
00992
00993 while ( indexColumnU[where] != iColumn ) {
00994 where++;
00995 }
00996 #if DEBUG_COIN
00997 if ( where >= end ) {
00998 abort ( );
00999 }
01000 #endif
01001 indexColumnU[where] = indexColumnU[end - 1];
01002 numberInRow[iRow]--;
01003 } else {
01004
01005 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01006 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01007
01008 temp2[word] = temp2[word] | ( 1 << bit );
01009 }
01010 }
01011 }
01012 numberInColumn[iColumn] = put - startColumn;
01013
01014 if ( positionLargest >= 0 ) {
01015 value = elementU[positionLargest];
01016 iRow = indexRowU[positionLargest];
01017 elementU[positionLargest] = elementU[startColumn];
01018 indexRowU[positionLargest] = indexRowU[startColumn];
01019 elementU[startColumn] = value;
01020 indexRowU[startColumn] = iRow;
01021 }
01022
01023 if ( nextCount[iColumn + numberRows_] != -2 ) {
01024
01025 deleteLink ( iColumn + numberRows_ );
01026 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01027 }
01028 temp2 += increment2;
01029 }
01030
01031 unsigned int *putBase = workArea2;
01032 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01033 int i = 0;
01034
01035
01036 while ( bigLoops ) {
01037 bigLoops--;
01038 int bit;
01039 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01040 unsigned int *putThis = putBase;
01041 int iRow = indexL[i];
01042
01043
01044 int number = 0;
01045 int jColumn;
01046
01047 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01048 unsigned int test = *putThis;
01049
01050 putThis += increment2;
01051 test = 1 - ( ( test >> bit ) & 1 );
01052 number += test;
01053 }
01054 int next = nextRow[iRow];
01055 CoinBigIndex space;
01056
01057 space = startRowU[next] - startRowU[iRow];
01058 number += numberInRow[iRow];
01059 if ( space < number ) {
01060 if ( !getRowSpace ( iRow, number ) ) {
01061 return false;
01062 }
01063 }
01064
01065 putThis = putBase;
01066 next = nextRow[iRow];
01067 number = numberInRow[iRow];
01068 CoinBigIndex end = startRowU[iRow] + number;
01069 int saveIndex = indexColumnU[startRowU[next]];
01070
01071
01072 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01073 unsigned int test = *putThis;
01074
01075 putThis += increment2;
01076 test = 1 - ( ( test >> bit ) & 1 );
01077 indexColumnU[end] = saveColumn[jColumn];
01078 end += test;
01079 }
01080
01081 indexColumnU[startRowU[next]] = saveIndex;
01082 markRow[iRow] = -1;
01083 number = end - startRowU[iRow];
01084 numberInRow[iRow] = number;
01085 deleteLink ( iRow );
01086 addLink ( iRow, number );
01087 }
01088 putBase++;
01089 }
01090 int bit;
01091
01092 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01093 unsigned int *putThis = putBase;
01094 int iRow = indexL[i];
01095
01096
01097 int number = 0;
01098 int jColumn;
01099
01100 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01101 unsigned int test = *putThis;
01102
01103 putThis += increment2;
01104 test = 1 - ( ( test >> bit ) & 1 );
01105 number += test;
01106 }
01107 int next = nextRow[iRow];
01108 CoinBigIndex space;
01109
01110 space = startRowU[next] - startRowU[iRow];
01111 number += numberInRow[iRow];
01112 if ( space < number ) {
01113 if ( !getRowSpace ( iRow, number ) ) {
01114 return false;
01115 }
01116 }
01117
01118 putThis = putBase;
01119 next = nextRow[iRow];
01120 number = numberInRow[iRow];
01121 CoinBigIndex end = startRowU[iRow] + number;
01122 int saveIndex;
01123
01124 saveIndex = indexColumnU[startRowU[next]];
01125
01126
01127 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01128 unsigned int test = *putThis;
01129
01130 putThis += increment2;
01131 test = 1 - ( ( test >> bit ) & 1 );
01132
01133 indexColumnU[end] = saveColumn[jColumn];
01134 end += test;
01135 }
01136 indexColumnU[startRowU[next]] = saveIndex;
01137 markRow[iRow] = -1;
01138 number = end - startRowU[iRow];
01139 numberInRow[iRow] = number;
01140 deleteLink ( iRow );
01141 addLink ( iRow, number );
01142 }
01143 markRow[pivotRow] = -1;
01144
01145 deleteLink ( pivotRow );
01146 deleteLink ( pivotColumn + numberRows_ );
01147 totalElements_ += added;
01148 return true;
01149 }
01150
01151
01153
01154 protected:
01155
01158
01159 double pivotTolerance_;
01161 double zeroTolerance_;
01162 #ifndef COIN_FAST_CODE
01164 double slackValue_;
01165 #else
01166 #ifndef slackValue_
01167 #define slackValue_ -1.0
01168 #endif
01169 #endif
01171 double areaFactor_;
01173 double relaxCheck_;
01175 int numberRows_;
01177 int numberRowsExtra_;
01179 int maximumRowsExtra_;
01181 int numberColumns_;
01183 int numberColumnsExtra_;
01185 int maximumColumnsExtra_;
01187 int numberGoodU_;
01189 int numberGoodL_;
01191 int maximumPivots_;
01193 int numberPivots_;
01196 CoinBigIndex totalElements_;
01198 CoinBigIndex factorElements_;
01200 CoinIntArrayWithLength pivotColumn_;
01202 CoinIntArrayWithLength permute_;
01204 CoinIntArrayWithLength permuteBack_;
01206 CoinIntArrayWithLength pivotColumnBack_;
01208 int status_;
01209
01214
01215
01217 int numberTrials_;
01219 CoinBigIndexArrayWithLength startRowU_;
01220
01222 CoinIntArrayWithLength numberInRow_;
01223
01225 CoinIntArrayWithLength numberInColumn_;
01226
01228 CoinIntArrayWithLength numberInColumnPlus_;
01229
01232 CoinIntArrayWithLength firstCount_;
01233
01235 CoinIntArrayWithLength nextCount_;
01236
01238 CoinIntArrayWithLength lastCount_;
01239
01241 CoinIntArrayWithLength nextColumn_;
01242
01244 CoinIntArrayWithLength lastColumn_;
01245
01247 CoinIntArrayWithLength nextRow_;
01248
01250 CoinIntArrayWithLength lastRow_;
01251
01253 CoinIntArrayWithLength saveColumn_;
01254
01256 CoinIntArrayWithLength markRow_;
01257
01259 int messageLevel_;
01260
01262 int biggerDimension_;
01263
01265 CoinIntArrayWithLength indexColumnU_;
01266
01268 CoinIntArrayWithLength pivotRowL_;
01269
01271 CoinDoubleArrayWithLength pivotRegion_;
01272
01274 int numberSlacks_;
01275
01277 int numberU_;
01278
01280 CoinBigIndex maximumU_;
01281
01283
01284
01286 CoinBigIndex lengthU_;
01287
01289 CoinBigIndex lengthAreaU_;
01290
01292 CoinDoubleArrayWithLength elementU_;
01293
01295 CoinIntArrayWithLength indexRowU_;
01296
01298 CoinBigIndexArrayWithLength startColumnU_;
01299
01301 CoinBigIndexArrayWithLength convertRowToColumnU_;
01302
01304 CoinBigIndex numberL_;
01305
01307 CoinBigIndex baseL_;
01308
01310 CoinBigIndex lengthL_;
01311
01313 CoinBigIndex lengthAreaL_;
01314
01316 CoinDoubleArrayWithLength elementL_;
01317
01319 CoinIntArrayWithLength indexRowL_;
01320
01322 CoinBigIndexArrayWithLength startColumnL_;
01323
01325 bool doForrestTomlin_;
01326
01328 int numberR_;
01329
01331 CoinBigIndex lengthR_;
01332
01334 CoinBigIndex lengthAreaR_;
01335
01337 double *elementR_;
01338
01340 int *indexRowR_;
01341
01343 CoinBigIndexArrayWithLength startColumnR_;
01344
01346 double * denseArea_;
01347
01349 int * densePermute_;
01350
01352 int numberDense_;
01353
01355 int denseThreshold_;
01356
01358 CoinDoubleArrayWithLength workArea_;
01359
01361 CoinUnsignedIntArrayWithLength workArea2_;
01362
01364 CoinBigIndex numberCompressions_;
01365
01367 mutable double ftranCountInput_;
01368 mutable double ftranCountAfterL_;
01369 mutable double ftranCountAfterR_;
01370 mutable double ftranCountAfterU_;
01371 mutable double btranCountInput_;
01372 mutable double btranCountAfterU_;
01373 mutable double btranCountAfterR_;
01374 mutable double btranCountAfterL_;
01375
01377 mutable int numberFtranCounts_;
01378 mutable int numberBtranCounts_;
01379
01381 double ftranAverageAfterL_;
01382 double ftranAverageAfterR_;
01383 double ftranAverageAfterU_;
01384 double btranAverageAfterU_;
01385 double btranAverageAfterR_;
01386 double btranAverageAfterL_;
01387
01389 mutable bool collectStatistics_;
01390
01392 int sparseThreshold_;
01393
01395 int sparseThreshold2_;
01396
01398 CoinBigIndexArrayWithLength startRowL_;
01399
01401 CoinIntArrayWithLength indexColumnL_;
01402
01404 CoinDoubleArrayWithLength elementByRowL_;
01405
01407 mutable CoinIntArrayWithLength sparse_;
01411 int biasLU_;
01417 int persistenceFlag_;
01419 };
01420
01421 #ifdef COIN_HAS_LAPACK
01422 #define DENSE_CODE 1
01423
01424 #ifndef ipfint
01425
01426 typedef int ipfint;
01427 typedef const int cipfint;
01428 #endif
01429 #endif
01430 #endif