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
00269
00270 inline double slackValue ( ) const {
00271 return slackValue_ ;
00272 }
00273 void slackValue ( double value );
00274 #endif
00275
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 assert (put<lengthAreaU_);
00997 elementU[put] = value;
00998 indexRowU[put] = indexL[j];
00999 if ( absValue > largest ) {
01000 largest = absValue;
01001 positionLargest = put;
01002 }
01003 put++;
01004 } else {
01005 work[j] = 0.0;
01006 added--;
01007 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01008 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01009
01010 if ( temp2[word] & ( 1 << bit ) ) {
01011
01012 iRow = indexL[j];
01013 CoinBigIndex start = startRowU[iRow];
01014 CoinBigIndex end = start + numberInRow[iRow];
01015 CoinBigIndex where = start;
01016
01017 while ( indexColumnU[where] != iColumn ) {
01018 where++;
01019 }
01020 #if DEBUG_COIN
01021 if ( where >= end ) {
01022 abort ( );
01023 }
01024 #endif
01025 indexColumnU[where] = indexColumnU[end - 1];
01026 numberInRow[iRow]--;
01027 } else {
01028
01029 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01030 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01031
01032 temp2[word] = temp2[word] | ( 1 << bit );
01033 }
01034 }
01035 }
01036 numberInColumn[iColumn] = put - startColumn;
01037
01038 if ( positionLargest >= 0 ) {
01039 value = elementU[positionLargest];
01040 iRow = indexRowU[positionLargest];
01041 elementU[positionLargest] = elementU[startColumn];
01042 indexRowU[positionLargest] = indexRowU[startColumn];
01043 elementU[startColumn] = value;
01044 indexRowU[startColumn] = iRow;
01045 }
01046
01047 if ( nextCount[iColumn + numberRows_] != -2 ) {
01048
01049 deleteLink ( iColumn + numberRows_ );
01050 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01051 }
01052 temp2 += increment2;
01053 }
01054
01055 unsigned int *putBase = workArea2;
01056 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01057 int i = 0;
01058
01059
01060 while ( bigLoops ) {
01061 bigLoops--;
01062 int bit;
01063 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01064 unsigned int *putThis = putBase;
01065 int iRow = indexL[i];
01066
01067
01068 int number = 0;
01069 int jColumn;
01070
01071 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01072 unsigned int test = *putThis;
01073
01074 putThis += increment2;
01075 test = 1 - ( ( test >> bit ) & 1 );
01076 number += test;
01077 }
01078 int next = nextRow[iRow];
01079 CoinBigIndex space;
01080
01081 space = startRowU[next] - startRowU[iRow];
01082 number += numberInRow[iRow];
01083 if ( space < number ) {
01084 if ( !getRowSpace ( iRow, number ) ) {
01085 return false;
01086 }
01087 }
01088
01089 putThis = putBase;
01090 next = nextRow[iRow];
01091 number = numberInRow[iRow];
01092 CoinBigIndex end = startRowU[iRow] + number;
01093 int saveIndex = indexColumnU[startRowU[next]];
01094
01095
01096 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01097 unsigned int test = *putThis;
01098
01099 putThis += increment2;
01100 test = 1 - ( ( test >> bit ) & 1 );
01101 indexColumnU[end] = saveColumn[jColumn];
01102 end += test;
01103 }
01104
01105 indexColumnU[startRowU[next]] = saveIndex;
01106 markRow[iRow] = static_cast<T>(largeInteger+1);
01107 number = end - startRowU[iRow];
01108 numberInRow[iRow] = number;
01109 deleteLink ( iRow );
01110 addLink ( iRow, number );
01111 }
01112 putBase++;
01113 }
01114 int bit;
01115
01116 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01117 unsigned int *putThis = putBase;
01118 int iRow = indexL[i];
01119
01120
01121 int number = 0;
01122 int jColumn;
01123
01124 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01125 unsigned int test = *putThis;
01126
01127 putThis += increment2;
01128 test = 1 - ( ( test >> bit ) & 1 );
01129 number += test;
01130 }
01131 int next = nextRow[iRow];
01132 CoinBigIndex space;
01133
01134 space = startRowU[next] - startRowU[iRow];
01135 number += numberInRow[iRow];
01136 if ( space < number ) {
01137 if ( !getRowSpace ( iRow, number ) ) {
01138 return false;
01139 }
01140 }
01141
01142 putThis = putBase;
01143 next = nextRow[iRow];
01144 number = numberInRow[iRow];
01145 CoinBigIndex end = startRowU[iRow] + number;
01146 int saveIndex;
01147
01148 saveIndex = indexColumnU[startRowU[next]];
01149
01150
01151 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01152 unsigned int test = *putThis;
01153
01154 putThis += increment2;
01155 test = 1 - ( ( test >> bit ) & 1 );
01156
01157 indexColumnU[end] = saveColumn[jColumn];
01158 end += test;
01159 }
01160 indexColumnU[startRowU[next]] = saveIndex;
01161 markRow[iRow] = static_cast<T>(largeInteger+1);
01162 number = end - startRowU[iRow];
01163 numberInRow[iRow] = number;
01164 deleteLink ( iRow );
01165 addLink ( iRow, number );
01166 }
01167 markRow[pivotRow] = static_cast<T>(largeInteger+1);
01168
01169 deleteLink ( pivotRow );
01170 deleteLink ( pivotColumn + numberRows_ );
01171 totalElements_ += added;
01172 return true;
01173 }
01174
01175
01177
01178 protected:
01179
01182
01183 double pivotTolerance_;
01185 double zeroTolerance_;
01186 #ifndef COIN_FAST_CODE
01187
01188 double slackValue_;
01189 #else
01190 #ifndef slackValue_
01191 #define slackValue_ -1.0
01192 #endif
01193 #endif
01194
01195 double areaFactor_;
01197 double relaxCheck_;
01199 int numberRows_;
01201 int numberRowsExtra_;
01203 int maximumRowsExtra_;
01205 int numberColumns_;
01207 int numberColumnsExtra_;
01209 int maximumColumnsExtra_;
01211 int numberGoodU_;
01213 int numberGoodL_;
01215 int maximumPivots_;
01217 int numberPivots_;
01220 CoinBigIndex totalElements_;
01222 CoinBigIndex factorElements_;
01224 CoinIntArrayWithLength pivotColumn_;
01226 CoinIntArrayWithLength permute_;
01228 CoinIntArrayWithLength permuteBack_;
01230 CoinIntArrayWithLength pivotColumnBack_;
01232 int status_;
01233
01238
01239
01241 int numberTrials_;
01243 CoinBigIndexArrayWithLength startRowU_;
01244
01246 CoinIntArrayWithLength numberInRow_;
01247
01249 CoinIntArrayWithLength numberInColumn_;
01250
01252 CoinIntArrayWithLength numberInColumnPlus_;
01253
01256 CoinIntArrayWithLength firstCount_;
01257
01259 CoinIntArrayWithLength nextCount_;
01260
01262 CoinIntArrayWithLength lastCount_;
01263
01265 CoinIntArrayWithLength nextColumn_;
01266
01268 CoinIntArrayWithLength lastColumn_;
01269
01271 CoinIntArrayWithLength nextRow_;
01272
01274 CoinIntArrayWithLength lastRow_;
01275
01277 CoinIntArrayWithLength saveColumn_;
01278
01280 CoinIntArrayWithLength markRow_;
01281
01283 int messageLevel_;
01284
01286 int biggerDimension_;
01287
01289 CoinIntArrayWithLength indexColumnU_;
01290
01292 CoinIntArrayWithLength pivotRowL_;
01293
01295 CoinFactorizationDoubleArrayWithLength pivotRegion_;
01296
01298 int numberSlacks_;
01299
01301 int numberU_;
01302
01304 CoinBigIndex maximumU_;
01305
01307
01308
01310 CoinBigIndex lengthU_;
01311
01313 CoinBigIndex lengthAreaU_;
01314
01316 CoinFactorizationDoubleArrayWithLength elementU_;
01317
01319 CoinIntArrayWithLength indexRowU_;
01320
01322 CoinBigIndexArrayWithLength startColumnU_;
01323
01325 CoinBigIndexArrayWithLength convertRowToColumnU_;
01326
01328 CoinBigIndex numberL_;
01329
01331 CoinBigIndex baseL_;
01332
01334 CoinBigIndex lengthL_;
01335
01337 CoinBigIndex lengthAreaL_;
01338
01340 CoinFactorizationDoubleArrayWithLength elementL_;
01341
01343 CoinIntArrayWithLength indexRowL_;
01344
01346 CoinBigIndexArrayWithLength startColumnL_;
01347
01349 bool doForrestTomlin_;
01350
01352 int numberR_;
01353
01355 CoinBigIndex lengthR_;
01356
01358 CoinBigIndex lengthAreaR_;
01359
01361 CoinFactorizationDouble *elementR_;
01362
01364 int *indexRowR_;
01365
01367 CoinBigIndexArrayWithLength startColumnR_;
01368
01370 double * denseArea_;
01371
01373 int * densePermute_;
01374
01376 int numberDense_;
01377
01379 int denseThreshold_;
01380
01382 CoinFactorizationDoubleArrayWithLength workArea_;
01383
01385 CoinUnsignedIntArrayWithLength workArea2_;
01386
01388 CoinBigIndex numberCompressions_;
01389
01391 mutable double ftranCountInput_;
01392 mutable double ftranCountAfterL_;
01393 mutable double ftranCountAfterR_;
01394 mutable double ftranCountAfterU_;
01395 mutable double btranCountInput_;
01396 mutable double btranCountAfterU_;
01397 mutable double btranCountAfterR_;
01398 mutable double btranCountAfterL_;
01399
01401 mutable int numberFtranCounts_;
01402 mutable int numberBtranCounts_;
01403
01405 double ftranAverageAfterL_;
01406 double ftranAverageAfterR_;
01407 double ftranAverageAfterU_;
01408 double btranAverageAfterU_;
01409 double btranAverageAfterR_;
01410 double btranAverageAfterL_;
01411
01413 mutable bool collectStatistics_;
01414
01416 int sparseThreshold_;
01417
01419 int sparseThreshold2_;
01420
01422 CoinBigIndexArrayWithLength startRowL_;
01423
01425 CoinIntArrayWithLength indexColumnL_;
01426
01428 CoinFactorizationDoubleArrayWithLength elementByRowL_;
01429
01431 mutable CoinIntArrayWithLength sparse_;
01435 int biasLU_;
01441 int persistenceFlag_;
01443 };
01444
01445 #ifdef COIN_HAS_LAPACK
01446 #define DENSE_CODE 1
01447
01448 #ifndef ipfint
01449
01450 typedef int ipfint;
01451 typedef const int cipfint;
01452 #endif
01453 #endif
01454 #endif
01455
01456 #ifdef UGLY_COIN_FACTOR_CODING
01457 #define FAC_UNSET (FAC_SET+1)
01458 {
01459 goodPivot=false;
01460
01461 CoinBigIndex startColumnThis = startColumn[iPivotColumn];
01462 CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
01463 int put = 0;
01464 CoinBigIndex startRowThis = startRow[iPivotRow];
01465 CoinBigIndex endRow = startRowThis + numberDoRow + 1;
01466 if ( pivotColumnPosition < 0 ) {
01467 for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01468 int iColumn = indexColumn[pivotColumnPosition];
01469 if ( iColumn != iPivotColumn ) {
01470 saveColumn[put++] = iColumn;
01471 } else {
01472 break;
01473 }
01474 }
01475 } else {
01476 for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
01477 saveColumn[put++] = indexColumn[i];
01478 }
01479 }
01480 assert (pivotColumnPosition<endRow);
01481 assert (indexColumn[pivotColumnPosition]==iPivotColumn);
01482 pivotColumnPosition++;
01483 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01484 saveColumn[put++] = indexColumn[pivotColumnPosition];
01485 }
01486
01487 int next = nextRow[iPivotRow];
01488 int last = lastRow[iPivotRow];
01489
01490 nextRow[last] = next;
01491 lastRow[next] = last;
01492 nextRow[iPivotRow] = numberGoodU_;
01493 lastRow[iPivotRow] = -2;
01494 numberInRow[iPivotRow] = 0;
01495
01496 CoinBigIndex l = lengthL_;
01497
01498 {
01499 if ( l + numberDoColumn > lengthAreaL_ ) {
01500
01501 if ((messageLevel_&4)!=0)
01502 printf("more memory needed in middle of invert\n");
01503 goto BAD_PIVOT;
01504 }
01505
01506 CoinBigIndex lSave = l;
01507
01508 CoinBigIndex * startColumnL = startColumnL_.array();
01509 startColumnL[numberGoodL_] = l;
01510 numberGoodL_++;
01511 startColumnL[numberGoodL_] = l + numberDoColumn;
01512 lengthL_ += numberDoColumn;
01513 if ( pivotRowPosition < 0 ) {
01514 for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01515 int iRow = indexRow[pivotRowPosition];
01516 if ( iRow != iPivotRow ) {
01517 indexRowL[l] = iRow;
01518 elementL[l] = element[pivotRowPosition];
01519 markRow[iRow] = l - lSave;
01520 l++;
01521
01522 CoinBigIndex start = startRow[iRow];
01523 CoinBigIndex end = start + numberInRow[iRow];
01524 CoinBigIndex where = start;
01525
01526 while ( indexColumn[where] != iPivotColumn ) {
01527 where++;
01528 }
01529 #if DEBUG_COIN
01530 if ( where >= end ) {
01531 abort ( );
01532 }
01533 #endif
01534 indexColumn[where] = indexColumn[end - 1];
01535 numberInRow[iRow]--;
01536 } else {
01537 break;
01538 }
01539 }
01540 } else {
01541 CoinBigIndex i;
01542
01543 for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
01544 int iRow = indexRow[i];
01545
01546 markRow[iRow] = l - lSave;
01547 indexRowL[l] = iRow;
01548 elementL[l] = element[i];
01549 l++;
01550
01551 CoinBigIndex start = startRow[iRow];
01552 CoinBigIndex end = start + numberInRow[iRow];
01553 CoinBigIndex where = start;
01554
01555 while ( indexColumn[where] != iPivotColumn ) {
01556 where++;
01557 }
01558 #if DEBUG_COIN
01559 if ( where >= end ) {
01560 abort ( );
01561 }
01562 #endif
01563 indexColumn[where] = indexColumn[end - 1];
01564 numberInRow[iRow]--;
01565 assert (numberInRow[iRow]>=0);
01566 }
01567 }
01568 assert (pivotRowPosition<endColumn);
01569 assert (indexRow[pivotRowPosition]==iPivotRow);
01570 CoinFactorizationDouble pivotElement = element[pivotRowPosition];
01571 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
01572
01573 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
01574 pivotRowPosition++;
01575 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01576 int iRow = indexRow[pivotRowPosition];
01577
01578 markRow[iRow] = l - lSave;
01579 indexRowL[l] = iRow;
01580 elementL[l] = element[pivotRowPosition];
01581 l++;
01582
01583 CoinBigIndex start = startRow[iRow];
01584 CoinBigIndex end = start + numberInRow[iRow];
01585 CoinBigIndex where = start;
01586
01587 while ( indexColumn[where] != iPivotColumn ) {
01588 where++;
01589 }
01590 #if DEBUG_COIN
01591 if ( where >= end ) {
01592 abort ( );
01593 }
01594 #endif
01595 indexColumn[where] = indexColumn[end - 1];
01596 numberInRow[iRow]--;
01597 assert (numberInRow[iRow]>=0);
01598 }
01599 markRow[iPivotRow] = FAC_SET;
01600
01601 numberInColumn[iPivotColumn] = 0;
01602
01603 int *indexL = &indexRowL[lSave];
01604 CoinFactorizationDouble *multipliersL = &elementL[lSave];
01605
01606
01607 int j;
01608
01609 for ( j = 0; j < numberDoColumn; j++ ) {
01610 multipliersL[j] *= pivotMultiplier;
01611 }
01612
01613 CoinBigIndex iErase;
01614 for ( iErase = 0; iErase < increment2 * numberDoRow;
01615 iErase++ ) {
01616 workArea2[iErase] = 0;
01617 }
01618 CoinBigIndex added = numberDoRow * numberDoColumn;
01619 unsigned int *temp2 = workArea2;
01620 int * nextColumn = nextColumn_.array();
01621
01622
01623 int jColumn;
01624 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01625 int iColumn = saveColumn[jColumn];
01626 CoinBigIndex startColumnThis = startColumn[iColumn];
01627 CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
01628 int iRow = indexRow[startColumnThis];
01629 CoinFactorizationDouble value = element[startColumnThis];
01630 double largest;
01631 CoinBigIndex put = startColumnThis;
01632 CoinBigIndex positionLargest = -1;
01633 CoinFactorizationDouble thisPivotValue = 0.0;
01634
01635
01636 bool checkLargest;
01637 int mark = markRow[iRow];
01638
01639 if ( mark == FAC_UNSET ) {
01640 largest = fabs ( value );
01641 positionLargest = put;
01642 put++;
01643 checkLargest = false;
01644 } else {
01645
01646 largest = 0.0;
01647 checkLargest = true;
01648 if ( mark != FAC_SET ) {
01649
01650 workArea[mark] = value;
01651 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01652 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01653
01654 temp2[word] = temp2[word] | ( 1 << bit );
01655 added--;
01656 } else {
01657 thisPivotValue = value;
01658 }
01659 }
01660 CoinBigIndex i;
01661 for ( i = startColumnThis + 1; i < endColumn; i++ ) {
01662 iRow = indexRow[i];
01663 value = element[i];
01664 int mark = markRow[iRow];
01665
01666 if ( mark == FAC_UNSET ) {
01667
01668 indexRow[put] = iRow;
01669 element[put] = value;
01670 if ( checkLargest ) {
01671 double absValue = fabs ( value );
01672
01673 if ( absValue > largest ) {
01674 largest = absValue;
01675 positionLargest = put;
01676 }
01677 }
01678 put++;
01679 } else if ( mark != FAC_SET ) {
01680
01681 workArea[mark] = value;
01682 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01683 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01684
01685 temp2[word] = temp2[word] | ( 1 << bit );
01686 added--;
01687 } else {
01688 thisPivotValue = value;
01689 }
01690 }
01691
01692 element[put] = element[startColumnThis];
01693 indexRow[put] = indexRow[startColumnThis];
01694 if ( positionLargest == startColumnThis ) {
01695 positionLargest = put;
01696 }
01697 put++;
01698 element[startColumnThis] = thisPivotValue;
01699 indexRow[startColumnThis] = iPivotRow;
01700
01701 startColumnThis++;
01702 numberInColumn[iColumn] = put - startColumnThis;
01703 int * numberInColumnPlus = numberInColumnPlus_.array();
01704 numberInColumnPlus[iColumn]++;
01705 startColumn[iColumn]++;
01706
01707 int next = nextColumn[iColumn];
01708 CoinBigIndex space;
01709
01710 space = startColumn[next] - put - numberInColumnPlus[next];
01711
01712 if ( numberDoColumn > space ) {
01713
01714 if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
01715 goto BAD_PIVOT;
01716 }
01717
01718 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
01719 startColumnThis = startColumn[iColumn];
01720 put = startColumnThis + numberInColumn[iColumn];
01721 }
01722 double tolerance = zeroTolerance_;
01723
01724 int *nextCount = nextCount_.array();
01725 for ( j = 0; j < numberDoColumn; j++ ) {
01726 value = workArea[j] - thisPivotValue * multipliersL[j];
01727 double absValue = fabs ( value );
01728
01729 if ( absValue > tolerance ) {
01730 workArea[j] = 0.0;
01731 element[put] = value;
01732 indexRow[put] = indexL[j];
01733 if ( absValue > largest ) {
01734 largest = absValue;
01735 positionLargest = put;
01736 }
01737 put++;
01738 } else {
01739 workArea[j] = 0.0;
01740 added--;
01741 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01742 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01743
01744 if ( temp2[word] & ( 1 << bit ) ) {
01745
01746 iRow = indexL[j];
01747 CoinBigIndex start = startRow[iRow];
01748 CoinBigIndex end = start + numberInRow[iRow];
01749 CoinBigIndex where = start;
01750
01751 while ( indexColumn[where] != iColumn ) {
01752 where++;
01753 }
01754 #if DEBUG_COIN
01755 if ( where >= end ) {
01756 abort ( );
01757 }
01758 #endif
01759 indexColumn[where] = indexColumn[end - 1];
01760 numberInRow[iRow]--;
01761 } else {
01762
01763 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01764 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01765
01766 temp2[word] = temp2[word] | ( 1 << bit );
01767 }
01768 }
01769 }
01770 numberInColumn[iColumn] = put - startColumnThis;
01771
01772 if ( positionLargest >= 0 ) {
01773 value = element[positionLargest];
01774 iRow = indexRow[positionLargest];
01775 element[positionLargest] = element[startColumnThis];
01776 indexRow[positionLargest] = indexRow[startColumnThis];
01777 element[startColumnThis] = value;
01778 indexRow[startColumnThis] = iRow;
01779 }
01780
01781 if ( nextCount[iColumn + numberRows_] != -2 ) {
01782
01783 deleteLink ( iColumn + numberRows_ );
01784 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01785 }
01786 temp2 += increment2;
01787 }
01788
01789 unsigned int *putBase = workArea2;
01790 int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01791 int i = 0;
01792
01793
01794 while ( bigLoops ) {
01795 bigLoops--;
01796 int bit;
01797 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01798 unsigned int *putThis = putBase;
01799 int iRow = indexL[i];
01800
01801
01802 int number = 0;
01803 int jColumn;
01804
01805 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01806 unsigned int test = *putThis;
01807
01808 putThis += increment2;
01809 test = 1 - ( ( test >> bit ) & 1 );
01810 number += test;
01811 }
01812 int next = nextRow[iRow];
01813 CoinBigIndex space;
01814
01815 space = startRow[next] - startRow[iRow];
01816 number += numberInRow[iRow];
01817 if ( space < number ) {
01818 if ( !getRowSpace ( iRow, number ) ) {
01819 goto BAD_PIVOT;
01820 }
01821 }
01822
01823 putThis = putBase;
01824 next = nextRow[iRow];
01825 number = numberInRow[iRow];
01826 CoinBigIndex end = startRow[iRow] + number;
01827 int saveIndex = indexColumn[startRow[next]];
01828
01829
01830 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01831 unsigned int test = *putThis;
01832
01833 putThis += increment2;
01834 test = 1 - ( ( test >> bit ) & 1 );
01835 indexColumn[end] = saveColumn[jColumn];
01836 end += test;
01837 }
01838
01839 indexColumn[startRow[next]] = saveIndex;
01840 markRow[iRow] = FAC_UNSET;
01841 number = end - startRow[iRow];
01842 numberInRow[iRow] = number;
01843 deleteLink ( iRow );
01844 addLink ( iRow, number );
01845 }
01846 putBase++;
01847 }
01848 int bit;
01849
01850 for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
01851 unsigned int *putThis = putBase;
01852 int iRow = indexL[i];
01853
01854
01855 int number = 0;
01856 int jColumn;
01857
01858 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01859 unsigned int test = *putThis;
01860
01861 putThis += increment2;
01862 test = 1 - ( ( test >> bit ) & 1 );
01863 number += test;
01864 }
01865 int next = nextRow[iRow];
01866 CoinBigIndex space;
01867
01868 space = startRow[next] - startRow[iRow];
01869 number += numberInRow[iRow];
01870 if ( space < number ) {
01871 if ( !getRowSpace ( iRow, number ) ) {
01872 goto BAD_PIVOT;
01873 }
01874 }
01875
01876 putThis = putBase;
01877 next = nextRow[iRow];
01878 number = numberInRow[iRow];
01879 CoinBigIndex end = startRow[iRow] + number;
01880 int saveIndex;
01881
01882 saveIndex = indexColumn[startRow[next]];
01883
01884
01885 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01886 unsigned int test = *putThis;
01887
01888 putThis += increment2;
01889 test = 1 - ( ( test >> bit ) & 1 );
01890
01891 indexColumn[end] = saveColumn[jColumn];
01892 end += test;
01893 }
01894 indexColumn[startRow[next]] = saveIndex;
01895 markRow[iRow] = FAC_UNSET;
01896 number = end - startRow[iRow];
01897 numberInRow[iRow] = number;
01898 deleteLink ( iRow );
01899 addLink ( iRow, number );
01900 }
01901 markRow[iPivotRow] = FAC_UNSET;
01902
01903 deleteLink ( iPivotRow );
01904 deleteLink ( iPivotColumn + numberRows_ );
01905 totalElements_ += added;
01906 goodPivot= true;
01907
01908 }
01909 BAD_PIVOT:
01910
01911 ;
01912 }
01913 #undef FAC_UNSET
01914 #endif