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