00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #ifndef CoinFactorization_H
00012 #define CoinFactorization_H
00013
00014
00015 #include <iostream>
00016 #include <string>
00017 #include <cassert>
00018 #include <cstdio>
00019 #include "CoinFinite.hpp"
00020 #include "CoinIndexedVector.hpp"
00021 class CoinPackedMatrix;
00047 class CoinFactorization {
00048 friend void CoinFactorizationUnitTest( const std::string & mpsDir );
00049
00050 public:
00051
00054
00055 CoinFactorization ( );
00057 CoinFactorization ( const CoinFactorization &other);
00058
00060 ~CoinFactorization ( );
00062 void almostDestructor();
00064 void show_self ( ) const;
00066 int saveFactorization (const char * file ) const;
00070 int restoreFactorization (const char * file , bool factor=false) ;
00072 void sort ( ) const;
00074 CoinFactorization & operator = ( const CoinFactorization & other );
00076
00086 int factorize ( const CoinPackedMatrix & matrix,
00087 int rowIsBasic[], int columnIsBasic[] ,
00088 double areaFactor = 0.0 );
00099 int factorize ( int numberRows,
00100 int numberColumns,
00101 CoinBigIndex numberElements,
00102 CoinBigIndex maximumL,
00103 CoinBigIndex maximumU,
00104 const int indicesRow[],
00105 const int indicesColumn[], const double elements[] ,
00106 int permutation[],
00107 double areaFactor = 0.0);
00112 int factorizePart1 ( int numberRows,
00113 int numberColumns,
00114 CoinBigIndex estimateNumberElements,
00115 int * indicesRow[],
00116 int * indicesColumn[],
00117 CoinFactorizationDouble * elements[],
00118 double areaFactor = 0.0);
00125 int factorizePart2 (int permutation[],int exactNumberElements);
00127 double conditionNumber() const;
00128
00130
00133
00134 inline int status ( ) const {
00135 return status_;
00136 }
00138 inline void setStatus ( int value)
00139 { status_=value; }
00141 inline int pivots ( ) const {
00142 return numberPivots_;
00143 }
00145 inline void setPivots ( int value )
00146 { numberPivots_=value; }
00148 inline int *permute ( ) const {
00149 return permute_.array();
00150 }
00152 inline int *pivotColumn ( ) const {
00153 return pivotColumn_.array();
00154 }
00156 inline CoinFactorizationDouble *pivotRegion ( ) const {
00157 return pivotRegion_.array();
00158 }
00160 inline int *permuteBack ( ) const {
00161 return permuteBack_.array();
00162 }
00165 inline int *pivotColumnBack ( ) const {
00166
00167 return pivotColumnBack_.array();
00168 }
00170 inline CoinBigIndex * startRowL() const
00171 { return startRowL_.array();}
00172
00174 inline CoinBigIndex * startColumnL() const
00175 { return startColumnL_.array();}
00176
00178 inline int * indexColumnL() const
00179 { return indexColumnL_.array();}
00180
00182 inline int * indexRowL() const
00183 { return indexRowL_.array();}
00184
00186 inline CoinFactorizationDouble * elementByRowL() const
00187 { return elementByRowL_.array();}
00188
00190 inline int numberRowsExtra ( ) const {
00191 return numberRowsExtra_;
00192 }
00194 inline void setNumberRows(int value)
00195 { numberRows_ = value; }
00197 inline int numberRows ( ) const {
00198 return numberRows_;
00199 }
00201 inline CoinBigIndex numberL() const
00202 { return numberL_;}
00203
00205 inline CoinBigIndex baseL() const
00206 { return baseL_;}
00208 inline int maximumRowsExtra ( ) const {
00209 return maximumRowsExtra_;
00210 }
00212 inline int numberColumns ( ) const {
00213 return numberColumns_;
00214 }
00216 inline int numberElements ( ) const {
00217 return totalElements_;
00218 }
00220 inline int numberForrestTomlin ( ) const {
00221 return numberInColumn_.array()[numberColumnsExtra_];
00222 }
00224 inline int numberGoodColumns ( ) const {
00225 return numberGoodU_;
00226 }
00228 inline double areaFactor ( ) const {
00229 return areaFactor_;
00230 }
00231 inline void areaFactor ( double value ) {
00232 areaFactor_=value;
00233 }
00235 double adjustedAreaFactor() const;
00237 inline void relaxAccuracyCheck(double value)
00238 { relaxCheck_ = value;}
00239 inline double getAccuracyCheck() const
00240 { return relaxCheck_;}
00242 inline int messageLevel ( ) const {
00243 return messageLevel_ ;
00244 }
00245 void messageLevel ( int value );
00247 inline int maximumPivots ( ) const {
00248 return maximumPivots_ ;
00249 }
00250 void maximumPivots ( int value );
00251
00253 inline int denseThreshold() const
00254 { return denseThreshold_;}
00256 inline void setDenseThreshold(int value)
00257 { denseThreshold_ = value;}
00259 inline double pivotTolerance ( ) const {
00260 return pivotTolerance_ ;
00261 }
00262 void pivotTolerance ( double value );
00264 inline double zeroTolerance ( ) const {
00265 return zeroTolerance_ ;
00266 }
00267 void zeroTolerance ( double value );
00268 #ifndef COIN_FAST_CODE
00270 inline double slackValue ( ) const {
00271 return slackValue_ ;
00272 }
00273 void slackValue ( double value );
00274 #endif
00276 double maximumCoefficient() const;
00278 inline bool forrestTomlin() const
00279 { return doForrestTomlin_;}
00280 inline void setForrestTomlin(bool value)
00281 { doForrestTomlin_=value;}
00283 inline bool spaceForForrestTomlin() const
00284 {
00285 CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
00286 CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
00287 return (space>=0)&&doForrestTomlin_;
00288 }
00290
00293
00295 inline int numberDense() const
00296 { return numberDense_;}
00297
00299 inline CoinBigIndex numberElementsU ( ) const {
00300 return lengthU_;
00301 }
00303 inline void setNumberElementsU(CoinBigIndex value)
00304 { lengthU_ = value; }
00306 inline CoinBigIndex lengthAreaU ( ) const {
00307 return lengthAreaU_;
00308 }
00310 inline CoinBigIndex numberElementsL ( ) const {
00311 return lengthL_;
00312 }
00314 inline CoinBigIndex lengthAreaL ( ) const {
00315 return lengthAreaL_;
00316 }
00318 inline CoinBigIndex numberElementsR ( ) const {
00319 return lengthR_;
00320 }
00322 inline CoinBigIndex numberCompressions() const
00323 { return numberCompressions_;}
00325 inline int * numberInRow() const
00326 { return numberInRow_.array();}
00328 inline int * numberInColumn() const
00329 { return numberInColumn_.array();}
00331 inline CoinFactorizationDouble * elementU() const
00332 { return elementU_.array();}
00334 inline int * indexRowU() const
00335 { return indexRowU_.array();}
00337 inline CoinBigIndex * startColumnU() const
00338 { return startColumnU_.array();}
00340 inline int maximumColumnsExtra()
00341 { return maximumColumnsExtra_;}
00345 inline int biasLU() const
00346 { return biasLU_;}
00347 inline void setBiasLU(int value)
00348 { biasLU_=value;}
00354 inline int persistenceFlag() const
00355 { return persistenceFlag_;}
00356 void setPersistenceFlag(int value);
00358
00361
00369 int replaceColumn ( CoinIndexedVector * regionSparse,
00370 int pivotRow,
00371 double pivotCheck ,
00372 bool checkBeforeModifying=false,
00373 double acceptablePivot=1.0e-8);
00378 void replaceColumnU ( CoinIndexedVector * regionSparse,
00379 CoinBigIndex * deleted,
00380 int internalPivotRow);
00382
00392 int updateColumnFT ( CoinIndexedVector * regionSparse,
00393 CoinIndexedVector * regionSparse2);
00396 int updateColumn ( CoinIndexedVector * regionSparse,
00397 CoinIndexedVector * regionSparse2,
00398 bool noPermute=false) const;
00404 int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
00405 CoinIndexedVector * regionSparse2,
00406 CoinIndexedVector * regionSparse3,
00407 bool noPermuteRegion3=false) ;
00412 int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00413 CoinIndexedVector * regionSparse2) const;
00415 void goSparse();
00417 inline int sparseThreshold ( ) const
00418 { return sparseThreshold_;}
00420 void sparseThreshold ( int value );
00422
00423
00427
00428 inline void clearArrays()
00429 { gutsOfDestructor();}
00431
00434
00437 int add ( CoinBigIndex numberElements,
00438 int indicesRow[],
00439 int indicesColumn[], double elements[] );
00440
00443 int addColumn ( CoinBigIndex numberElements,
00444 int indicesRow[], double elements[] );
00445
00448 int addRow ( CoinBigIndex numberElements,
00449 int indicesColumn[], double elements[] );
00450
00452 int deleteColumn ( int Row );
00454 int deleteRow ( int Row );
00455
00459 int replaceRow ( int whichRow, int numberElements,
00460 const int indicesColumn[], const double elements[] );
00462 void emptyRows(int numberToEmpty, const int which[]);
00464
00465
00466 void checkSparse();
00468 inline bool collectStatistics() const
00469 { return collectStatistics_;}
00471 inline void setCollectStatistics(bool onOff) const
00472 { collectStatistics_ = onOff;}
00474 void gutsOfDestructor(int type=1);
00476 void gutsOfInitialize(int type);
00477 void gutsOfCopy(const CoinFactorization &other);
00478
00480 void resetStatistics();
00481
00482
00484
00486
00487 void getAreas ( int numberRows,
00488 int numberColumns,
00489 CoinBigIndex maximumL,
00490 CoinBigIndex maximumU );
00491
00494 void preProcess ( int state,
00495 int possibleDuplicates = -1 );
00497 int factor ( );
00498 protected:
00501 int factorSparse ( );
00504 int factorSparseSmall ( );
00507 int factorSparseLarge ( );
00510 int factorDense ( );
00511
00513 bool pivotOneOtherRow ( int pivotRow,
00514 int pivotColumn );
00516 bool pivotRowSingleton ( int pivotRow,
00517 int pivotColumn );
00519 bool pivotColumnSingleton ( int pivotRow,
00520 int pivotColumn );
00521
00526 bool getColumnSpace ( int iColumn,
00527 int extraNeeded );
00528
00531 bool reorderU();
00535 bool getColumnSpaceIterateR ( int iColumn, double value,
00536 int iRow);
00542 CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00543 int iRow);
00547 bool getRowSpace ( int iRow, int extraNeeded );
00548
00552 bool getRowSpaceIterate ( int iRow,
00553 int extraNeeded );
00555 void checkConsistency ( );
00557 inline void addLink ( int index, int count ) {
00558 int *nextCount = nextCount_.array();
00559 int *firstCount = firstCount_.array();
00560 int *lastCount = lastCount_.array();
00561 int next = firstCount[count];
00562 lastCount[index] = -2 - count;
00563 if ( next < 0 ) {
00564
00565 firstCount[count] = index;
00566 nextCount[index] = -1;
00567 } else {
00568 firstCount[count] = index;
00569 nextCount[index] = next;
00570 lastCount[next] = index;
00571 }}
00573 inline void deleteLink ( int index ) {
00574 int *nextCount = nextCount_.array();
00575 int *firstCount = firstCount_.array();
00576 int *lastCount = lastCount_.array();
00577 int next = nextCount[index];
00578 int last = lastCount[index];
00579 if ( last >= 0 ) {
00580 nextCount[last] = next;
00581 } else {
00582 int count = -last - 2;
00583
00584 firstCount[count] = next;
00585 }
00586 if ( next >= 0 ) {
00587 lastCount[next] = last;
00588 }
00589 nextCount[index] = -2;
00590 lastCount[index] = -2;
00591 return;
00592 }
00594 void separateLinks(int count,bool rowsFirst);
00596 void cleanup ( );
00597
00599 void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00601 void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00603 void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00605 void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00606
00608 void updateColumnR ( CoinIndexedVector * region ) const;
00611 void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00612
00614 void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00615
00617 void updateColumnUSparse ( CoinIndexedVector * regionSparse,
00618 int * indexIn) const;
00620 void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
00621 int * indexIn) const;
00623 int updateColumnUDensish ( double * COIN_RESTRICT region,
00624 int * COIN_RESTRICT regionIndex) const;
00626 void updateTwoColumnsUDensish (
00627 int & numberNonZero1,
00628 double * COIN_RESTRICT region1,
00629 int * COIN_RESTRICT index1,
00630 int & numberNonZero2,
00631 double * COIN_RESTRICT region2,
00632 int * COIN_RESTRICT index2) const;
00634 void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
00636 void permuteBack ( CoinIndexedVector * regionSparse,
00637 CoinIndexedVector * outVector) const;
00638
00640 void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00643 void updateColumnTransposeU ( CoinIndexedVector * region,
00644 int smallestIndex) const;
00647 void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00648 int smallestIndex) const;
00651 void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00652 int smallestIndex) const;
00655 void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00658 void updateColumnTransposeUByColumn ( CoinIndexedVector * region,
00659 int smallestIndex) const;
00660
00662 void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00664 void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00666 void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00667
00669 void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00671 void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00673 void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00675 void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00677 void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00678 public:
00683 int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00684 int pivotRow, double alpha);
00685 protected:
00688 int checkPivot(double saveFromU, double oldPivot) const;
00689
00690 #ifdef INT_IS_8
00691 #define COINFACTORIZATION_BITS_PER_INT 64
00692 #define COINFACTORIZATION_SHIFT_PER_INT 6
00693 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00694 #else
00695 #define COINFACTORIZATION_BITS_PER_INT 32
00696 #define COINFACTORIZATION_SHIFT_PER_INT 5
00697 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00698 #endif
00699 template <class T> inline bool
00700 pivot ( int pivotRow,
00701 int pivotColumn,
00702 CoinBigIndex pivotRowPosition,
00703 CoinBigIndex pivotColumnPosition,
00704 CoinFactorizationDouble work[],
00705 unsigned int workArea2[],
00706 int increment2,
00707 T markRow[] ,
00708 int largeInteger)
00709 {
00710 int *indexColumnU = indexColumnU_.array();
00711 CoinBigIndex *startColumnU = startColumnU_.array();
00712 int *numberInColumn = numberInColumn_.array();
00713 CoinFactorizationDouble *elementU = elementU_.array();
00714 int *indexRowU = indexRowU_.array();
00715 CoinBigIndex *startRowU = startRowU_.array();
00716 int *numberInRow = numberInRow_.array();
00717 CoinFactorizationDouble *elementL = elementL_.array();
00718 int *indexRowL = indexRowL_.array();
00719 int *saveColumn = saveColumn_.array();
00720 int *nextRow = nextRow_.array();
00721 int *lastRow = lastRow_.array() ;
00722
00723
00724 int numberInPivotRow = numberInRow[pivotRow] - 1;
00725 CoinBigIndex startColumn = startColumnU[pivotColumn];
00726 int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00727 CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00728 int put = 0;
00729 CoinBigIndex startRow = startRowU[pivotRow];
00730 CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00731
00732 if ( pivotColumnPosition < 0 ) {
00733 for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00734 int iColumn = indexColumnU[pivotColumnPosition];
00735 if ( iColumn != pivotColumn ) {
00736 saveColumn[put++] = iColumn;
00737 } else {
00738 break;
00739 }
00740 }
00741 } else {
00742 for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00743 saveColumn[put++] = indexColumnU[i];
00744 }
00745 }
00746 assert (pivotColumnPosition<endRow);
00747 assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00748 pivotColumnPosition++;
00749 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00750 saveColumn[put++] = indexColumnU[pivotColumnPosition];
00751 }
00752
00753 int next = nextRow[pivotRow];
00754 int last = lastRow[pivotRow];
00755
00756 nextRow[last] = next;
00757 lastRow[next] = last;
00758 nextRow[pivotRow] = numberGoodU_;
00759 lastRow[pivotRow] = -2;
00760 numberInRow[pivotRow] = 0;
00761
00762 CoinBigIndex l = lengthL_;
00763
00764 if ( l + numberInPivotColumn > lengthAreaL_ ) {
00765
00766 if ((messageLevel_&4)!=0)
00767 printf("more memory needed in middle of invert\n");
00768 return false;
00769 }
00770
00771 CoinBigIndex lSave = l;
00772
00773 CoinBigIndex * startColumnL = startColumnL_.array();
00774 startColumnL[numberGoodL_] = l;
00775 numberGoodL_++;
00776 startColumnL[numberGoodL_] = l + numberInPivotColumn;
00777 lengthL_ += numberInPivotColumn;
00778 if ( pivotRowPosition < 0 ) {
00779 for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00780 int iRow = indexRowU[pivotRowPosition];
00781 if ( iRow != pivotRow ) {
00782 indexRowL[l] = iRow;
00783 elementL[l] = elementU[pivotRowPosition];
00784 markRow[iRow] = static_cast<T>(l - lSave);
00785 l++;
00786
00787 CoinBigIndex start = startRowU[iRow];
00788 CoinBigIndex end = start + numberInRow[iRow];
00789 CoinBigIndex where = start;
00790
00791 while ( indexColumnU[where] != pivotColumn ) {
00792 where++;
00793 }
00794 #if DEBUG_COIN
00795 if ( where >= end ) {
00796 abort ( );
00797 }
00798 #endif
00799 indexColumnU[where] = indexColumnU[end - 1];
00800 numberInRow[iRow]--;
00801 } else {
00802 break;
00803 }
00804 }
00805 } else {
00806 CoinBigIndex i;
00807
00808 for ( i = startColumn; i < pivotRowPosition; i++ ) {
00809 int iRow = indexRowU[i];
00810
00811 markRow[iRow] = static_cast<T>(l - lSave);
00812 indexRowL[l] = iRow;
00813 elementL[l] = elementU[i];
00814 l++;
00815
00816 CoinBigIndex start = startRowU[iRow];
00817 CoinBigIndex end = start + numberInRow[iRow];
00818 CoinBigIndex where = start;
00819
00820 while ( indexColumnU[where] != pivotColumn ) {
00821 where++;
00822 }
00823 #if DEBUG_COIN
00824 if ( where >= end ) {
00825 abort ( );
00826 }
00827 #endif
00828 indexColumnU[where] = indexColumnU[end - 1];
00829 numberInRow[iRow]--;
00830 assert (numberInRow[iRow]>=0);
00831 }
00832 }
00833 assert (pivotRowPosition<endColumn);
00834 assert (indexRowU[pivotRowPosition]==pivotRow);
00835 CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
00836 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
00837
00838 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00839 pivotRowPosition++;
00840 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00841 int iRow = indexRowU[pivotRowPosition];
00842
00843 markRow[iRow] = static_cast<T>(l - lSave);
00844 indexRowL[l] = iRow;
00845 elementL[l] = elementU[pivotRowPosition];
00846 l++;
00847
00848 CoinBigIndex start = startRowU[iRow];
00849 CoinBigIndex end = start + numberInRow[iRow];
00850 CoinBigIndex where = start;
00851
00852 while ( indexColumnU[where] != pivotColumn ) {
00853 where++;
00854 }
00855 #if DEBUG_COIN
00856 if ( where >= end ) {
00857 abort ( );
00858 }
00859 #endif
00860 indexColumnU[where] = indexColumnU[end - 1];
00861 numberInRow[iRow]--;
00862 assert (numberInRow[iRow]>=0);
00863 }
00864 markRow[pivotRow] = static_cast<T>(largeInteger);
00865
00866 numberInColumn[pivotColumn] = 0;
00867
00868 int *indexL = &indexRowL[lSave];
00869 CoinFactorizationDouble *multipliersL = &elementL[lSave];
00870
00871
00872 int j;
00873
00874 for ( j = 0; j < numberInPivotColumn; j++ ) {
00875 multipliersL[j] *= pivotMultiplier;
00876 }
00877
00878 CoinBigIndex iErase;
00879 for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00880 iErase++ ) {
00881 workArea2[iErase] = 0;
00882 }
00883 CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00884 unsigned int *temp2 = workArea2;
00885 int * nextColumn = nextColumn_.array();
00886
00887
00888 int jColumn;
00889 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
00890 int iColumn = saveColumn[jColumn];
00891 CoinBigIndex startColumn = startColumnU[iColumn];
00892 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
00893 int iRow = indexRowU[startColumn];
00894 CoinFactorizationDouble value = elementU[startColumn];
00895 double largest;
00896 CoinBigIndex put = startColumn;
00897 CoinBigIndex positionLargest = -1;
00898 CoinFactorizationDouble thisPivotValue = 0.0;
00899
00900
00901 bool checkLargest;
00902 int mark = markRow[iRow];
00903
00904 if ( mark == largeInteger+1 ) {
00905 largest = fabs ( value );
00906 positionLargest = put;
00907 put++;
00908 checkLargest = false;
00909 } else {
00910
00911 largest = 0.0;
00912 checkLargest = true;
00913 if ( mark != largeInteger ) {
00914
00915 work[mark] = value;
00916 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00917 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00918
00919 temp2[word] = temp2[word] | ( 1 << bit );
00920 added--;
00921 } else {
00922 thisPivotValue = value;
00923 }
00924 }
00925 CoinBigIndex i;
00926 for ( i = startColumn + 1; i < endColumn; i++ ) {
00927 iRow = indexRowU[i];
00928 value = elementU[i];
00929 int mark = markRow[iRow];
00930
00931 if ( mark == largeInteger+1 ) {
00932
00933 indexRowU[put] = iRow;
00934 elementU[put] = value;
00935 if ( checkLargest ) {
00936 double absValue = fabs ( value );
00937
00938 if ( absValue > largest ) {
00939 largest = absValue;
00940 positionLargest = put;
00941 }
00942 }
00943 put++;
00944 } else if ( mark != largeInteger ) {
00945
00946 work[mark] = value;
00947 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00948 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00949
00950 temp2[word] = temp2[word] | ( 1 << bit );
00951 added--;
00952 } else {
00953 thisPivotValue = value;
00954 }
00955 }
00956
00957 elementU[put] = elementU[startColumn];
00958 indexRowU[put] = indexRowU[startColumn];
00959 if ( positionLargest == startColumn ) {
00960 positionLargest = put;
00961 }
00962 put++;
00963 elementU[startColumn] = thisPivotValue;
00964 indexRowU[startColumn] = pivotRow;
00965
00966 startColumn++;
00967 numberInColumn[iColumn] = put - startColumn;
00968 int * numberInColumnPlus = numberInColumnPlus_.array();
00969 numberInColumnPlus[iColumn]++;
00970 startColumnU[iColumn]++;
00971
00972 int next = nextColumn[iColumn];
00973 CoinBigIndex space;
00974
00975 space = startColumnU[next] - put - numberInColumnPlus[next];
00976
00977 if ( numberInPivotColumn > space ) {
00978
00979 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
00980 return false;
00981 }
00982
00983 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
00984 startColumn = startColumnU[iColumn];
00985 put = startColumn + numberInColumn[iColumn];
00986 }
00987 double tolerance = zeroTolerance_;
00988
00989 int *nextCount = nextCount_.array();
00990 for ( j = 0; j < numberInPivotColumn; j++ ) {
00991 value = work[j] - thisPivotValue * multipliersL[j];
00992 double absValue = fabs ( value );
00993
00994 if ( absValue > tolerance ) {
00995 work[j] = 0.0;
00996 elementU[put] = value;
00997 indexRowU[put] = indexL[j];
00998 if ( absValue > largest ) {
00999 largest = absValue;
01000 positionLargest = put;
01001 }
01002 put++;
01003 } else {
01004 work[j] = 0.0;
01005 added--;
01006 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01007 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01008
01009 if ( temp2[word] & ( 1 << bit ) ) {
01010
01011 iRow = indexL[j];
01012 CoinBigIndex start = startRowU[iRow];
01013 CoinBigIndex end = start + numberInRow[iRow];
01014 CoinBigIndex where = start;
01015
01016 while ( indexColumnU[where] != iColumn ) {
01017 where++;
01018 }
01019 #if DEBUG_COIN
01020 if ( where >= end ) {
01021 abort ( );
01022 }
01023 #endif
01024 indexColumnU[where] = indexColumnU[end - 1];
01025 numberInRow[iRow]--;
01026 } else {
01027
01028 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01029 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01030
01031 temp2[word] = temp2[word] | ( 1 << bit );
01032 }
01033 }
01034 }
01035 numberInColumn[iColumn] = put - startColumn;
01036
01037 if ( positionLargest >= 0 ) {
01038 value = elementU[positionLargest];
01039 iRow = indexRowU[positionLargest];
01040 elementU[positionLargest] = elementU[startColumn];
01041 indexRowU[positionLargest] = indexRowU[startColumn];
01042 elementU[startColumn] = value;
01043 indexRowU[startColumn] = iRow;
01044 }
01045
01046 if ( nextCount[iColumn + numberRows_] != -2 ) {
01047
01048 deleteLink ( iColumn + numberRows_ );
01049 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01050 }
01051 temp2 += increment2;
01052 }
01053
01054 unsigned int *putBase = workArea2;
01055 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01056 int i = 0;
01057
01058
01059 while ( bigLoops ) {
01060 bigLoops--;
01061 int bit;
01062 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01063 unsigned int *putThis = putBase;
01064 int iRow = indexL[i];
01065
01066
01067 int number = 0;
01068 int jColumn;
01069
01070 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01071 unsigned int test = *putThis;
01072
01073 putThis += increment2;
01074 test = 1 - ( ( test >> bit ) & 1 );
01075 number += test;
01076 }
01077 int next = nextRow[iRow];
01078 CoinBigIndex space;
01079
01080 space = startRowU[next] - startRowU[iRow];
01081 number += numberInRow[iRow];
01082 if ( space < number ) {
01083 if ( !getRowSpace ( iRow, number ) ) {
01084 return false;
01085 }
01086 }
01087
01088 putThis = putBase;
01089 next = nextRow[iRow];
01090 number = numberInRow[iRow];
01091 CoinBigIndex end = startRowU[iRow] + number;
01092 int saveIndex = indexColumnU[startRowU[next]];
01093
01094
01095 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01096 unsigned int test = *putThis;
01097
01098 putThis += increment2;
01099 test = 1 - ( ( test >> bit ) & 1 );
01100 indexColumnU[end] = saveColumn[jColumn];
01101 end += test;
01102 }
01103
01104 indexColumnU[startRowU[next]] = saveIndex;
01105 markRow[iRow] = static_cast<T>(largeInteger+1);
01106 number = end - startRowU[iRow];
01107 numberInRow[iRow] = number;
01108 deleteLink ( iRow );
01109 addLink ( iRow, number );
01110 }
01111 putBase++;
01112 }
01113 int bit;
01114
01115 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01116 unsigned int *putThis = putBase;
01117 int iRow = indexL[i];
01118
01119
01120 int number = 0;
01121 int jColumn;
01122
01123 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01124 unsigned int test = *putThis;
01125
01126 putThis += increment2;
01127 test = 1 - ( ( test >> bit ) & 1 );
01128 number += test;
01129 }
01130 int next = nextRow[iRow];
01131 CoinBigIndex space;
01132
01133 space = startRowU[next] - startRowU[iRow];
01134 number += numberInRow[iRow];
01135 if ( space < number ) {
01136 if ( !getRowSpace ( iRow, number ) ) {
01137 return false;
01138 }
01139 }
01140
01141 putThis = putBase;
01142 next = nextRow[iRow];
01143 number = numberInRow[iRow];
01144 CoinBigIndex end = startRowU[iRow] + number;
01145 int saveIndex;
01146
01147 saveIndex = indexColumnU[startRowU[next]];
01148
01149
01150 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01151 unsigned int test = *putThis;
01152
01153 putThis += increment2;
01154 test = 1 - ( ( test >> bit ) & 1 );
01155
01156 indexColumnU[end] = saveColumn[jColumn];
01157 end += test;
01158 }
01159 indexColumnU[startRowU[next]] = saveIndex;
01160 markRow[iRow] = static_cast<T>(largeInteger+1);
01161 number = end - startRowU[iRow];
01162 numberInRow[iRow] = number;
01163 deleteLink ( iRow );
01164 addLink ( iRow, number );
01165 }
01166 markRow[pivotRow] = static_cast<T>(largeInteger+1);
01167
01168 deleteLink ( pivotRow );
01169 deleteLink ( pivotColumn + numberRows_ );
01170 totalElements_ += added;
01171 return true;
01172 }
01173
01174
01176
01177 protected:
01178
01181
01182 double pivotTolerance_;
01184 double zeroTolerance_;
01185 #ifndef COIN_FAST_CODE
01187 double slackValue_;
01188 #else
01189 #ifndef slackValue_
01190 #define slackValue_ -1.0
01191 #endif
01192 #endif
01194 double areaFactor_;
01196 double relaxCheck_;
01198 int numberRows_;
01200 int numberRowsExtra_;
01202 int maximumRowsExtra_;
01204 int numberColumns_;
01206 int numberColumnsExtra_;
01208 int maximumColumnsExtra_;
01210 int numberGoodU_;
01212 int numberGoodL_;
01214 int maximumPivots_;
01216 int numberPivots_;
01219 CoinBigIndex totalElements_;
01221 CoinBigIndex factorElements_;
01223 CoinIntArrayWithLength pivotColumn_;
01225 CoinIntArrayWithLength permute_;
01227 CoinIntArrayWithLength permuteBack_;
01229 CoinIntArrayWithLength pivotColumnBack_;
01231 int status_;
01232
01237
01238
01240 int numberTrials_;
01242 CoinBigIndexArrayWithLength startRowU_;
01243
01245 CoinIntArrayWithLength numberInRow_;
01246
01248 CoinIntArrayWithLength numberInColumn_;
01249
01251 CoinIntArrayWithLength numberInColumnPlus_;
01252
01255 CoinIntArrayWithLength firstCount_;
01256
01258 CoinIntArrayWithLength nextCount_;
01259
01261 CoinIntArrayWithLength lastCount_;
01262
01264 CoinIntArrayWithLength nextColumn_;
01265
01267 CoinIntArrayWithLength lastColumn_;
01268
01270 CoinIntArrayWithLength nextRow_;
01271
01273 CoinIntArrayWithLength lastRow_;
01274
01276 CoinIntArrayWithLength saveColumn_;
01277
01279 CoinIntArrayWithLength markRow_;
01280
01282 int messageLevel_;
01283
01285 int biggerDimension_;
01286
01288 CoinIntArrayWithLength indexColumnU_;
01289
01291 CoinIntArrayWithLength pivotRowL_;
01292
01294 CoinFactorizationDoubleArrayWithLength pivotRegion_;
01295
01297 int numberSlacks_;
01298
01300 int numberU_;
01301
01303 CoinBigIndex maximumU_;
01304
01306
01307
01309 CoinBigIndex lengthU_;
01310
01312 CoinBigIndex lengthAreaU_;
01313
01315 CoinFactorizationDoubleArrayWithLength elementU_;
01316
01318 CoinIntArrayWithLength indexRowU_;
01319
01321 CoinBigIndexArrayWithLength startColumnU_;
01322
01324 CoinBigIndexArrayWithLength convertRowToColumnU_;
01325
01327 CoinBigIndex numberL_;
01328
01330 CoinBigIndex baseL_;
01331
01333 CoinBigIndex lengthL_;
01334
01336 CoinBigIndex lengthAreaL_;
01337
01339 CoinFactorizationDoubleArrayWithLength elementL_;
01340
01342 CoinIntArrayWithLength indexRowL_;
01343
01345 CoinBigIndexArrayWithLength startColumnL_;
01346
01348 bool doForrestTomlin_;
01349
01351 int numberR_;
01352
01354 CoinBigIndex lengthR_;
01355
01357 CoinBigIndex lengthAreaR_;
01358
01360 CoinFactorizationDouble *elementR_;
01361
01363 int *indexRowR_;
01364
01366 CoinBigIndexArrayWithLength startColumnR_;
01367
01369 double * denseArea_;
01370
01372 int * densePermute_;
01373
01375 int numberDense_;
01376
01378 int denseThreshold_;
01379
01381 CoinFactorizationDoubleArrayWithLength workArea_;
01382
01384 CoinUnsignedIntArrayWithLength workArea2_;
01385
01387 CoinBigIndex numberCompressions_;
01388
01390 mutable double ftranCountInput_;
01391 mutable double ftranCountAfterL_;
01392 mutable double ftranCountAfterR_;
01393 mutable double ftranCountAfterU_;
01394 mutable double btranCountInput_;
01395 mutable double btranCountAfterU_;
01396 mutable double btranCountAfterR_;
01397 mutable double btranCountAfterL_;
01398
01400 mutable int numberFtranCounts_;
01401 mutable int numberBtranCounts_;
01402
01404 double ftranAverageAfterL_;
01405 double ftranAverageAfterR_;
01406 double ftranAverageAfterU_;
01407 double btranAverageAfterU_;
01408 double btranAverageAfterR_;
01409 double btranAverageAfterL_;
01410
01412 mutable bool collectStatistics_;
01413
01415 int sparseThreshold_;
01416
01418 int sparseThreshold2_;
01419
01421 CoinBigIndexArrayWithLength startRowL_;
01422
01424 CoinIntArrayWithLength indexColumnL_;
01425
01427 CoinFactorizationDoubleArrayWithLength elementByRowL_;
01428
01430 mutable CoinIntArrayWithLength sparse_;
01434 int biasLU_;
01440 int persistenceFlag_;
01442 };
01443
01444 #ifdef COIN_HAS_LAPACK
01445 #define DENSE_CODE 1
01446
01447 #ifndef ipfint
01448
01449 typedef int ipfint;
01450 typedef const int cipfint;
01451 #endif
01452 #endif
01453 #endif
01454
01455 #ifdef UGLY_COIN_FACTOR_CODING
01456 #define FAC_UNSET (FAC_SET+1)
01457 {
01458 goodPivot=false;
01459
01460 CoinBigIndex startColumnThis = startColumn[iPivotColumn];
01461 CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
01462 int put = 0;
01463 CoinBigIndex startRowThis = startRow[iPivotRow];
01464 CoinBigIndex endRow = startRowThis + numberDoRow + 1;
01465 if ( pivotColumnPosition < 0 ) {
01466 for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01467 int iColumn = indexColumn[pivotColumnPosition];
01468 if ( iColumn != iPivotColumn ) {
01469 saveColumn[put++] = iColumn;
01470 } else {
01471 break;
01472 }
01473 }
01474 } else {
01475 for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
01476 saveColumn[put++] = indexColumn[i];
01477 }
01478 }
01479 assert (pivotColumnPosition<endRow);
01480 assert (indexColumn[pivotColumnPosition]==iPivotColumn);
01481 pivotColumnPosition++;
01482 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01483 saveColumn[put++] = indexColumn[pivotColumnPosition];
01484 }
01485
01486 int next = nextRow[iPivotRow];
01487 int last = lastRow[iPivotRow];
01488
01489 nextRow[last] = next;
01490 lastRow[next] = last;
01491 nextRow[iPivotRow] = numberGoodU_;
01492 lastRow[iPivotRow] = -2;
01493 numberInRow[iPivotRow] = 0;
01494
01495 CoinBigIndex l = lengthL_;
01496
01497 {
01498 if ( l + numberDoColumn > lengthAreaL_ ) {
01499
01500 if ((messageLevel_&4)!=0)
01501 printf("more memory needed in middle of invert\n");
01502 goto BAD_PIVOT;
01503 }
01504
01505 CoinBigIndex lSave = l;
01506
01507 CoinBigIndex * startColumnL = startColumnL_.array();
01508 startColumnL[numberGoodL_] = l;
01509 numberGoodL_++;
01510 startColumnL[numberGoodL_] = l + numberDoColumn;
01511 lengthL_ += numberDoColumn;
01512 if ( pivotRowPosition < 0 ) {
01513 for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01514 int iRow = indexRow[pivotRowPosition];
01515 if ( iRow != iPivotRow ) {
01516 indexRowL[l] = iRow;
01517 elementL[l] = element[pivotRowPosition];
01518 markRow[iRow] = l - lSave;
01519 l++;
01520
01521 CoinBigIndex start = startRow[iRow];
01522 CoinBigIndex end = start + numberInRow[iRow];
01523 CoinBigIndex where = start;
01524
01525 while ( indexColumn[where] != iPivotColumn ) {
01526 where++;
01527 }
01528 #if DEBUG_COIN
01529 if ( where >= end ) {
01530 abort ( );
01531 }
01532 #endif
01533 indexColumn[where] = indexColumn[end - 1];
01534 numberInRow[iRow]--;
01535 } else {
01536 break;
01537 }
01538 }
01539 } else {
01540 CoinBigIndex i;
01541
01542 for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
01543 int iRow = indexRow[i];
01544
01545 markRow[iRow] = l - lSave;
01546 indexRowL[l] = iRow;
01547 elementL[l] = element[i];
01548 l++;
01549
01550 CoinBigIndex start = startRow[iRow];
01551 CoinBigIndex end = start + numberInRow[iRow];
01552 CoinBigIndex where = start;
01553
01554 while ( indexColumn[where] != iPivotColumn ) {
01555 where++;
01556 }
01557 #if DEBUG_COIN
01558 if ( where >= end ) {
01559 abort ( );
01560 }
01561 #endif
01562 indexColumn[where] = indexColumn[end - 1];
01563 numberInRow[iRow]--;
01564 assert (numberInRow[iRow]>=0);
01565 }
01566 }
01567 assert (pivotRowPosition<endColumn);
01568 assert (indexRow[pivotRowPosition]==iPivotRow);
01569 CoinFactorizationDouble pivotElement = element[pivotRowPosition];
01570 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
01571
01572 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
01573 pivotRowPosition++;
01574 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01575 int iRow = indexRow[pivotRowPosition];
01576
01577 markRow[iRow] = l - lSave;
01578 indexRowL[l] = iRow;
01579 elementL[l] = element[pivotRowPosition];
01580 l++;
01581
01582 CoinBigIndex start = startRow[iRow];
01583 CoinBigIndex end = start + numberInRow[iRow];
01584 CoinBigIndex where = start;
01585
01586 while ( indexColumn[where] != iPivotColumn ) {
01587 where++;
01588 }
01589 #if DEBUG_COIN
01590 if ( where >= end ) {
01591 abort ( );
01592 }
01593 #endif
01594 indexColumn[where] = indexColumn[end - 1];
01595 numberInRow[iRow]--;
01596 assert (numberInRow[iRow]>=0);
01597 }
01598 markRow[iPivotRow] = FAC_SET;
01599
01600 numberInColumn[iPivotColumn] = 0;
01601
01602 int *indexL = &indexRowL[lSave];
01603 CoinFactorizationDouble *multipliersL = &elementL[lSave];
01604
01605
01606 int j;
01607
01608 for ( j = 0; j < numberDoColumn; j++ ) {
01609 multipliersL[j] *= pivotMultiplier;
01610 }
01611
01612 CoinBigIndex iErase;
01613 for ( iErase = 0; iErase < increment2 * numberDoRow;
01614 iErase++ ) {
01615 workArea2[iErase] = 0;
01616 }
01617 CoinBigIndex added = numberDoRow * numberDoColumn;
01618 unsigned int *temp2 = workArea2;
01619 int * nextColumn = nextColumn_.array();
01620
01621
01622 int jColumn;
01623 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01624 int iColumn = saveColumn[jColumn];
01625 CoinBigIndex startColumnThis = startColumn[iColumn];
01626 CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
01627 int iRow = indexRow[startColumnThis];
01628 CoinFactorizationDouble value = element[startColumnThis];
01629 double largest;
01630 CoinBigIndex put = startColumnThis;
01631 CoinBigIndex positionLargest = -1;
01632 CoinFactorizationDouble thisPivotValue = 0.0;
01633
01634
01635 bool checkLargest;
01636 int mark = markRow[iRow];
01637
01638 if ( mark == FAC_UNSET ) {
01639 largest = fabs ( value );
01640 positionLargest = put;
01641 put++;
01642 checkLargest = false;
01643 } else {
01644
01645 largest = 0.0;
01646 checkLargest = true;
01647 if ( mark != FAC_SET ) {
01648
01649 workArea[mark] = value;
01650 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01651 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01652
01653 temp2[word] = temp2[word] | ( 1 << bit );
01654 added--;
01655 } else {
01656 thisPivotValue = value;
01657 }
01658 }
01659 CoinBigIndex i;
01660 for ( i = startColumnThis + 1; i < endColumn; i++ ) {
01661 iRow = indexRow[i];
01662 value = element[i];
01663 int mark = markRow[iRow];
01664
01665 if ( mark == FAC_UNSET ) {
01666
01667 indexRow[put] = iRow;
01668 element[put] = value;
01669 if ( checkLargest ) {
01670 double absValue = fabs ( value );
01671
01672 if ( absValue > largest ) {
01673 largest = absValue;
01674 positionLargest = put;
01675 }
01676 }
01677 put++;
01678 } else if ( mark != FAC_SET ) {
01679
01680 workArea[mark] = value;
01681 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01682 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01683
01684 temp2[word] = temp2[word] | ( 1 << bit );
01685 added--;
01686 } else {
01687 thisPivotValue = value;
01688 }
01689 }
01690
01691 element[put] = element[startColumnThis];
01692 indexRow[put] = indexRow[startColumnThis];
01693 if ( positionLargest == startColumnThis ) {
01694 positionLargest = put;
01695 }
01696 put++;
01697 element[startColumnThis] = thisPivotValue;
01698 indexRow[startColumnThis] = iPivotRow;
01699
01700 startColumnThis++;
01701 numberInColumn[iColumn] = put - startColumnThis;
01702 int * numberInColumnPlus = numberInColumnPlus_.array();
01703 numberInColumnPlus[iColumn]++;
01704 startColumn[iColumn]++;
01705
01706 int next = nextColumn[iColumn];
01707 CoinBigIndex space;
01708
01709 space = startColumn[next] - put - numberInColumnPlus[next];
01710
01711 if ( numberDoColumn > space ) {
01712
01713 if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
01714 goto BAD_PIVOT;
01715 }
01716
01717 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
01718 startColumnThis = startColumn[iColumn];
01719 put = startColumnThis + numberInColumn[iColumn];
01720 }
01721 double tolerance = zeroTolerance_;
01722
01723 int *nextCount = nextCount_.array();
01724 for ( j = 0; j < numberDoColumn; j++ ) {
01725 value = workArea[j] - thisPivotValue * multipliersL[j];
01726 double absValue = fabs ( value );
01727
01728 if ( absValue > tolerance ) {
01729 workArea[j] = 0.0;
01730 element[put] = value;
01731 indexRow[put] = indexL[j];
01732 if ( absValue > largest ) {
01733 largest = absValue;
01734 positionLargest = put;
01735 }
01736 put++;
01737 } else {
01738 workArea[j] = 0.0;
01739 added--;
01740 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01741 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01742
01743 if ( temp2[word] & ( 1 << bit ) ) {
01744
01745 iRow = indexL[j];
01746 CoinBigIndex start = startRow[iRow];
01747 CoinBigIndex end = start + numberInRow[iRow];
01748 CoinBigIndex where = start;
01749
01750 while ( indexColumn[where] != iColumn ) {
01751 where++;
01752 }
01753 #if DEBUG_COIN
01754 if ( where >= end ) {
01755 abort ( );
01756 }
01757 #endif
01758 indexColumn[where] = indexColumn[end - 1];
01759 numberInRow[iRow]--;
01760 } else {
01761
01762 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01763 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01764
01765 temp2[word] = temp2[word] | ( 1 << bit );
01766 }
01767 }
01768 }
01769 numberInColumn[iColumn] = put - startColumnThis;
01770
01771 if ( positionLargest >= 0 ) {
01772 value = element[positionLargest];
01773 iRow = indexRow[positionLargest];
01774 element[positionLargest] = element[startColumnThis];
01775 indexRow[positionLargest] = indexRow[startColumnThis];
01776 element[startColumnThis] = value;
01777 indexRow[startColumnThis] = iRow;
01778 }
01779
01780 if ( nextCount[iColumn + numberRows_] != -2 ) {
01781
01782 deleteLink ( iColumn + numberRows_ );
01783 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01784 }
01785 temp2 += increment2;
01786 }
01787
01788 unsigned int *putBase = workArea2;
01789 int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01790 int i = 0;
01791
01792
01793 while ( bigLoops ) {
01794 bigLoops--;
01795 int bit;
01796 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01797 unsigned int *putThis = putBase;
01798 int iRow = indexL[i];
01799
01800
01801 int number = 0;
01802 int jColumn;
01803
01804 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01805 unsigned int test = *putThis;
01806
01807 putThis += increment2;
01808 test = 1 - ( ( test >> bit ) & 1 );
01809 number += test;
01810 }
01811 int next = nextRow[iRow];
01812 CoinBigIndex space;
01813
01814 space = startRow[next] - startRow[iRow];
01815 number += numberInRow[iRow];
01816 if ( space < number ) {
01817 if ( !getRowSpace ( iRow, number ) ) {
01818 goto BAD_PIVOT;
01819 }
01820 }
01821
01822 putThis = putBase;
01823 next = nextRow[iRow];
01824 number = numberInRow[iRow];
01825 CoinBigIndex end = startRow[iRow] + number;
01826 int saveIndex = indexColumn[startRow[next]];
01827
01828
01829 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01830 unsigned int test = *putThis;
01831
01832 putThis += increment2;
01833 test = 1 - ( ( test >> bit ) & 1 );
01834 indexColumn[end] = saveColumn[jColumn];
01835 end += test;
01836 }
01837
01838 indexColumn[startRow[next]] = saveIndex;
01839 markRow[iRow] = FAC_UNSET;
01840 number = end - startRow[iRow];
01841 numberInRow[iRow] = number;
01842 deleteLink ( iRow );
01843 addLink ( iRow, number );
01844 }
01845 putBase++;
01846 }
01847 int bit;
01848
01849 for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
01850 unsigned int *putThis = putBase;
01851 int iRow = indexL[i];
01852
01853
01854 int number = 0;
01855 int jColumn;
01856
01857 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01858 unsigned int test = *putThis;
01859
01860 putThis += increment2;
01861 test = 1 - ( ( test >> bit ) & 1 );
01862 number += test;
01863 }
01864 int next = nextRow[iRow];
01865 CoinBigIndex space;
01866
01867 space = startRow[next] - startRow[iRow];
01868 number += numberInRow[iRow];
01869 if ( space < number ) {
01870 if ( !getRowSpace ( iRow, number ) ) {
01871 goto BAD_PIVOT;
01872 }
01873 }
01874
01875 putThis = putBase;
01876 next = nextRow[iRow];
01877 number = numberInRow[iRow];
01878 CoinBigIndex end = startRow[iRow] + number;
01879 int saveIndex;
01880
01881 saveIndex = indexColumn[startRow[next]];
01882
01883
01884 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01885 unsigned int test = *putThis;
01886
01887 putThis += increment2;
01888 test = 1 - ( ( test >> bit ) & 1 );
01889
01890 indexColumn[end] = saveColumn[jColumn];
01891 end += test;
01892 }
01893 indexColumn[startRow[next]] = saveIndex;
01894 markRow[iRow] = FAC_UNSET;
01895 number = end - startRow[iRow];
01896 numberInRow[iRow] = number;
01897 deleteLink ( iRow );
01898 addLink ( iRow, number );
01899 }
01900 markRow[iPivotRow] = FAC_UNSET;
01901
01902 deleteLink ( iPivotRow );
01903 deleteLink ( iPivotColumn + numberRows_ );
01904 totalElements_ += added;
01905 goodPivot= true;
01906
01907 }
01908 BAD_PIVOT:
01909
01910 ;
01911 }
01912 #undef FAC_UNSET
01913 #endif