00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #ifndef CoinFactorization_H
00011 #define CoinFactorization_H
00012
00013 #include <iostream>
00014 #include <string>
00015 #include <cassert>
00016 #include <cstdio>
00017 #include "CoinFinite.hpp"
00018 #include "CoinIndexedVector.hpp"
00019 class CoinPackedMatrix;
00046 class CoinFactorization {
00047 friend void CoinFactorizationUnitTest( const std::string & mpsDir );
00048
00049 public:
00050
00053
00054 CoinFactorization ( );
00056 CoinFactorization ( const CoinFactorization &other);
00057
00059 ~CoinFactorization ( );
00061 void almostDestructor();
00063 void show_self ( ) const;
00065 int saveFactorization (const char * file ) const;
00069 int restoreFactorization (const char * file , bool factor=false) ;
00071 void sort ( ) const;
00073 CoinFactorization & operator = ( const CoinFactorization & other );
00075
00085 int factorize ( const CoinPackedMatrix & matrix,
00086 int rowIsBasic[], int columnIsBasic[] ,
00087 double areaFactor = 0.0 );
00098 int factorize ( int numberRows,
00099 int numberColumns,
00100 CoinBigIndex numberElements,
00101 CoinBigIndex maximumL,
00102 CoinBigIndex maximumU,
00103 const int indicesRow[],
00104 const int indicesColumn[], const double elements[] ,
00105 int permutation[],
00106 double areaFactor = 0.0);
00111 int factorizePart1 ( int numberRows,
00112 int numberColumns,
00113 CoinBigIndex estimateNumberElements,
00114 int * indicesRow[],
00115 int * indicesColumn[],
00116 double * elements[],
00117 double areaFactor = 0.0);
00124 int factorizePart2 (int permutation[],int exactNumberElements);
00126 double conditionNumber() const;
00127
00129
00132
00133 inline int status ( ) const {
00134 return status_;
00135 }
00137 inline void setStatus ( int value)
00138 { status_=value; }
00140 inline int pivots ( ) const {
00141 return numberPivots_;
00142 }
00144 inline void setPivots ( int value )
00145 { numberPivots_=value; }
00147 inline int *permute ( ) const {
00148 return permute_.array();
00149 }
00151 inline int *pivotColumn ( ) const {
00152 return pivotColumn_.array();
00153 }
00155 inline double *pivotRegion ( ) const {
00156 return pivotRegion_.array();
00157 }
00159 inline int *permuteBack ( ) const {
00160 return permuteBack_.array();
00161 }
00163 inline int *pivotColumnBack ( ) const {
00164 return pivotColumnBack_.array();
00165 }
00167 inline CoinBigIndex * startRowL() const
00168 { return startRowL_.array();}
00169
00171 inline CoinBigIndex * startColumnL() const
00172 { return startColumnL_.array();}
00173
00175 inline int * indexColumnL() const
00176 { return indexColumnL_.array();}
00177
00179 inline int * indexRowL() const
00180 { return indexRowL_.array();}
00181
00183 inline double * elementByRowL() const
00184 { return elementByRowL_.array();}
00185
00187 inline int numberRowsExtra ( ) const {
00188 return numberRowsExtra_;
00189 }
00191 inline void setNumberRows(int value)
00192 { numberRows_ = value; }
00194 inline int numberRows ( ) const {
00195 return numberRows_;
00196 }
00198 inline CoinBigIndex numberL() const
00199 { return numberL_;}
00200
00202 inline CoinBigIndex baseL() const
00203 { return baseL_;}
00205 inline int maximumRowsExtra ( ) const {
00206 return maximumRowsExtra_;
00207 }
00209 inline int numberColumns ( ) const {
00210 return numberColumns_;
00211 }
00213 inline int numberElements ( ) const {
00214 return totalElements_;
00215 }
00217 inline int numberForrestTomlin ( ) const {
00218 return numberInColumn_.array()[numberColumnsExtra_];
00219 }
00221 inline int numberGoodColumns ( ) const {
00222 return numberGoodU_;
00223 }
00225 inline double areaFactor ( ) const {
00226 return areaFactor_;
00227 }
00228 inline void areaFactor ( double value ) {
00229 areaFactor_=value;
00230 }
00232 double adjustedAreaFactor() const;
00234 inline void relaxAccuracyCheck(double value)
00235 { relaxCheck_ = value;}
00236 inline double getAccuracyCheck() const
00237 { return relaxCheck_;}
00239 inline int messageLevel ( ) const {
00240 return messageLevel_ ;
00241 }
00242 void messageLevel ( int value );
00244 inline int maximumPivots ( ) const {
00245 return maximumPivots_ ;
00246 }
00247 void maximumPivots ( int value );
00248
00250 inline int denseThreshold() const
00251 { return denseThreshold_;}
00253 inline void setDenseThreshold(int value)
00254 { denseThreshold_ = value;}
00256 inline double pivotTolerance ( ) const {
00257 return pivotTolerance_ ;
00258 }
00259 void pivotTolerance ( double value );
00261 inline double zeroTolerance ( ) const {
00262 return zeroTolerance_ ;
00263 }
00264 void zeroTolerance ( double value );
00265 #ifndef COIN_FAST_CODE
00267 inline double slackValue ( ) const {
00268 return slackValue_ ;
00269 }
00270 void slackValue ( double value );
00271 #endif
00273 double maximumCoefficient() const;
00275 inline bool forrestTomlin() const
00276 { return doForrestTomlin_;}
00277 inline void setForrestTomlin(bool value)
00278 { doForrestTomlin_=value;}
00280 inline bool spaceForForrestTomlin() const
00281 {
00282 CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
00283 CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
00284 return (space>=0)&&doForrestTomlin_;
00285 }
00287
00290
00292 inline int numberDense() const
00293 { return numberDense_;}
00294
00296 inline CoinBigIndex numberElementsU ( ) const {
00297 return lengthU_;
00298 }
00300 inline void setNumberElementsU(CoinBigIndex value)
00301 { lengthU_ = value; }
00303 inline CoinBigIndex lengthAreaU ( ) const {
00304 return lengthAreaU_;
00305 }
00307 inline CoinBigIndex numberElementsL ( ) const {
00308 return lengthL_;
00309 }
00311 inline CoinBigIndex lengthAreaL ( ) const {
00312 return lengthAreaL_;
00313 }
00315 inline CoinBigIndex numberElementsR ( ) const {
00316 return lengthR_;
00317 }
00319 inline CoinBigIndex numberCompressions() const
00320 { return numberCompressions_;}
00322 inline int * numberInRow() const
00323 { return numberInRow_.array();}
00325 inline int * numberInColumn() const
00326 { return numberInColumn_.array();}
00328 inline double * elementU() const
00329 { return elementU_.array();}
00331 inline int * indexRowU() const
00332 { return indexRowU_.array();}
00334 inline CoinBigIndex * startColumnU() const
00335 { return startColumnU_.array();}
00337 inline int maximumColumnsExtra()
00338 { return maximumColumnsExtra_;}
00342 inline int biasLU() const
00343 { return biasLU_;}
00344 inline void setBiasLU(int value)
00345 { biasLU_=value;}
00351 inline int persistenceFlag() const
00352 { return persistenceFlag_;}
00353 void setPersistenceFlag(int value);
00355
00358
00366 int replaceColumn ( CoinIndexedVector * regionSparse,
00367 int pivotRow,
00368 double pivotCheck ,
00369 bool checkBeforeModifying=false);
00371
00381 int updateColumnFT ( CoinIndexedVector * regionSparse,
00382 CoinIndexedVector * regionSparse2);
00385 int updateColumn ( CoinIndexedVector * regionSparse,
00386 CoinIndexedVector * regionSparse2,
00387 bool noPermute=false) const;
00393 int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
00394 CoinIndexedVector * regionSparse2,
00395 CoinIndexedVector * regionSparse3,
00396 bool noPermuteRegion3=false) ;
00401 int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00402 CoinIndexedVector * regionSparse2) const;
00404 void goSparse();
00406 inline int sparseThreshold ( ) const
00407 { return sparseThreshold_;}
00409 void sparseThreshold ( int value );
00411
00412
00416
00417 inline void clearArrays()
00418 { gutsOfDestructor();}
00420
00423
00426 int add ( CoinBigIndex numberElements,
00427 int indicesRow[],
00428 int indicesColumn[], double elements[] );
00429
00432 int addColumn ( CoinBigIndex numberElements,
00433 int indicesRow[], double elements[] );
00434
00437 int addRow ( CoinBigIndex numberElements,
00438 int indicesColumn[], double elements[] );
00439
00441 int deleteColumn ( int Row );
00443 int deleteRow ( int Row );
00444
00448 int replaceRow ( int whichRow, int numberElements,
00449 const int indicesColumn[], const double elements[] );
00451 void emptyRows(int numberToEmpty, const int which[]);
00453
00454
00455 void checkSparse();
00457 inline bool collectStatistics() const
00458 { return collectStatistics_;}
00460 inline void setCollectStatistics(bool onOff) const
00461 { collectStatistics_ = onOff;}
00463 void gutsOfDestructor(int type=1);
00465 void gutsOfInitialize(int type);
00466 void gutsOfCopy(const CoinFactorization &other);
00467
00469 void resetStatistics();
00470
00471
00473
00475
00476 void getAreas ( int numberRows,
00477 int numberColumns,
00478 CoinBigIndex maximumL,
00479 CoinBigIndex maximumU );
00480
00483 void preProcess ( int state,
00484 int possibleDuplicates = -1 );
00486 int factor ( );
00487 protected:
00490 int factorSparse ( );
00493 int factorSparseSmall ( );
00496 int factorSparseLarge ( );
00499 int factorDense ( );
00500
00502 bool pivotOneOtherRow ( int pivotRow,
00503 int pivotColumn );
00505 bool pivotRowSingleton ( int pivotRow,
00506 int pivotColumn );
00508 bool pivotColumnSingleton ( int pivotRow,
00509 int pivotColumn );
00510
00515 bool getColumnSpace ( int iColumn,
00516 int extraNeeded );
00517
00521 bool getColumnSpaceIterateR ( int iColumn, double value,
00522 int iRow);
00528 CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00529 int iRow);
00533 bool getRowSpace ( int iRow, int extraNeeded );
00534
00538 bool getRowSpaceIterate ( int iRow,
00539 int extraNeeded );
00541 void checkConsistency ( );
00543 inline void addLink ( int index, int count ) {
00544 int *nextCount = nextCount_.array();
00545 int *firstCount = firstCount_.array();
00546 int *lastCount = lastCount_.array();
00547 int next = firstCount[count];
00548 lastCount[index] = -2 - count;
00549 if ( next < 0 ) {
00550
00551 firstCount[count] = index;
00552 nextCount[index] = -1;
00553 } else {
00554 firstCount[count] = index;
00555 nextCount[index] = next;
00556 lastCount[next] = index;
00557 }}
00559 inline void deleteLink ( int index ) {
00560 int *nextCount = nextCount_.array();
00561 int *firstCount = firstCount_.array();
00562 int *lastCount = lastCount_.array();
00563 int next = nextCount[index];
00564 int last = lastCount[index];
00565 if ( last >= 0 ) {
00566 nextCount[last] = next;
00567 } else {
00568 int count = -last - 2;
00569
00570 firstCount[count] = next;
00571 }
00572 if ( next >= 0 ) {
00573 lastCount[next] = last;
00574 }
00575 nextCount[index] = -2;
00576 lastCount[index] = -2;
00577 return;
00578 }
00580 void separateLinks(int count,bool rowsFirst);
00582 void cleanup ( );
00583
00585 void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00587 void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00589 void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00591 void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00592
00594 void updateColumnR ( CoinIndexedVector * region ) const;
00597 void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00598
00600 void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00601
00603 void updateColumnUSparse ( CoinIndexedVector * regionSparse,
00604 int * indexIn) const;
00606 void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
00607 int * indexIn) const;
00609 int updateColumnUDensish ( double * COIN_RESTRICT region,
00610 int * COIN_RESTRICT regionIndex) const;
00612 void updateTwoColumnsUDensish (
00613 int & numberNonZero1,
00614 double * COIN_RESTRICT region1,
00615 int * COIN_RESTRICT index1,
00616 int & numberNonZero2,
00617 double * COIN_RESTRICT region2,
00618 int * COIN_RESTRICT index2) const;
00620 void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
00622 void permuteBack ( CoinIndexedVector * regionSparse,
00623 CoinIndexedVector * outVector) const;
00624
00626 void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00629 void updateColumnTransposeU ( CoinIndexedVector * region,
00630 int smallestIndex) const;
00633 void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00634 int smallestIndex) const;
00637 void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00638 int smallestIndex) const;
00641 void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00642
00644 void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00646 void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00648 void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00649
00651 void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00653 void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00655 void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00657 void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00659 void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00660 public:
00665 int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00666 int pivotRow, double alpha);
00667 protected:
00670 int checkPivot(double saveFromU, double oldPivot) const;
00671
00672 #ifdef INT_IS_8
00673 #define COINFACTORIZATION_BITS_PER_INT 64
00674 #define COINFACTORIZATION_SHIFT_PER_INT 6
00675 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00676 #else
00677 #define COINFACTORIZATION_BITS_PER_INT 32
00678 #define COINFACTORIZATION_SHIFT_PER_INT 5
00679 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00680 #endif
00681 template <class T> inline bool
00682 pivot ( int pivotRow,
00683 int pivotColumn,
00684 CoinBigIndex pivotRowPosition,
00685 CoinBigIndex pivotColumnPosition,
00686 double work[],
00687 unsigned int workArea2[],
00688 int increment,
00689 int increment2,
00690 T markRow[] ,
00691 int largeInteger)
00692 {
00693 int *indexColumnU = indexColumnU_.array();
00694 CoinBigIndex *startColumnU = startColumnU_.array();
00695 int *numberInColumn = numberInColumn_.array();
00696 double *elementU = elementU_.array();
00697 int *indexRowU = indexRowU_.array();
00698 CoinBigIndex *startRowU = startRowU_.array();
00699 int *numberInRow = numberInRow_.array();
00700 double *elementL = elementL_.array();
00701 int *indexRowL = indexRowL_.array();
00702 int *saveColumn = saveColumn_.array();
00703 int *nextRow = nextRow_.array();
00704 int *lastRow = lastRow_.array() ;
00705
00706
00707 int numberInPivotRow = numberInRow[pivotRow] - 1;
00708 CoinBigIndex startColumn = startColumnU[pivotColumn];
00709 int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00710 CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00711 int put = 0;
00712 CoinBigIndex startRow = startRowU[pivotRow];
00713 CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00714
00715 if ( pivotColumnPosition < 0 ) {
00716 for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00717 int iColumn = indexColumnU[pivotColumnPosition];
00718 if ( iColumn != pivotColumn ) {
00719 saveColumn[put++] = iColumn;
00720 } else {
00721 break;
00722 }
00723 }
00724 } else {
00725 for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00726 saveColumn[put++] = indexColumnU[i];
00727 }
00728 }
00729 assert (pivotColumnPosition<endRow);
00730 assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00731 pivotColumnPosition++;
00732 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00733 saveColumn[put++] = indexColumnU[pivotColumnPosition];
00734 }
00735
00736 int next = nextRow[pivotRow];
00737 int last = lastRow[pivotRow];
00738
00739 nextRow[last] = next;
00740 lastRow[next] = last;
00741 nextRow[pivotRow] = numberGoodU_;
00742 lastRow[pivotRow] = -2;
00743 numberInRow[pivotRow] = 0;
00744
00745 CoinBigIndex l = lengthL_;
00746
00747 if ( l + numberInPivotColumn > lengthAreaL_ ) {
00748
00749 printf("more memory needed in middle of invert\n");
00750 return false;
00751 }
00752
00753 CoinBigIndex lSave = l;
00754
00755 pivotRowL_.array()[numberGoodL_] = pivotRow;
00756 CoinBigIndex * startColumnL = startColumnL_.array();
00757 startColumnL[numberGoodL_] = l;
00758 numberGoodL_++;
00759 startColumnL[numberGoodL_] = l + numberInPivotColumn;
00760 lengthL_ += numberInPivotColumn;
00761 if ( pivotRowPosition < 0 ) {
00762 for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00763 int iRow = indexRowU[pivotRowPosition];
00764 if ( iRow != pivotRow ) {
00765 indexRowL[l] = iRow;
00766 elementL[l] = elementU[pivotRowPosition];
00767 markRow[iRow] = static_cast<T>(l - lSave);
00768 l++;
00769
00770 CoinBigIndex start = startRowU[iRow];
00771 CoinBigIndex end = start + numberInRow[iRow];
00772 CoinBigIndex where = start;
00773
00774 while ( indexColumnU[where] != pivotColumn ) {
00775 where++;
00776 }
00777 #if DEBUG_COIN
00778 if ( where >= end ) {
00779 abort ( );
00780 }
00781 #endif
00782 indexColumnU[where] = indexColumnU[end - 1];
00783 numberInRow[iRow]--;
00784 } else {
00785 break;
00786 }
00787 }
00788 } else {
00789 CoinBigIndex i;
00790
00791 for ( i = startColumn; i < pivotRowPosition; i++ ) {
00792 int iRow = indexRowU[i];
00793
00794 markRow[iRow] = static_cast<T>(l - lSave);
00795 indexRowL[l] = iRow;
00796 elementL[l] = elementU[i];
00797 l++;
00798
00799 CoinBigIndex start = startRowU[iRow];
00800 CoinBigIndex end = start + numberInRow[iRow];
00801 CoinBigIndex where = start;
00802
00803 while ( indexColumnU[where] != pivotColumn ) {
00804 where++;
00805 }
00806 #if DEBUG_COIN
00807 if ( where >= end ) {
00808 abort ( );
00809 }
00810 #endif
00811 indexColumnU[where] = indexColumnU[end - 1];
00812 numberInRow[iRow]--;
00813 assert (numberInRow[iRow]>=0);
00814 }
00815 }
00816 assert (pivotRowPosition<endColumn);
00817 assert (indexRowU[pivotRowPosition]==pivotRow);
00818 double pivotElement = elementU[pivotRowPosition];
00819 double pivotMultiplier = 1.0 / pivotElement;
00820
00821 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00822 pivotRowPosition++;
00823 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00824 int iRow = indexRowU[pivotRowPosition];
00825
00826 markRow[iRow] = static_cast<T>(l - lSave);
00827 indexRowL[l] = iRow;
00828 elementL[l] = elementU[pivotRowPosition];
00829 l++;
00830
00831 CoinBigIndex start = startRowU[iRow];
00832 CoinBigIndex end = start + numberInRow[iRow];
00833 CoinBigIndex where = start;
00834
00835 while ( indexColumnU[where] != pivotColumn ) {
00836 where++;
00837 }
00838 #if DEBUG_COIN
00839 if ( where >= end ) {
00840 abort ( );
00841 }
00842 #endif
00843 indexColumnU[where] = indexColumnU[end - 1];
00844 numberInRow[iRow]--;
00845 assert (numberInRow[iRow]>=0);
00846 }
00847 markRow[pivotRow] = static_cast<T>(largeInteger);
00848
00849 numberInColumn[pivotColumn] = 0;
00850
00851 int *indexL = &indexRowL[lSave];
00852 double *multipliersL = &elementL[lSave];
00853
00854
00855 int j;
00856
00857 for ( j = 0; j < numberInPivotColumn; j++ ) {
00858 multipliersL[j] *= pivotMultiplier;
00859 }
00860
00861 CoinBigIndex iErase;
00862 for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00863 iErase++ ) {
00864 workArea2[iErase] = 0;
00865 }
00866 CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00867 unsigned int *temp2 = workArea2;
00868 int * nextColumn = nextColumn_.array();
00869
00870
00871 int jColumn;
00872 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
00873 int iColumn = saveColumn[jColumn];
00874 CoinBigIndex startColumn = startColumnU[iColumn];
00875 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
00876 int iRow = indexRowU[startColumn];
00877 double value = elementU[startColumn];
00878 double largest;
00879 CoinBigIndex put = startColumn;
00880 CoinBigIndex positionLargest = -1;
00881 double thisPivotValue = 0.0;
00882
00883
00884 bool checkLargest;
00885 int mark = markRow[iRow];
00886
00887 if ( mark == largeInteger+1 ) {
00888 largest = fabs ( value );
00889 positionLargest = put;
00890 put++;
00891 checkLargest = false;
00892 } else {
00893
00894 largest = 0.0;
00895 checkLargest = true;
00896 if ( mark != largeInteger ) {
00897
00898 work[mark] = value;
00899 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00900 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00901
00902 temp2[word] = temp2[word] | ( 1 << bit );
00903 added--;
00904 } else {
00905 thisPivotValue = value;
00906 }
00907 }
00908 CoinBigIndex i;
00909 for ( i = startColumn + 1; i < endColumn; i++ ) {
00910 iRow = indexRowU[i];
00911 value = elementU[i];
00912 int mark = markRow[iRow];
00913
00914 if ( mark == largeInteger+1 ) {
00915
00916 indexRowU[put] = iRow;
00917 elementU[put] = value;
00918 if ( checkLargest ) {
00919 double absValue = fabs ( value );
00920
00921 if ( absValue > largest ) {
00922 largest = absValue;
00923 positionLargest = put;
00924 }
00925 }
00926 put++;
00927 } else if ( mark != largeInteger ) {
00928
00929 work[mark] = value;
00930 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00931 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00932
00933 temp2[word] = temp2[word] | ( 1 << bit );
00934 added--;
00935 } else {
00936 thisPivotValue = value;
00937 }
00938 }
00939
00940 elementU[put] = elementU[startColumn];
00941 indexRowU[put] = indexRowU[startColumn];
00942 if ( positionLargest == startColumn ) {
00943 positionLargest = put;
00944 }
00945 put++;
00946 elementU[startColumn] = thisPivotValue;
00947 indexRowU[startColumn] = pivotRow;
00948
00949 startColumn++;
00950 numberInColumn[iColumn] = put - startColumn;
00951 int * numberInColumnPlus = numberInColumnPlus_.array();
00952 numberInColumnPlus[iColumn]++;
00953 startColumnU[iColumn]++;
00954
00955 int next = nextColumn[iColumn];
00956 CoinBigIndex space;
00957
00958 space = startColumnU[next] - put - numberInColumnPlus[next];
00959
00960 if ( numberInPivotColumn > space ) {
00961
00962 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
00963 return false;
00964 }
00965
00966 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
00967 startColumn = startColumnU[iColumn];
00968 put = startColumn + numberInColumn[iColumn];
00969 }
00970 double tolerance = zeroTolerance_;
00971
00972 int *nextCount = nextCount_.array();
00973 for ( j = 0; j < numberInPivotColumn; j++ ) {
00974 value = work[j] - thisPivotValue * multipliersL[j];
00975 double absValue = fabs ( value );
00976
00977 if ( absValue > tolerance ) {
00978 work[j] = 0.0;
00979 elementU[put] = value;
00980 indexRowU[put] = indexL[j];
00981 if ( absValue > largest ) {
00982 largest = absValue;
00983 positionLargest = put;
00984 }
00985 put++;
00986 } else {
00987 work[j] = 0.0;
00988 added--;
00989 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
00990 int bit = j & COINFACTORIZATION_MASK_PER_INT;
00991
00992 if ( temp2[word] & ( 1 << bit ) ) {
00993
00994 iRow = indexL[j];
00995 CoinBigIndex start = startRowU[iRow];
00996 CoinBigIndex end = start + numberInRow[iRow];
00997 CoinBigIndex where = start;
00998
00999 while ( indexColumnU[where] != iColumn ) {
01000 where++;
01001 }
01002 #if DEBUG_COIN
01003 if ( where >= end ) {
01004 abort ( );
01005 }
01006 #endif
01007 indexColumnU[where] = indexColumnU[end - 1];
01008 numberInRow[iRow]--;
01009 } else {
01010
01011 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01012 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01013
01014 temp2[word] = temp2[word] | ( 1 << bit );
01015 }
01016 }
01017 }
01018 numberInColumn[iColumn] = put - startColumn;
01019
01020 if ( positionLargest >= 0 ) {
01021 value = elementU[positionLargest];
01022 iRow = indexRowU[positionLargest];
01023 elementU[positionLargest] = elementU[startColumn];
01024 indexRowU[positionLargest] = indexRowU[startColumn];
01025 elementU[startColumn] = value;
01026 indexRowU[startColumn] = iRow;
01027 }
01028
01029 if ( nextCount[iColumn + numberRows_] != -2 ) {
01030
01031 deleteLink ( iColumn + numberRows_ );
01032 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01033 }
01034 temp2 += increment2;
01035 }
01036
01037 unsigned int *putBase = workArea2;
01038 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01039 int i = 0;
01040
01041
01042 while ( bigLoops ) {
01043 bigLoops--;
01044 int bit;
01045 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01046 unsigned int *putThis = putBase;
01047 int iRow = indexL[i];
01048
01049
01050 int number = 0;
01051 int jColumn;
01052
01053 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01054 unsigned int test = *putThis;
01055
01056 putThis += increment2;
01057 test = 1 - ( ( test >> bit ) & 1 );
01058 number += test;
01059 }
01060 int next = nextRow[iRow];
01061 CoinBigIndex space;
01062
01063 space = startRowU[next] - startRowU[iRow];
01064 number += numberInRow[iRow];
01065 if ( space < number ) {
01066 if ( !getRowSpace ( iRow, number ) ) {
01067 return false;
01068 }
01069 }
01070
01071 putThis = putBase;
01072 next = nextRow[iRow];
01073 number = numberInRow[iRow];
01074 CoinBigIndex end = startRowU[iRow] + number;
01075 int saveIndex = indexColumnU[startRowU[next]];
01076
01077
01078 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01079 unsigned int test = *putThis;
01080
01081 putThis += increment2;
01082 test = 1 - ( ( test >> bit ) & 1 );
01083 indexColumnU[end] = saveColumn[jColumn];
01084 end += test;
01085 }
01086
01087 indexColumnU[startRowU[next]] = saveIndex;
01088 markRow[iRow] = static_cast<T>(largeInteger+1);
01089 number = end - startRowU[iRow];
01090 numberInRow[iRow] = number;
01091 deleteLink ( iRow );
01092 addLink ( iRow, number );
01093 }
01094 putBase++;
01095 }
01096 int bit;
01097
01098 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01099 unsigned int *putThis = putBase;
01100 int iRow = indexL[i];
01101
01102
01103 int number = 0;
01104 int jColumn;
01105
01106 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01107 unsigned int test = *putThis;
01108
01109 putThis += increment2;
01110 test = 1 - ( ( test >> bit ) & 1 );
01111 number += test;
01112 }
01113 int next = nextRow[iRow];
01114 CoinBigIndex space;
01115
01116 space = startRowU[next] - startRowU[iRow];
01117 number += numberInRow[iRow];
01118 if ( space < number ) {
01119 if ( !getRowSpace ( iRow, number ) ) {
01120 return false;
01121 }
01122 }
01123
01124 putThis = putBase;
01125 next = nextRow[iRow];
01126 number = numberInRow[iRow];
01127 CoinBigIndex end = startRowU[iRow] + number;
01128 int saveIndex;
01129
01130 saveIndex = indexColumnU[startRowU[next]];
01131
01132
01133 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01134 unsigned int test = *putThis;
01135
01136 putThis += increment2;
01137 test = 1 - ( ( test >> bit ) & 1 );
01138
01139 indexColumnU[end] = saveColumn[jColumn];
01140 end += test;
01141 }
01142 indexColumnU[startRowU[next]] = saveIndex;
01143 markRow[iRow] = static_cast<T>(largeInteger+1);
01144 number = end - startRowU[iRow];
01145 numberInRow[iRow] = number;
01146 deleteLink ( iRow );
01147 addLink ( iRow, number );
01148 }
01149 markRow[pivotRow] = static_cast<T>(largeInteger+1);
01150
01151 deleteLink ( pivotRow );
01152 deleteLink ( pivotColumn + numberRows_ );
01153 totalElements_ += added;
01154 return true;
01155 }
01156
01157
01159
01160 protected:
01161
01164
01165 double pivotTolerance_;
01167 double zeroTolerance_;
01168 #ifndef COIN_FAST_CODE
01170 double slackValue_;
01171 #else
01172 #ifndef slackValue_
01173 #define slackValue_ -1.0
01174 #endif
01175 #endif
01177 double areaFactor_;
01179 double relaxCheck_;
01181 int numberRows_;
01183 int numberRowsExtra_;
01185 int maximumRowsExtra_;
01187 int numberColumns_;
01189 int numberColumnsExtra_;
01191 int maximumColumnsExtra_;
01193 int numberGoodU_;
01195 int numberGoodL_;
01197 int maximumPivots_;
01199 int numberPivots_;
01202 CoinBigIndex totalElements_;
01204 CoinBigIndex factorElements_;
01206 CoinIntArrayWithLength pivotColumn_;
01208 CoinIntArrayWithLength permute_;
01210 CoinIntArrayWithLength permuteBack_;
01212 CoinIntArrayWithLength pivotColumnBack_;
01214 int status_;
01215
01220
01221
01223 int numberTrials_;
01225 CoinBigIndexArrayWithLength startRowU_;
01226
01228 CoinIntArrayWithLength numberInRow_;
01229
01231 CoinIntArrayWithLength numberInColumn_;
01232
01234 CoinIntArrayWithLength numberInColumnPlus_;
01235
01238 CoinIntArrayWithLength firstCount_;
01239
01241 CoinIntArrayWithLength nextCount_;
01242
01244 CoinIntArrayWithLength lastCount_;
01245
01247 CoinIntArrayWithLength nextColumn_;
01248
01250 CoinIntArrayWithLength lastColumn_;
01251
01253 CoinIntArrayWithLength nextRow_;
01254
01256 CoinIntArrayWithLength lastRow_;
01257
01259 CoinIntArrayWithLength saveColumn_;
01260
01262 CoinIntArrayWithLength markRow_;
01263
01265 int messageLevel_;
01266
01268 int biggerDimension_;
01269
01271 CoinIntArrayWithLength indexColumnU_;
01272
01274 CoinIntArrayWithLength pivotRowL_;
01275
01277 CoinDoubleArrayWithLength pivotRegion_;
01278
01280 int numberSlacks_;
01281
01283 int numberU_;
01284
01286 CoinBigIndex maximumU_;
01287
01289
01290
01292 CoinBigIndex lengthU_;
01293
01295 CoinBigIndex lengthAreaU_;
01296
01298 CoinDoubleArrayWithLength elementU_;
01299
01301 CoinIntArrayWithLength indexRowU_;
01302
01304 CoinBigIndexArrayWithLength startColumnU_;
01305
01307 CoinBigIndexArrayWithLength convertRowToColumnU_;
01308
01310 CoinBigIndex numberL_;
01311
01313 CoinBigIndex baseL_;
01314
01316 CoinBigIndex lengthL_;
01317
01319 CoinBigIndex lengthAreaL_;
01320
01322 CoinDoubleArrayWithLength elementL_;
01323
01325 CoinIntArrayWithLength indexRowL_;
01326
01328 CoinBigIndexArrayWithLength startColumnL_;
01329
01331 bool doForrestTomlin_;
01332
01334 int numberR_;
01335
01337 CoinBigIndex lengthR_;
01338
01340 CoinBigIndex lengthAreaR_;
01341
01343 double *elementR_;
01344
01346 int *indexRowR_;
01347
01349 CoinBigIndexArrayWithLength startColumnR_;
01350
01352 double * denseArea_;
01353
01355 int * densePermute_;
01356
01358 int numberDense_;
01359
01361 int denseThreshold_;
01362
01364 CoinDoubleArrayWithLength workArea_;
01365
01367 CoinUnsignedIntArrayWithLength workArea2_;
01368
01370 CoinBigIndex numberCompressions_;
01371
01373 mutable double ftranCountInput_;
01374 mutable double ftranCountAfterL_;
01375 mutable double ftranCountAfterR_;
01376 mutable double ftranCountAfterU_;
01377 mutable double btranCountInput_;
01378 mutable double btranCountAfterU_;
01379 mutable double btranCountAfterR_;
01380 mutable double btranCountAfterL_;
01381
01383 mutable int numberFtranCounts_;
01384 mutable int numberBtranCounts_;
01385
01387 double ftranAverageAfterL_;
01388 double ftranAverageAfterR_;
01389 double ftranAverageAfterU_;
01390 double btranAverageAfterU_;
01391 double btranAverageAfterR_;
01392 double btranAverageAfterL_;
01393
01395 mutable bool collectStatistics_;
01396
01398 int sparseThreshold_;
01399
01401 int sparseThreshold2_;
01402
01404 CoinBigIndexArrayWithLength startRowL_;
01405
01407 CoinIntArrayWithLength indexColumnL_;
01408
01410 CoinDoubleArrayWithLength elementByRowL_;
01411
01413 mutable CoinIntArrayWithLength sparse_;
01417 int biasLU_;
01423 int persistenceFlag_;
01425 };
01426
01427 #ifdef COIN_HAS_LAPACK
01428 #define DENSE_CODE 1
01429
01430 #ifndef ipfint
01431
01432 typedef int ipfint;
01433 typedef const int cipfint;
01434 #endif
01435 #endif
01436 #endif
01437
01438 #ifdef UGLY_COIN_FACTOR_CODING
01439 #define FAC_UNSET (FAC_SET+1)
01440 {
01441 goodPivot=false;
01442
01443 CoinBigIndex startColumnThis = startColumn[iPivotColumn];
01444 CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
01445 int put = 0;
01446 CoinBigIndex startRowThis = startRow[iPivotRow];
01447 CoinBigIndex endRow = startRowThis + numberDoRow + 1;
01448 if ( pivotColumnPosition < 0 ) {
01449 for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01450 int iColumn = indexColumn[pivotColumnPosition];
01451 if ( iColumn != iPivotColumn ) {
01452 saveColumn[put++] = iColumn;
01453 } else {
01454 break;
01455 }
01456 }
01457 } else {
01458 for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
01459 saveColumn[put++] = indexColumn[i];
01460 }
01461 }
01462 assert (pivotColumnPosition<endRow);
01463 assert (indexColumn[pivotColumnPosition]==iPivotColumn);
01464 pivotColumnPosition++;
01465 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01466 saveColumn[put++] = indexColumn[pivotColumnPosition];
01467 }
01468
01469 int next = nextRow[iPivotRow];
01470 int last = lastRow[iPivotRow];
01471
01472 nextRow[last] = next;
01473 lastRow[next] = last;
01474 nextRow[iPivotRow] = numberGoodU_;
01475 lastRow[iPivotRow] = -2;
01476 numberInRow[iPivotRow] = 0;
01477
01478 CoinBigIndex l = lengthL_;
01479
01480 {
01481 if ( l + numberDoColumn > lengthAreaL_ ) {
01482
01483 printf("more memory needed in middle of invert\n");
01484 goto BAD_PIVOT;
01485 }
01486
01487 CoinBigIndex lSave = l;
01488
01489 pivotRowL_.array()[numberGoodL_] = iPivotRow;
01490 CoinBigIndex * startColumnL = startColumnL_.array();
01491 startColumnL[numberGoodL_] = l;
01492 numberGoodL_++;
01493 startColumnL[numberGoodL_] = l + numberDoColumn;
01494 lengthL_ += numberDoColumn;
01495 if ( pivotRowPosition < 0 ) {
01496 for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01497 int iRow = indexRow[pivotRowPosition];
01498 if ( iRow != iPivotRow ) {
01499 indexRowL[l] = iRow;
01500 elementL[l] = element[pivotRowPosition];
01501 markRow[iRow] = l - lSave;
01502 l++;
01503
01504 CoinBigIndex start = startRow[iRow];
01505 CoinBigIndex end = start + numberInRow[iRow];
01506 CoinBigIndex where = start;
01507
01508 while ( indexColumn[where] != iPivotColumn ) {
01509 where++;
01510 }
01511 #if DEBUG_COIN
01512 if ( where >= end ) {
01513 abort ( );
01514 }
01515 #endif
01516 indexColumn[where] = indexColumn[end - 1];
01517 numberInRow[iRow]--;
01518 } else {
01519 break;
01520 }
01521 }
01522 } else {
01523 CoinBigIndex i;
01524
01525 for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
01526 int iRow = indexRow[i];
01527
01528 markRow[iRow] = l - lSave;
01529 indexRowL[l] = iRow;
01530 elementL[l] = element[i];
01531 l++;
01532
01533 CoinBigIndex start = startRow[iRow];
01534 CoinBigIndex end = start + numberInRow[iRow];
01535 CoinBigIndex where = start;
01536
01537 while ( indexColumn[where] != iPivotColumn ) {
01538 where++;
01539 }
01540 #if DEBUG_COIN
01541 if ( where >= end ) {
01542 abort ( );
01543 }
01544 #endif
01545 indexColumn[where] = indexColumn[end - 1];
01546 numberInRow[iRow]--;
01547 assert (numberInRow[iRow]>=0);
01548 }
01549 }
01550 assert (pivotRowPosition<endColumn);
01551 assert (indexRow[pivotRowPosition]==iPivotRow);
01552 double pivotElement = element[pivotRowPosition];
01553 double pivotMultiplier = 1.0 / pivotElement;
01554
01555 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
01556 pivotRowPosition++;
01557 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01558 int iRow = indexRow[pivotRowPosition];
01559
01560 markRow[iRow] = l - lSave;
01561 indexRowL[l] = iRow;
01562 elementL[l] = element[pivotRowPosition];
01563 l++;
01564
01565 CoinBigIndex start = startRow[iRow];
01566 CoinBigIndex end = start + numberInRow[iRow];
01567 CoinBigIndex where = start;
01568
01569 while ( indexColumn[where] != iPivotColumn ) {
01570 where++;
01571 }
01572 #if DEBUG_COIN
01573 if ( where >= end ) {
01574 abort ( );
01575 }
01576 #endif
01577 indexColumn[where] = indexColumn[end - 1];
01578 numberInRow[iRow]--;
01579 assert (numberInRow[iRow]>=0);
01580 }
01581 markRow[iPivotRow] = FAC_SET;
01582
01583 numberInColumn[iPivotColumn] = 0;
01584
01585 int *indexL = &indexRowL[lSave];
01586 double *multipliersL = &elementL[lSave];
01587
01588
01589 int j;
01590
01591 for ( j = 0; j < numberDoColumn; j++ ) {
01592 multipliersL[j] *= pivotMultiplier;
01593 }
01594
01595 CoinBigIndex iErase;
01596 for ( iErase = 0; iErase < increment2 * numberDoRow;
01597 iErase++ ) {
01598 workArea2[iErase] = 0;
01599 }
01600 CoinBigIndex added = numberDoRow * numberDoColumn;
01601 unsigned int *temp2 = workArea2;
01602 int * nextColumn = nextColumn_.array();
01603
01604
01605 int jColumn;
01606 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01607 int iColumn = saveColumn[jColumn];
01608 CoinBigIndex startColumnThis = startColumn[iColumn];
01609 CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
01610 int iRow = indexRow[startColumnThis];
01611 double value = element[startColumnThis];
01612 double largest;
01613 CoinBigIndex put = startColumnThis;
01614 CoinBigIndex positionLargest = -1;
01615 double thisPivotValue = 0.0;
01616
01617
01618 bool checkLargest;
01619 int mark = markRow[iRow];
01620
01621 if ( mark == FAC_UNSET ) {
01622 largest = fabs ( value );
01623 positionLargest = put;
01624 put++;
01625 checkLargest = false;
01626 } else {
01627
01628 largest = 0.0;
01629 checkLargest = true;
01630 if ( mark != FAC_SET ) {
01631
01632 workArea[mark] = value;
01633 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01634 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01635
01636 temp2[word] = temp2[word] | ( 1 << bit );
01637 added--;
01638 } else {
01639 thisPivotValue = value;
01640 }
01641 }
01642 CoinBigIndex i;
01643 for ( i = startColumnThis + 1; i < endColumn; i++ ) {
01644 iRow = indexRow[i];
01645 value = element[i];
01646 int mark = markRow[iRow];
01647
01648 if ( mark == FAC_UNSET ) {
01649
01650 indexRow[put] = iRow;
01651 element[put] = value;
01652 if ( checkLargest ) {
01653 double absValue = fabs ( value );
01654
01655 if ( absValue > largest ) {
01656 largest = absValue;
01657 positionLargest = put;
01658 }
01659 }
01660 put++;
01661 } else if ( mark != FAC_SET ) {
01662
01663 workArea[mark] = value;
01664 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01665 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01666
01667 temp2[word] = temp2[word] | ( 1 << bit );
01668 added--;
01669 } else {
01670 thisPivotValue = value;
01671 }
01672 }
01673
01674 element[put] = element[startColumnThis];
01675 indexRow[put] = indexRow[startColumnThis];
01676 if ( positionLargest == startColumnThis ) {
01677 positionLargest = put;
01678 }
01679 put++;
01680 element[startColumnThis] = thisPivotValue;
01681 indexRow[startColumnThis] = iPivotRow;
01682
01683 startColumnThis++;
01684 numberInColumn[iColumn] = put - startColumnThis;
01685 int * numberInColumnPlus = numberInColumnPlus_.array();
01686 numberInColumnPlus[iColumn]++;
01687 startColumn[iColumn]++;
01688
01689 int next = nextColumn[iColumn];
01690 CoinBigIndex space;
01691
01692 space = startColumn[next] - put - numberInColumnPlus[next];
01693
01694 if ( numberDoColumn > space ) {
01695
01696 if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
01697 goto BAD_PIVOT;
01698 }
01699
01700 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
01701 startColumnThis = startColumn[iColumn];
01702 put = startColumnThis + numberInColumn[iColumn];
01703 }
01704 double tolerance = zeroTolerance_;
01705
01706 int *nextCount = nextCount_.array();
01707 for ( j = 0; j < numberDoColumn; j++ ) {
01708 value = workArea[j] - thisPivotValue * multipliersL[j];
01709 double absValue = fabs ( value );
01710
01711 if ( absValue > tolerance ) {
01712 workArea[j] = 0.0;
01713 element[put] = value;
01714 indexRow[put] = indexL[j];
01715 if ( absValue > largest ) {
01716 largest = absValue;
01717 positionLargest = put;
01718 }
01719 put++;
01720 } else {
01721 workArea[j] = 0.0;
01722 added--;
01723 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01724 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01725
01726 if ( temp2[word] & ( 1 << bit ) ) {
01727
01728 iRow = indexL[j];
01729 CoinBigIndex start = startRow[iRow];
01730 CoinBigIndex end = start + numberInRow[iRow];
01731 CoinBigIndex where = start;
01732
01733 while ( indexColumn[where] != iColumn ) {
01734 where++;
01735 }
01736 #if DEBUG_COIN
01737 if ( where >= end ) {
01738 abort ( );
01739 }
01740 #endif
01741 indexColumn[where] = indexColumn[end - 1];
01742 numberInRow[iRow]--;
01743 } else {
01744
01745 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01746 int bit = j & COINFACTORIZATION_MASK_PER_INT;
01747
01748 temp2[word] = temp2[word] | ( 1 << bit );
01749 }
01750 }
01751 }
01752 numberInColumn[iColumn] = put - startColumnThis;
01753
01754 if ( positionLargest >= 0 ) {
01755 value = element[positionLargest];
01756 iRow = indexRow[positionLargest];
01757 element[positionLargest] = element[startColumnThis];
01758 indexRow[positionLargest] = indexRow[startColumnThis];
01759 element[startColumnThis] = value;
01760 indexRow[startColumnThis] = iRow;
01761 }
01762
01763 if ( nextCount[iColumn + numberRows_] != -2 ) {
01764
01765 deleteLink ( iColumn + numberRows_ );
01766 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01767 }
01768 temp2 += increment2;
01769 }
01770
01771 unsigned int *putBase = workArea2;
01772 int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01773 int i = 0;
01774
01775
01776 while ( bigLoops ) {
01777 bigLoops--;
01778 int bit;
01779 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01780 unsigned int *putThis = putBase;
01781 int iRow = indexL[i];
01782
01783
01784 int number = 0;
01785 int jColumn;
01786
01787 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01788 unsigned int test = *putThis;
01789
01790 putThis += increment2;
01791 test = 1 - ( ( test >> bit ) & 1 );
01792 number += test;
01793 }
01794 int next = nextRow[iRow];
01795 CoinBigIndex space;
01796
01797 space = startRow[next] - startRow[iRow];
01798 number += numberInRow[iRow];
01799 if ( space < number ) {
01800 if ( !getRowSpace ( iRow, number ) ) {
01801 goto BAD_PIVOT;
01802 }
01803 }
01804
01805 putThis = putBase;
01806 next = nextRow[iRow];
01807 number = numberInRow[iRow];
01808 CoinBigIndex end = startRow[iRow] + number;
01809 int saveIndex = indexColumn[startRow[next]];
01810
01811
01812 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01813 unsigned int test = *putThis;
01814
01815 putThis += increment2;
01816 test = 1 - ( ( test >> bit ) & 1 );
01817 indexColumn[end] = saveColumn[jColumn];
01818 end += test;
01819 }
01820
01821 indexColumn[startRow[next]] = saveIndex;
01822 markRow[iRow] = FAC_UNSET;
01823 number = end - startRow[iRow];
01824 numberInRow[iRow] = number;
01825 deleteLink ( iRow );
01826 addLink ( iRow, number );
01827 }
01828 putBase++;
01829 }
01830 int bit;
01831
01832 for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
01833 unsigned int *putThis = putBase;
01834 int iRow = indexL[i];
01835
01836
01837 int number = 0;
01838 int jColumn;
01839
01840 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01841 unsigned int test = *putThis;
01842
01843 putThis += increment2;
01844 test = 1 - ( ( test >> bit ) & 1 );
01845 number += test;
01846 }
01847 int next = nextRow[iRow];
01848 CoinBigIndex space;
01849
01850 space = startRow[next] - startRow[iRow];
01851 number += numberInRow[iRow];
01852 if ( space < number ) {
01853 if ( !getRowSpace ( iRow, number ) ) {
01854 goto BAD_PIVOT;
01855 }
01856 }
01857
01858 putThis = putBase;
01859 next = nextRow[iRow];
01860 number = numberInRow[iRow];
01861 CoinBigIndex end = startRow[iRow] + number;
01862 int saveIndex;
01863
01864 saveIndex = indexColumn[startRow[next]];
01865
01866
01867 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01868 unsigned int test = *putThis;
01869
01870 putThis += increment2;
01871 test = 1 - ( ( test >> bit ) & 1 );
01872
01873 indexColumn[end] = saveColumn[jColumn];
01874 end += test;
01875 }
01876 indexColumn[startRow[next]] = saveIndex;
01877 markRow[iRow] = FAC_UNSET;
01878 number = end - startRow[iRow];
01879 numberInRow[iRow] = number;
01880 deleteLink ( iRow );
01881 addLink ( iRow, number );
01882 }
01883 markRow[iPivotRow] = FAC_UNSET;
01884
01885 deleteLink ( iPivotRow );
01886 deleteLink ( iPivotColumn + numberRows_ );
01887 totalElements_ += added;
01888 goodPivot= true;
01889
01890 }
01891 BAD_PIVOT:
01892
01893 ;
01894 }
01895 #undef FAC_UNSET
01896 #endif