00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef CoinFactorization_H
00013 #define CoinFactorization_H
00014
00015
00016 #include <iostream>
00017 #include <string>
00018 #include <cassert>
00019 #include <cstdio>
00020 #include <cmath>
00021 #include "CoinTypes.hpp"
00022 #include "CoinIndexedVector.hpp"
00023
00024 class CoinPackedMatrix;
00050 class CoinFactorization {
00051 friend void CoinFactorizationUnitTest( const std::string & mpsDir );
00052
00053 public:
00054
00057
00058 CoinFactorization ( );
00060 CoinFactorization ( const CoinFactorization &other);
00061
00063 ~CoinFactorization ( );
00065 void almostDestructor();
00067 void show_self ( ) const;
00069 int saveFactorization (const char * file ) const;
00073 int restoreFactorization (const char * file , bool factor=false) ;
00075 void sort ( ) const;
00077 CoinFactorization & operator = ( const CoinFactorization & other );
00079
00089 int factorize ( const CoinPackedMatrix & matrix,
00090 int rowIsBasic[], int columnIsBasic[] ,
00091 double areaFactor = 0.0 );
00102 int factorize ( int numberRows,
00103 int numberColumns,
00104 CoinBigIndex numberElements,
00105 CoinBigIndex maximumL,
00106 CoinBigIndex maximumU,
00107 const int indicesRow[],
00108 const int indicesColumn[], const double elements[] ,
00109 int permutation[],
00110 double areaFactor = 0.0);
00115 int factorizePart1 ( int numberRows,
00116 int numberColumns,
00117 CoinBigIndex estimateNumberElements,
00118 int * indicesRow[],
00119 int * indicesColumn[],
00120 CoinFactorizationDouble * elements[],
00121 double areaFactor = 0.0);
00128 int factorizePart2 (int permutation[],int exactNumberElements);
00130 double conditionNumber() const;
00131
00133
00136
00137 inline int status ( ) const {
00138 return status_;
00139 }
00141 inline void setStatus ( int value)
00142 { status_=value; }
00144 inline int pivots ( ) const {
00145 return numberPivots_;
00146 }
00148 inline void setPivots ( int value )
00149 { numberPivots_=value; }
00151 inline int *permute ( ) const {
00152 return permute_.array();
00153 }
00155 inline int *pivotColumn ( ) const {
00156 return pivotColumn_.array();
00157 }
00159 inline CoinFactorizationDouble *pivotRegion ( ) const {
00160 return pivotRegion_.array();
00161 }
00163 inline int *permuteBack ( ) const {
00164 return permuteBack_.array();
00165 }
00168 inline int *pivotColumnBack ( ) const {
00169
00170 return pivotColumnBack_.array();
00171 }
00173 inline CoinBigIndex * startRowL() const
00174 { return startRowL_.array();}
00175
00177 inline CoinBigIndex * startColumnL() const
00178 { return startColumnL_.array();}
00179
00181 inline int * indexColumnL() const
00182 { return indexColumnL_.array();}
00183
00185 inline int * indexRowL() const
00186 { return indexRowL_.array();}
00187
00189 inline CoinFactorizationDouble * elementByRowL() const
00190 { return elementByRowL_.array();}
00191
00193 inline int numberRowsExtra ( ) const {
00194 return numberRowsExtra_;
00195 }
00197 inline void setNumberRows(int value)
00198 { numberRows_ = value; }
00200 inline int numberRows ( ) const {
00201 return numberRows_;
00202 }
00204 inline CoinBigIndex numberL() const
00205 { return numberL_;}
00206
00208 inline CoinBigIndex baseL() const
00209 { return baseL_;}
00211 inline int maximumRowsExtra ( ) const {
00212 return maximumRowsExtra_;
00213 }
00215 inline int numberColumns ( ) const {
00216 return numberColumns_;
00217 }
00219 inline int numberElements ( ) const {
00220 return totalElements_;
00221 }
00223 inline int numberForrestTomlin ( ) const {
00224 return numberInColumn_.array()[numberColumnsExtra_];
00225 }
00227 inline int numberGoodColumns ( ) const {
00228 return numberGoodU_;
00229 }
00231 inline double areaFactor ( ) const {
00232 return areaFactor_;
00233 }
00234 inline void areaFactor ( double value ) {
00235 areaFactor_=value;
00236 }
00238 double adjustedAreaFactor() const;
00240 inline void relaxAccuracyCheck(double value)
00241 { relaxCheck_ = value;}
00242 inline double getAccuracyCheck() const
00243 { return relaxCheck_;}
00245 inline int messageLevel ( ) const {
00246 return messageLevel_ ;
00247 }
00248 void messageLevel ( int value );
00250 inline int maximumPivots ( ) const {
00251 return maximumPivots_ ;
00252 }
00253 void maximumPivots ( int value );
00254
00256 inline int denseThreshold() const
00257 { return denseThreshold_;}
00259 inline void setDenseThreshold(int value)
00260 { denseThreshold_ = value;}
00262 inline double pivotTolerance ( ) const {
00263 return pivotTolerance_ ;
00264 }
00265 void pivotTolerance ( double value );
00267 inline double zeroTolerance ( ) const {
00268 return zeroTolerance_ ;
00269 }
00270 void zeroTolerance ( double value );
00271 #ifndef COIN_FAST_CODE
00273 inline double slackValue ( ) const {
00274 return slackValue_ ;
00275 }
00276 void slackValue ( double value );
00277 #endif
00279 double maximumCoefficient() const;
00281 inline bool forrestTomlin() const
00282 { return doForrestTomlin_;}
00283 inline void setForrestTomlin(bool value)
00284 { doForrestTomlin_=value;}
00286 inline bool spaceForForrestTomlin() const
00287 {
00288 CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
00289 CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
00290 return (space>=0)&&doForrestTomlin_;
00291 }
00293
00296
00298 inline int numberDense() const
00299 { return numberDense_;}
00300
00302 inline CoinBigIndex numberElementsU ( ) const {
00303 return lengthU_;
00304 }
00306 inline void setNumberElementsU(CoinBigIndex value)
00307 { lengthU_ = value; }
00309 inline CoinBigIndex lengthAreaU ( ) const {
00310 return lengthAreaU_;
00311 }
00313 inline CoinBigIndex numberElementsL ( ) const {
00314 return lengthL_;
00315 }
00317 inline CoinBigIndex lengthAreaL ( ) const {
00318 return lengthAreaL_;
00319 }
00321 inline CoinBigIndex numberElementsR ( ) const {
00322 return lengthR_;
00323 }
00325 inline CoinBigIndex numberCompressions() const
00326 { return numberCompressions_;}
00328 inline int * numberInRow() const
00329 { return numberInRow_.array();}
00331 inline int * numberInColumn() const
00332 { return numberInColumn_.array();}
00334 inline CoinFactorizationDouble * elementU() const
00335 { return elementU_.array();}
00337 inline int * indexRowU() const
00338 { return indexRowU_.array();}
00340 inline CoinBigIndex * startColumnU() const
00341 { return startColumnU_.array();}
00343 inline int maximumColumnsExtra()
00344 { return maximumColumnsExtra_;}
00348 inline int biasLU() const
00349 { return biasLU_;}
00350 inline void setBiasLU(int value)
00351 { biasLU_=value;}
00357 inline int persistenceFlag() const
00358 { return persistenceFlag_;}
00359 void setPersistenceFlag(int value);
00361
00364
00372 int replaceColumn ( CoinIndexedVector * regionSparse,
00373 int pivotRow,
00374 double pivotCheck ,
00375 bool checkBeforeModifying=false,
00376 double acceptablePivot=1.0e-8);
00381 void replaceColumnU ( CoinIndexedVector * regionSparse,
00382 CoinBigIndex * deleted,
00383 int internalPivotRow);
00385
00395 int updateColumnFT ( CoinIndexedVector * regionSparse,
00396 CoinIndexedVector * regionSparse2);
00399 int updateColumn ( CoinIndexedVector * regionSparse,
00400 CoinIndexedVector * regionSparse2,
00401 bool noPermute=false) const;
00407 int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
00408 CoinIndexedVector * regionSparse2,
00409 CoinIndexedVector * regionSparse3,
00410 bool noPermuteRegion3=false) ;
00415 int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00416 CoinIndexedVector * regionSparse2) const;
00418 void goSparse();
00420 inline int sparseThreshold ( ) const
00421 { return sparseThreshold_;}
00423 void sparseThreshold ( int value );
00425
00426
00430
00431 inline void clearArrays()
00432 { gutsOfDestructor();}
00434
00437
00440 int add ( CoinBigIndex numberElements,
00441 int indicesRow[],
00442 int indicesColumn[], double elements[] );
00443
00446 int addColumn ( CoinBigIndex numberElements,
00447 int indicesRow[], double elements[] );
00448
00451 int addRow ( CoinBigIndex numberElements,
00452 int indicesColumn[], double elements[] );
00453
00455 int deleteColumn ( int Row );
00457 int deleteRow ( int Row );
00458
00462 int replaceRow ( int whichRow, int numberElements,
00463 const int indicesColumn[], const double elements[] );
00465 void emptyRows(int numberToEmpty, const int which[]);
00467
00468
00469 void checkSparse();
00471 inline bool collectStatistics() const
00472 { return collectStatistics_;}
00474 inline void setCollectStatistics(bool onOff) const
00475 { collectStatistics_ = onOff;}
00477 void gutsOfDestructor(int type=1);
00479 void gutsOfInitialize(int type);
00480 void gutsOfCopy(const CoinFactorization &other);
00481
00483 void resetStatistics();
00484
00485
00487
00489
00490 void getAreas ( int numberRows,
00491 int numberColumns,
00492 CoinBigIndex maximumL,
00493 CoinBigIndex maximumU );
00494
00497 void preProcess ( int state,
00498 int possibleDuplicates = -1 );
00500 int factor ( );
00501 protected:
00504 int factorSparse ( );
00507 int factorSparseSmall ( );
00510 int factorSparseLarge ( );
00513 int factorDense ( );
00514
00516 bool pivotOneOtherRow ( int pivotRow,
00517 int pivotColumn );
00519 bool pivotRowSingleton ( int pivotRow,
00520 int pivotColumn );
00522 bool pivotColumnSingleton ( int pivotRow,
00523 int pivotColumn );
00524
00529 bool getColumnSpace ( int iColumn,
00530 int extraNeeded );
00531
00534 bool reorderU();
00538 bool getColumnSpaceIterateR ( int iColumn, double value,
00539 int iRow);
00545 CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00546 int iRow);
00550 bool getRowSpace ( int iRow, int extraNeeded );
00551
00555 bool getRowSpaceIterate ( int iRow,
00556 int extraNeeded );
00558 void checkConsistency ( );
00560 inline void addLink ( int index, int count ) {
00561 int *nextCount = nextCount_.array();
00562 int *firstCount = firstCount_.array();
00563 int *lastCount = lastCount_.array();
00564 int next = firstCount[count];
00565 lastCount[index] = -2 - count;
00566 if ( next < 0 ) {
00567
00568 firstCount[count] = index;
00569 nextCount[index] = -1;
00570 } else {
00571 firstCount[count] = index;
00572 nextCount[index] = next;
00573 lastCount[next] = index;
00574 }}
00576 inline void deleteLink ( int index ) {
00577 int *nextCount = nextCount_.array();
00578 int *firstCount = firstCount_.array();
00579 int *lastCount = lastCount_.array();
00580 int next = nextCount[index];
00581 int last = lastCount[index];
00582 if ( last >= 0 ) {
00583 nextCount[last] = next;
00584 } else {
00585 int count = -last - 2;
00586
00587 firstCount[count] = next;
00588 }
00589 if ( next >= 0 ) {
00590 lastCount[next] = last;
00591 }
00592 nextCount[index] = -2;
00593 lastCount[index] = -2;
00594 return;
00595 }
00597 void separateLinks(int count,bool rowsFirst);
00599 void cleanup ( );
00600
00602 void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00604 void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00606 void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00608 void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00609
00611 void updateColumnR ( CoinIndexedVector * region ) const;
00614 void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00615
00617 void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00618
00620 void updateColumnUSparse ( CoinIndexedVector * regionSparse,
00621 int * indexIn) const;
00623 void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
00624 int * indexIn) const;
00626 int updateColumnUDensish ( double * COIN_RESTRICT region,
00627 int * COIN_RESTRICT regionIndex) const;
00629 void updateTwoColumnsUDensish (
00630 int & numberNonZero1,
00631 double * COIN_RESTRICT region1,
00632 int * COIN_RESTRICT index1,
00633 int & numberNonZero2,
00634 double * COIN_RESTRICT region2,
00635 int * COIN_RESTRICT index2) const;
00637 void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
00639 void permuteBack ( CoinIndexedVector * regionSparse,
00640 CoinIndexedVector * outVector) const;
00641
00643 void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00646 void updateColumnTransposeU ( CoinIndexedVector * region,
00647 int smallestIndex) const;
00650 void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00651 int smallestIndex) const;
00654 void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00655 int smallestIndex) const;
00658 void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00661 void updateColumnTransposeUByColumn ( CoinIndexedVector * region,
00662 int smallestIndex) const;
00663
00665 void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00667 void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00669 void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00670
00672 void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00674 void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00676 void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00678 void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00680 void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00681 public:
00686 int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00687 int pivotRow, double alpha);
00688 protected:
00691 int checkPivot(double saveFromU, double oldPivot) const;
00692
00693 #ifdef INT_IS_8
00694 #define COINFACTORIZATION_BITS_PER_INT 64
00695 #define COINFACTORIZATION_SHIFT_PER_INT 6
00696 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00697 #else
00698 #define COINFACTORIZATION_BITS_PER_INT 32
00699 #define COINFACTORIZATION_SHIFT_PER_INT 5
00700 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00701 #endif
00702 template <class T> inline bool
00703 pivot ( int pivotRow,
00704 int pivotColumn,
00705 CoinBigIndex pivotRowPosition,
00706 CoinBigIndex pivotColumnPosition,
00707 CoinFactorizationDouble work[],
00708 unsigned int workArea2[],
00709 int increment2,
00710 T markRow[] ,
00711 int largeInteger)
00712 {
00713 int *indexColumnU = indexColumnU_.array();
00714 CoinBigIndex *startColumnU = startColumnU_.array();
00715 int *numberInColumn = numberInColumn_.array();
00716 CoinFactorizationDouble *elementU = elementU_.array();
00717 int *indexRowU = indexRowU_.array();
00718 CoinBigIndex *startRowU = startRowU_.array();
00719 int *numberInRow = numberInRow_.array();
00720 CoinFactorizationDouble *elementL = elementL_.array();
00721 int *indexRowL = indexRowL_.array();
00722 int *saveColumn = saveColumn_.array();
00723 int *nextRow = nextRow_.array();
00724 int *lastRow = lastRow_.array() ;
00725
00726
00727 int numberInPivotRow = numberInRow[pivotRow] - 1;
00728 CoinBigIndex startColumn = startColumnU[pivotColumn];
00729 int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00730 CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00731 int put = 0;
00732 CoinBigIndex startRow = startRowU[pivotRow];
00733 CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00734
00735 if ( pivotColumnPosition < 0 ) {
00736 for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00737 int iColumn = indexColumnU[pivotColumnPosition];
00738 if ( iColumn != pivotColumn ) {
00739 saveColumn[put++] = iColumn;
00740 } else {
00741 break;
00742 }
00743 }
00744 } else {
00745 for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00746 saveColumn[put++] = indexColumnU[i];
00747 }
00748 }
00749 assert (pivotColumnPosition<endRow);
00750 assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00751 pivotColumnPosition++;
00752 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00753 saveColumn[put++] = indexColumnU[pivotColumnPosition];
00754 }
00755
00756 int next = nextRow[pivotRow];
00757 int last = lastRow[pivotRow];
00758
00759 nextRow[last] = next;
00760 lastRow[next] = last;
00761 nextRow[pivotRow] = numberGoodU_;
00762 lastRow[pivotRow] = -2;
00763 numberInRow[pivotRow] = 0;
00764
00765 CoinBigIndex l = lengthL_;
00766
00767 if ( l + numberInPivotColumn > lengthAreaL_ ) {
00768
00769 if ((messageLevel_&4)!=0)
00770 printf("more memory needed in middle of invert\n");
00771 return false;
00772 }
00773
00774 CoinBigIndex lSave = l;
00775
00776 CoinBigIndex * startColumnL = startColumnL_.array();
00777 startColumnL[numberGoodL_] = l;
00778 numberGoodL_++;
00779 startColumnL[numberGoodL_] = l + numberInPivotColumn;
00780 lengthL_ += numberInPivotColumn;
00781 if ( pivotRowPosition < 0 ) {
00782 for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00783 int iRow = indexRowU[pivotRowPosition];
00784 if ( iRow != pivotRow ) {
00785 indexRowL[l] = iRow;
00786 elementL[l] = elementU[pivotRowPosition];
00787 markRow[iRow] = static_cast<T>(l - lSave);
00788 l++;
00789
00790 CoinBigIndex start = startRowU[iRow];
00791 CoinBigIndex end = start + numberInRow[iRow];
00792 CoinBigIndex where = start;
00793
00794 while ( indexColumnU[where] != pivotColumn ) {
00795 where++;
00796 }
00797 #if DEBUG_COIN
00798 if ( where >= end ) {
00799 abort ( );
00800 }
00801 #endif
00802 indexColumnU[where] = indexColumnU[end - 1];
00803 numberInRow[iRow]--;
00804 } else {
00805 break;
00806 }
00807 }
00808 } else {
00809 CoinBigIndex i;
00810
00811 for ( i = startColumn; i < pivotRowPosition; i++ ) {
00812 int iRow = indexRowU[i];
00813
00814 markRow[iRow] = static_cast<T>(l - lSave);
00815 indexRowL[l] = iRow;
00816 elementL[l] = elementU[i];
00817 l++;
00818
00819 CoinBigIndex start = startRowU[iRow];
00820 CoinBigIndex end = start + numberInRow[iRow];
00821 CoinBigIndex where = start;
00822
00823 while ( indexColumnU[where] != pivotColumn ) {
00824 where++;
00825 }
00826 #if DEBUG_COIN
00827 if ( where >= end ) {
00828 abort ( );
00829 }
00830 #endif
00831 indexColumnU[where] = indexColumnU[end - 1];
00832 numberInRow[iRow]--;
00833 assert (numberInRow[iRow]>=0);
00834 }
00835 }
00836 assert (pivotRowPosition<endColumn);
00837 assert (indexRowU[pivotRowPosition]==pivotRow);
00838 CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
00839 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
00840
00841 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00842 pivotRowPosition++;
00843 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00844 int iRow = indexRowU[pivotRowPosition];
00845
00846 markRow[iRow] = static_cast<T>(l - lSave);
00847 indexRowL[l] = iRow;
00848 elementL[l] = elementU[pivotRowPosition];
00849 l++;
00850
00851 CoinBigIndex start = startRowU[iRow];
00852 CoinBigIndex end = start + numberInRow[iRow];
00853 CoinBigIndex where = start;
00854
00855 while ( indexColumnU[where] != pivotColumn ) {
00856 where++;
00857 }
00858 #if DEBUG_COIN
00859 if ( where >= end ) {
00860 abort ( );
00861 }
00862 #endif
00863 indexColumnU[where] = indexColumnU[end - 1];
00864 numberInRow[iRow]--;
00865 assert (numberInRow[iRow]>=0);
00866 }
00867 markRow[pivotRow] = static_cast<T>(largeInteger);
00868
00869 numberInColumn[pivotColumn] = 0;
00870
00871 int *indexL = &indexRowL[lSave];
00872 CoinFactorizationDouble *multipliersL = &elementL[lSave];
00873
00874
00875 int j;
00876
00877 for ( j = 0; j < numberInPivotColumn; j++ ) {
00878 multipliersL[j] *= pivotMultiplier;
00879 }
00880
00881 CoinBigIndex iErase;
00882 for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00883 iErase++ ) {
00884 workArea2[iErase] = 0;
00885 }
00886 CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00887 unsigned int *temp2 = workArea2;
00888 int * nextColumn = nextColumn_.array();
00889
00890
00891 int jColumn;
00892 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
00893 int iColumn = saveColumn[jColumn];
00894 CoinBigIndex startColumn = startColumnU[iColumn];
00895 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
00896 int iRow = indexRowU[startColumn];
00897 CoinFactorizationDouble value = elementU[startColumn];
00898 double largest;
00899 CoinBigIndex put = startColumn;
00900 CoinBigIndex positionLargest = -1;
00901 CoinFactorizationDouble thisPivotValue = 0.0;
00902
00903
00904 bool checkLargest;
00905 int mark = markRow[iRow];
00906
00907 if ( mark == largeInteger+1 ) {
00908 largest = fabs ( value );
00909 positionLargest = put;
00910 put++;
00911 checkLargest = false;
00912 } else {
00913
00914 largest = 0.0;
00915 checkLargest = true;
00916 if ( mark != largeInteger ) {
00917
00918 work[mark] = value;
00919 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00920 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00921
00922 temp2[word] = temp2[word] | ( 1 << bit );
00923 added--;
00924 } else {
00925 thisPivotValue = value;
00926 }
00927 }
00928 CoinBigIndex i;
00929 for ( i = startColumn + 1; i < endColumn; i++ ) {
00930 iRow = indexRowU[i];
00931 value = elementU[i];
00932 int mark = markRow[iRow];
00933
00934 if ( mark == largeInteger+1 ) {
00935
00936 indexRowU[put] = iRow;
00937 elementU[put] = value;
00938 if ( checkLargest ) {
00939 double absValue = fabs ( value );
00940
00941 if ( absValue > largest ) {
00942 largest = absValue;
00943 positionLargest = put;
00944 }
00945 }
00946 put++;
00947 } else if ( mark != largeInteger ) {
00948
00949 work[mark] = value;
00950 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00951 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00952
00953 temp2[word] = temp2[word] | ( 1 << bit );
00954 added--;
00955 } else {
00956 thisPivotValue = value;
00957 }
00958 }
00959
00960 elementU[put] = elementU[startColumn];
00961 indexRowU[put] = indexRowU[startColumn];
00962 if ( positionLargest == startColumn ) {
00963 positionLargest = put;
00964 }
00965 put++;
00966 elementU[startColumn] = thisPivotValue;
00967 indexRowU[startColumn] = pivotRow;
00968
00969 startColumn++;
00970 numberInColumn[iColumn] = put - startColumn;
00971 int * numberInColumnPlus = numberInColumnPlus_.array();
00972 numberInColumnPlus[iColumn]++;
00973 startColumnU[iColumn]++;
00974
00975 int next = nextColumn[iColumn];
00976 CoinBigIndex space;
00977
00978 space = startColumnU[next] - put - numberInColumnPlus[next];
00979
00980 if ( numberInPivotColumn > space ) {
00981
00982 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
00983 return false;
00984 }
00985
00986 if (positionLargest >= 0)
00987 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
00988 startColumn = startColumnU[iColumn];
00989 put = startColumn + numberInColumn[iColumn];
00990 }
00991 double tolerance = zeroTolerance_;
00992
00993 int *nextCount = nextCount_.array();
00994 for ( j = 0; j < numberInPivotColumn; j++ ) {
00995 value = work[j] - thisPivotValue * multipliersL[j];
00996 double absValue = fabs ( value );
00997
00998 if ( absValue > tolerance ) {
00999 work[j] = 0.0;
01000 assert (put<lengthAreaU_);
01001 elementU[put] = value;
01002 indexRowU[put] = indexL[j];
01003 if ( absValue > largest ) {
01004 largest = absValue;
01005 positionLargest = put;
01006 }
01007 put++;
01008 } else {
01009 work[j] = 0.0;
01010 added--;
01011 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01012 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01013
01014 if ( temp2[word] & ( 1 << bit ) ) {
01015
01016 iRow = indexL[j];
01017 CoinBigIndex start = startRowU[iRow];
01018 CoinBigIndex end = start + numberInRow[iRow];
01019 CoinBigIndex where = start;
01020
01021 while ( indexColumnU[where] != iColumn ) {
01022 where++;
01023 }
01024 #if DEBUG_COIN
01025 if ( where >= end ) {
01026 abort ( );
01027 }
01028 #endif
01029 indexColumnU[where] = indexColumnU[end - 1];
01030 numberInRow[iRow]--;
01031 } else {
01032
01033 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01034 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01035
01036 temp2[word] = temp2[word] | ( 1 << bit );
01037 }
01038 }
01039 }
01040 numberInColumn[iColumn] = put - startColumn;
01041
01042 if ( positionLargest >= 0 ) {
01043 value = elementU[positionLargest];
01044 iRow = indexRowU[positionLargest];
01045 elementU[positionLargest] = elementU[startColumn];
01046 indexRowU[positionLargest] = indexRowU[startColumn];
01047 elementU[startColumn] = value;
01048 indexRowU[startColumn] = iRow;
01049 }
01050
01051 if ( nextCount[iColumn + numberRows_] != -2 ) {
01052
01053 deleteLink ( iColumn + numberRows_ );
01054 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01055 }
01056 temp2 += increment2;
01057 }
01058
01059 unsigned int *putBase = workArea2;
01060 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01061 int i = 0;
01062
01063
01064 while ( bigLoops ) {
01065 bigLoops--;
01066 int bit;
01067 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01068 unsigned int *putThis = putBase;
01069 int iRow = indexL[i];
01070
01071
01072 int number = 0;
01073 int jColumn;
01074
01075 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01076 unsigned int test = *putThis;
01077
01078 putThis += increment2;
01079 test = 1 - ( ( test >> bit ) & 1 );
01080 number += test;
01081 }
01082 int next = nextRow[iRow];
01083 CoinBigIndex space;
01084
01085 space = startRowU[next] - startRowU[iRow];
01086 number += numberInRow[iRow];
01087 if ( space < number ) {
01088 if ( !getRowSpace ( iRow, number ) ) {
01089 return false;
01090 }
01091 }
01092
01093 putThis = putBase;
01094 next = nextRow[iRow];
01095 number = numberInRow[iRow];
01096 CoinBigIndex end = startRowU[iRow] + number;
01097 int saveIndex = indexColumnU[startRowU[next]];
01098
01099
01100 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01101 unsigned int test = *putThis;
01102
01103 putThis += increment2;
01104 test = 1 - ( ( test >> bit ) & 1 );
01105 indexColumnU[end] = saveColumn[jColumn];
01106 end += test;
01107 }
01108
01109 indexColumnU[startRowU[next]] = saveIndex;
01110 markRow[iRow] = static_cast<T>(largeInteger+1);
01111 number = end - startRowU[iRow];
01112 numberInRow[iRow] = number;
01113 deleteLink ( iRow );
01114 addLink ( iRow, number );
01115 }
01116 putBase++;
01117 }
01118 int bit;
01119
01120 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01121 unsigned int *putThis = putBase;
01122 int iRow = indexL[i];
01123
01124
01125 int number = 0;
01126 int jColumn;
01127
01128 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01129 unsigned int test = *putThis;
01130
01131 putThis += increment2;
01132 test = 1 - ( ( test >> bit ) & 1 );
01133 number += test;
01134 }
01135 int next = nextRow[iRow];
01136 CoinBigIndex space;
01137
01138 space = startRowU[next] - startRowU[iRow];
01139 number += numberInRow[iRow];
01140 if ( space < number ) {
01141 if ( !getRowSpace ( iRow, number ) ) {
01142 return false;
01143 }
01144 }
01145
01146 putThis = putBase;
01147 next = nextRow[iRow];
01148 number = numberInRow[iRow];
01149 CoinBigIndex end = startRowU[iRow] + number;
01150 int saveIndex;
01151
01152 saveIndex = indexColumnU[startRowU[next]];
01153
01154
01155 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01156 unsigned int test = *putThis;
01157
01158 putThis += increment2;
01159 test = 1 - ( ( test >> bit ) & 1 );
01160
01161 indexColumnU[end] = saveColumn[jColumn];
01162 end += test;
01163 }
01164 indexColumnU[startRowU[next]] = saveIndex;
01165 markRow[iRow] = static_cast<T>(largeInteger+1);
01166 number = end - startRowU[iRow];
01167 numberInRow[iRow] = number;
01168 deleteLink ( iRow );
01169 addLink ( iRow, number );
01170 }
01171 markRow[pivotRow] = static_cast<T>(largeInteger+1);
01172
01173 deleteLink ( pivotRow );
01174 deleteLink ( pivotColumn + numberRows_ );
01175 totalElements_ += added;
01176 return true;
01177 }
01178
01179
01181
01182 protected:
01183
01186
01187 double pivotTolerance_;
01189 double zeroTolerance_;
01190 #ifndef COIN_FAST_CODE
01192 double slackValue_;
01193 #else
01194 #ifndef slackValue_
01195 #define slackValue_ -1.0
01196 #endif
01197 #endif
01199 double areaFactor_;
01201 double relaxCheck_;
01203 int numberRows_;
01205 int numberRowsExtra_;
01207 int maximumRowsExtra_;
01209 int numberColumns_;
01211 int numberColumnsExtra_;
01213 int maximumColumnsExtra_;
01215 int numberGoodU_;
01217 int numberGoodL_;
01219 int maximumPivots_;
01221 int numberPivots_;
01224 CoinBigIndex totalElements_;
01226 CoinBigIndex factorElements_;
01228 CoinIntArrayWithLength pivotColumn_;
01230 CoinIntArrayWithLength permute_;
01232 CoinIntArrayWithLength permuteBack_;
01234 CoinIntArrayWithLength pivotColumnBack_;
01236 int status_;
01237
01242
01243
01245 int numberTrials_;
01247 CoinBigIndexArrayWithLength startRowU_;
01248
01250 CoinIntArrayWithLength numberInRow_;
01251
01253 CoinIntArrayWithLength numberInColumn_;
01254
01256 CoinIntArrayWithLength numberInColumnPlus_;
01257
01260 CoinIntArrayWithLength firstCount_;
01261
01263 CoinIntArrayWithLength nextCount_;
01264
01266 CoinIntArrayWithLength lastCount_;
01267
01269 CoinIntArrayWithLength nextColumn_;
01270
01272 CoinIntArrayWithLength lastColumn_;
01273
01275 CoinIntArrayWithLength nextRow_;
01276
01278 CoinIntArrayWithLength lastRow_;
01279
01281 CoinIntArrayWithLength saveColumn_;
01282
01284 CoinIntArrayWithLength markRow_;
01285
01287 int messageLevel_;
01288
01290 int biggerDimension_;
01291
01293 CoinIntArrayWithLength indexColumnU_;
01294
01296 CoinIntArrayWithLength pivotRowL_;
01297
01299 CoinFactorizationDoubleArrayWithLength pivotRegion_;
01300
01302 int numberSlacks_;
01303
01305 int numberU_;
01306
01308 CoinBigIndex maximumU_;
01309
01311
01312
01314 CoinBigIndex lengthU_;
01315
01317 CoinBigIndex lengthAreaU_;
01318
01320 CoinFactorizationDoubleArrayWithLength elementU_;
01321
01323 CoinIntArrayWithLength indexRowU_;
01324
01326 CoinBigIndexArrayWithLength startColumnU_;
01327
01329 CoinBigIndexArrayWithLength convertRowToColumnU_;
01330
01332 CoinBigIndex numberL_;
01333
01335 CoinBigIndex baseL_;
01336
01338 CoinBigIndex lengthL_;
01339
01341 CoinBigIndex lengthAreaL_;
01342
01344 CoinFactorizationDoubleArrayWithLength elementL_;
01345
01347 CoinIntArrayWithLength indexRowL_;
01348
01350 CoinBigIndexArrayWithLength startColumnL_;
01351
01353 bool doForrestTomlin_;
01354
01356 int numberR_;
01357
01359 CoinBigIndex lengthR_;
01360
01362 CoinBigIndex lengthAreaR_;
01363
01365 CoinFactorizationDouble *elementR_;
01366
01368 int *indexRowR_;
01369
01371 CoinBigIndexArrayWithLength startColumnR_;
01372
01374 double * denseArea_;
01375
01377 int * densePermute_;
01378
01380 int numberDense_;
01381
01383 int denseThreshold_;
01384
01386 CoinFactorizationDoubleArrayWithLength workArea_;
01387
01389 CoinUnsignedIntArrayWithLength workArea2_;
01390
01392 CoinBigIndex numberCompressions_;
01393
01395 mutable double ftranCountInput_;
01396 mutable double ftranCountAfterL_;
01397 mutable double ftranCountAfterR_;
01398 mutable double ftranCountAfterU_;
01399 mutable double btranCountInput_;
01400 mutable double btranCountAfterU_;
01401 mutable double btranCountAfterR_;
01402 mutable double btranCountAfterL_;
01403
01405 mutable int numberFtranCounts_;
01406 mutable int numberBtranCounts_;
01407
01409 double ftranAverageAfterL_;
01410 double ftranAverageAfterR_;
01411 double ftranAverageAfterU_;
01412 double btranAverageAfterU_;
01413 double btranAverageAfterR_;
01414 double btranAverageAfterL_;
01415
01417 mutable bool collectStatistics_;
01418
01420 int sparseThreshold_;
01421
01423 int sparseThreshold2_;
01424
01426 CoinBigIndexArrayWithLength startRowL_;
01427
01429 CoinIntArrayWithLength indexColumnL_;
01430
01432 CoinFactorizationDoubleArrayWithLength elementByRowL_;
01433
01435 mutable CoinIntArrayWithLength sparse_;
01439 int biasLU_;
01445 int persistenceFlag_;
01447 };
01448
01449 #ifdef COIN_HAS_LAPACK
01450 #define DENSE_CODE 1
01451
01452 #ifndef ipfint
01453
01454 typedef int ipfint;
01455 typedef const int cipfint;
01456 #endif
01457 #endif
01458 #endif
01459
01460 #ifdef UGLY_COIN_FACTOR_CODING
01461 #define FAC_UNSET (FAC_SET+1)
01462 {
01463 goodPivot=false;
01464
01465 CoinBigIndex startColumnThis = startColumn[iPivotColumn];
01466 CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
01467 int put = 0;
01468 CoinBigIndex startRowThis = startRow[iPivotRow];
01469 CoinBigIndex endRow = startRowThis + numberDoRow + 1;
01470 if ( pivotColumnPosition < 0 ) {
01471 for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01472 int iColumn = indexColumn[pivotColumnPosition];
01473 if ( iColumn != iPivotColumn ) {
01474 saveColumn[put++] = iColumn;
01475 } else {
01476 break;
01477 }
01478 }
01479 } else {
01480 for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
01481 saveColumn[put++] = indexColumn[i];
01482 }
01483 }
01484 assert (pivotColumnPosition<endRow);
01485 assert (indexColumn[pivotColumnPosition]==iPivotColumn);
01486 pivotColumnPosition++;
01487 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01488 saveColumn[put++] = indexColumn[pivotColumnPosition];
01489 }
01490
01491 int next = nextRow[iPivotRow];
01492 int last = lastRow[iPivotRow];
01493
01494 nextRow[last] = next;
01495 lastRow[next] = last;
01496 nextRow[iPivotRow] = numberGoodU_;
01497 lastRow[iPivotRow] = -2;
01498 numberInRow[iPivotRow] = 0;
01499
01500 CoinBigIndex l = lengthL_;
01501
01502 {
01503 if ( l + numberDoColumn > lengthAreaL_ ) {
01504
01505 if ((messageLevel_&4)!=0)
01506 printf("more memory needed in middle of invert\n");
01507 goto BAD_PIVOT;
01508 }
01509
01510 CoinBigIndex lSave = l;
01511
01512 CoinBigIndex * startColumnL = startColumnL_.array();
01513 startColumnL[numberGoodL_] = l;
01514 numberGoodL_++;
01515 startColumnL[numberGoodL_] = l + numberDoColumn;
01516 lengthL_ += numberDoColumn;
01517 if ( pivotRowPosition < 0 ) {
01518 for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01519 int iRow = indexRow[pivotRowPosition];
01520 if ( iRow != iPivotRow ) {
01521 indexRowL[l] = iRow;
01522 elementL[l] = element[pivotRowPosition];
01523 markRow[iRow] = l - lSave;
01524 l++;
01525
01526 CoinBigIndex start = startRow[iRow];
01527 CoinBigIndex end = start + numberInRow[iRow];
01528 CoinBigIndex where = start;
01529
01530 while ( indexColumn[where] != iPivotColumn ) {
01531 where++;
01532 }
01533 #if DEBUG_COIN
01534 if ( where >= end ) {
01535 abort ( );
01536 }
01537 #endif
01538 indexColumn[where] = indexColumn[end - 1];
01539 numberInRow[iRow]--;
01540 } else {
01541 break;
01542 }
01543 }
01544 } else {
01545 CoinBigIndex i;
01546
01547 for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
01548 int iRow = indexRow[i];
01549
01550 markRow[iRow] = l - lSave;
01551 indexRowL[l] = iRow;
01552 elementL[l] = element[i];
01553 l++;
01554
01555 CoinBigIndex start = startRow[iRow];
01556 CoinBigIndex end = start + numberInRow[iRow];
01557 CoinBigIndex where = start;
01558
01559 while ( indexColumn[where] != iPivotColumn ) {
01560 where++;
01561 }
01562 #if DEBUG_COIN
01563 if ( where >= end ) {
01564 abort ( );
01565 }
01566 #endif
01567 indexColumn[where] = indexColumn[end - 1];
01568 numberInRow[iRow]--;
01569 assert (numberInRow[iRow]>=0);
01570 }
01571 }
01572 assert (pivotRowPosition<endColumn);
01573 assert (indexRow[pivotRowPosition]==iPivotRow);
01574 CoinFactorizationDouble pivotElement = element[pivotRowPosition];
01575 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
01576
01577 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
01578 pivotRowPosition++;
01579 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01580 int iRow = indexRow[pivotRowPosition];
01581
01582 markRow[iRow] = l - lSave;
01583 indexRowL[l] = iRow;
01584 elementL[l] = element[pivotRowPosition];
01585 l++;
01586
01587 CoinBigIndex start = startRow[iRow];
01588 CoinBigIndex end = start + numberInRow[iRow];
01589 CoinBigIndex where = start;
01590
01591 while ( indexColumn[where] != iPivotColumn ) {
01592 where++;
01593 }
01594 #if DEBUG_COIN
01595 if ( where >= end ) {
01596 abort ( );
01597 }
01598 #endif
01599 indexColumn[where] = indexColumn[end - 1];
01600 numberInRow[iRow]--;
01601 assert (numberInRow[iRow]>=0);
01602 }
01603 markRow[iPivotRow] = FAC_SET;
01604
01605 numberInColumn[iPivotColumn] = 0;
01606
01607 int *indexL = &indexRowL[lSave];
01608 CoinFactorizationDouble *multipliersL = &elementL[lSave];
01609
01610
01611 int j;
01612
01613 for ( j = 0; j < numberDoColumn; j++ ) {
01614 multipliersL[j] *= pivotMultiplier;
01615 }
01616
01617 CoinBigIndex iErase;
01618 for ( iErase = 0; iErase < increment2 * numberDoRow;
01619 iErase++ ) {
01620 workArea2[iErase] = 0;
01621 }
01622 CoinBigIndex added = numberDoRow * numberDoColumn;
01623 unsigned int *temp2 = workArea2;
01624 int * nextColumn = nextColumn_.array();
01625
01626
01627 int jColumn;
01628 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01629 int iColumn = saveColumn[jColumn];
01630 CoinBigIndex startColumnThis = startColumn[iColumn];
01631 CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
01632 int iRow = indexRow[startColumnThis];
01633 CoinFactorizationDouble value = element[startColumnThis];
01634 double largest;
01635 CoinBigIndex put = startColumnThis;
01636 CoinBigIndex positionLargest = -1;
01637 CoinFactorizationDouble thisPivotValue = 0.0;
01638
01639
01640 bool checkLargest;
01641 int mark = markRow[iRow];
01642
01643 if ( mark == FAC_UNSET ) {
01644 largest = fabs ( value );
01645 positionLargest = put;
01646 put++;
01647 checkLargest = false;
01648 } else {
01649
01650 largest = 0.0;
01651 checkLargest = true;
01652 if ( mark != FAC_SET ) {
01653
01654 workArea[mark] = value;
01655 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01656 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01657
01658 temp2[word] = temp2[word] | ( 1 << bit );
01659 added--;
01660 } else {
01661 thisPivotValue = value;
01662 }
01663 }
01664 CoinBigIndex i;
01665 for ( i = startColumnThis + 1; i < endColumn; i++ ) {
01666 iRow = indexRow[i];
01667 value = element[i];
01668 int mark = markRow[iRow];
01669
01670 if ( mark == FAC_UNSET ) {
01671
01672 indexRow[put] = iRow;
01673 element[put] = value;
01674 if ( checkLargest ) {
01675 double absValue = fabs ( value );
01676
01677 if ( absValue > largest ) {
01678 largest = absValue;
01679 positionLargest = put;
01680 }
01681 }
01682 put++;
01683 } else if ( mark != FAC_SET ) {
01684
01685 workArea[mark] = value;
01686 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01687 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01688
01689 temp2[word] = temp2[word] | ( 1 << bit );
01690 added--;
01691 } else {
01692 thisPivotValue = value;
01693 }
01694 }
01695
01696 element[put] = element[startColumnThis];
01697 indexRow[put] = indexRow[startColumnThis];
01698 if ( positionLargest == startColumnThis ) {
01699 positionLargest = put;
01700 }
01701 put++;
01702 element[startColumnThis] = thisPivotValue;
01703 indexRow[startColumnThis] = iPivotRow;
01704
01705 startColumnThis++;
01706 numberInColumn[iColumn] = put - startColumnThis;
01707 int * numberInColumnPlus = numberInColumnPlus_.array();
01708 numberInColumnPlus[iColumn]++;
01709 startColumn[iColumn]++;
01710
01711 int next = nextColumn[iColumn];
01712 CoinBigIndex space;
01713
01714 space = startColumn[next] - put - numberInColumnPlus[next];
01715
01716 if ( numberDoColumn > space ) {
01717
01718 if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
01719 goto BAD_PIVOT;
01720 }
01721
01722 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
01723 startColumnThis = startColumn[iColumn];
01724 put = startColumnThis + numberInColumn[iColumn];
01725 }
01726 double tolerance = zeroTolerance_;
01727
01728 int *nextCount = nextCount_.array();
01729 for ( j = 0; j < numberDoColumn; j++ ) {
01730 value = workArea[j] - thisPivotValue * multipliersL[j];
01731 double absValue = fabs ( value );
01732
01733 if ( absValue > tolerance ) {
01734 workArea[j] = 0.0;
01735 element[put] = value;
01736 indexRow[put] = indexL[j];
01737 if ( absValue > largest ) {
01738 largest = absValue;
01739 positionLargest = put;
01740 }
01741 put++;
01742 } else {
01743 workArea[j] = 0.0;
01744 added--;
01745 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01746 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01747
01748 if ( temp2[word] & ( 1 << bit ) ) {
01749
01750 iRow = indexL[j];
01751 CoinBigIndex start = startRow[iRow];
01752 CoinBigIndex end = start + numberInRow[iRow];
01753 CoinBigIndex where = start;
01754
01755 while ( indexColumn[where] != iColumn ) {
01756 where++;
01757 }
01758 #if DEBUG_COIN
01759 if ( where >= end ) {
01760 abort ( );
01761 }
01762 #endif
01763 indexColumn[where] = indexColumn[end - 1];
01764 numberInRow[iRow]--;
01765 } else {
01766
01767 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01768 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01769
01770 temp2[word] = temp2[word] | ( 1 << bit );
01771 }
01772 }
01773 }
01774 numberInColumn[iColumn] = put - startColumnThis;
01775
01776 if ( positionLargest >= 0 ) {
01777 value = element[positionLargest];
01778 iRow = indexRow[positionLargest];
01779 element[positionLargest] = element[startColumnThis];
01780 indexRow[positionLargest] = indexRow[startColumnThis];
01781 element[startColumnThis] = value;
01782 indexRow[startColumnThis] = iRow;
01783 }
01784
01785 if ( nextCount[iColumn + numberRows_] != -2 ) {
01786
01787 deleteLink ( iColumn + numberRows_ );
01788 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01789 }
01790 temp2 += increment2;
01791 }
01792
01793 unsigned int *putBase = workArea2;
01794 int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01795 int i = 0;
01796
01797
01798 while ( bigLoops ) {
01799 bigLoops--;
01800 int bit;
01801 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01802 unsigned int *putThis = putBase;
01803 int iRow = indexL[i];
01804
01805
01806 int number = 0;
01807 int jColumn;
01808
01809 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01810 unsigned int test = *putThis;
01811
01812 putThis += increment2;
01813 test = 1 - ( ( test >> bit ) & 1 );
01814 number += test;
01815 }
01816 int next = nextRow[iRow];
01817 CoinBigIndex space;
01818
01819 space = startRow[next] - startRow[iRow];
01820 number += numberInRow[iRow];
01821 if ( space < number ) {
01822 if ( !getRowSpace ( iRow, number ) ) {
01823 goto BAD_PIVOT;
01824 }
01825 }
01826
01827 putThis = putBase;
01828 next = nextRow[iRow];
01829 number = numberInRow[iRow];
01830 CoinBigIndex end = startRow[iRow] + number;
01831 int saveIndex = indexColumn[startRow[next]];
01832
01833
01834 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01835 unsigned int test = *putThis;
01836
01837 putThis += increment2;
01838 test = 1 - ( ( test >> bit ) & 1 );
01839 indexColumn[end] = saveColumn[jColumn];
01840 end += test;
01841 }
01842
01843 indexColumn[startRow[next]] = saveIndex;
01844 markRow[iRow] = FAC_UNSET;
01845 number = end - startRow[iRow];
01846 numberInRow[iRow] = number;
01847 deleteLink ( iRow );
01848 addLink ( iRow, number );
01849 }
01850 putBase++;
01851 }
01852 int bit;
01853
01854 for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
01855 unsigned int *putThis = putBase;
01856 int iRow = indexL[i];
01857
01858
01859 int number = 0;
01860 int jColumn;
01861
01862 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01863 unsigned int test = *putThis;
01864
01865 putThis += increment2;
01866 test = 1 - ( ( test >> bit ) & 1 );
01867 number += test;
01868 }
01869 int next = nextRow[iRow];
01870 CoinBigIndex space;
01871
01872 space = startRow[next] - startRow[iRow];
01873 number += numberInRow[iRow];
01874 if ( space < number ) {
01875 if ( !getRowSpace ( iRow, number ) ) {
01876 goto BAD_PIVOT;
01877 }
01878 }
01879
01880 putThis = putBase;
01881 next = nextRow[iRow];
01882 number = numberInRow[iRow];
01883 CoinBigIndex end = startRow[iRow] + number;
01884 int saveIndex;
01885
01886 saveIndex = indexColumn[startRow[next]];
01887
01888
01889 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01890 unsigned int test = *putThis;
01891
01892 putThis += increment2;
01893 test = 1 - ( ( test >> bit ) & 1 );
01894
01895 indexColumn[end] = saveColumn[jColumn];
01896 end += test;
01897 }
01898 indexColumn[startRow[next]] = saveIndex;
01899 markRow[iRow] = FAC_UNSET;
01900 number = end - startRow[iRow];
01901 numberInRow[iRow] = number;
01902 deleteLink ( iRow );
01903 addLink ( iRow, number );
01904 }
01905 markRow[iPivotRow] = FAC_UNSET;
01906
01907 deleteLink ( iPivotRow );
01908 deleteLink ( iPivotColumn + numberRows_ );
01909 totalElements_ += added;
01910 goodPivot= true;
01911
01912 }
01913 BAD_PIVOT:
01914
01915 ;
01916 }
01917 #undef FAC_UNSET
01918 #endif