CoinUtils trunk
|
00001 /* $Id$ */ 00002 // Copyright (C) 2002, International Business Machines 00003 // Corporation and others. All Rights Reserved. 00004 // This code is licensed under the terms of the Eclipse Public License (EPL). 00005 00006 /* 00007 Authors 00008 00009 John Forrest 00010 00011 */ 00012 #ifndef CoinFactorization_H 00013 #define CoinFactorization_H 00014 //#define COIN_ONE_ETA_COPY 100 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 //return firstCount_.array(); 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 00272 00273 inline double slackValue ( ) const { 00274 return slackValue_ ; 00275 } 00276 void slackValue ( double value ); 00277 #endif 00278 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 //first with that count 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 /********************************* START LARGE TEMPLATE ********/ 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 //store pivot columns (so can easily compress) 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 //take out this bit of indexColumnU 00756 int next = nextRow[pivotRow]; 00757 int last = lastRow[pivotRow]; 00758 00759 nextRow[last] = next; 00760 lastRow[next] = last; 00761 nextRow[pivotRow] = numberGoodU_; //use for permute 00762 lastRow[pivotRow] = -2; 00763 numberInRow[pivotRow] = 0; 00764 //store column in L, compress in U and take column out 00765 CoinBigIndex l = lengthL_; 00766 00767 if ( l + numberInPivotColumn > lengthAreaL_ ) { 00768 //need more memory 00769 if ((messageLevel_&4)!=0) 00770 printf("more memory needed in middle of invert\n"); 00771 return false; 00772 } 00773 //l+=currentAreaL_->elementByColumn-elementL; 00774 CoinBigIndex lSave = l; 00775 00776 CoinBigIndex * startColumnL = startColumnL_.array(); 00777 startColumnL[numberGoodL_] = l; //for luck and first time 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 //take out of row list 00790 CoinBigIndex start = startRowU[iRow]; 00791 CoinBigIndex end = start + numberInRow[iRow]; 00792 CoinBigIndex where = start; 00793 00794 while ( indexColumnU[where] != pivotColumn ) { 00795 where++; 00796 } /* endwhile */ 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 //take out of row list 00819 CoinBigIndex start = startRowU[iRow]; 00820 CoinBigIndex end = start + numberInRow[iRow]; 00821 CoinBigIndex where = start; 00822 00823 while ( indexColumnU[where] != pivotColumn ) { 00824 where++; 00825 } /* endwhile */ 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 //take out of row list 00851 CoinBigIndex start = startRowU[iRow]; 00852 CoinBigIndex end = start + numberInRow[iRow]; 00853 CoinBigIndex where = start; 00854 00855 while ( indexColumnU[where] != pivotColumn ) { 00856 where++; 00857 } /* endwhile */ 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 //compress pivot column (move pivot to front including saved) 00869 numberInColumn[pivotColumn] = 0; 00870 //use end of L for temporary space 00871 int *indexL = &indexRowL[lSave]; 00872 CoinFactorizationDouble *multipliersL = &elementL[lSave]; 00873 00874 //adjust 00875 int j; 00876 00877 for ( j = 0; j < numberInPivotColumn; j++ ) { 00878 multipliersL[j] *= pivotMultiplier; 00879 } 00880 //zero out fill 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 //pack down and move to work 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 //compress column and find largest not updated 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 //need to find largest 00914 largest = 0.0; 00915 checkLargest = true; 00916 if ( mark != largeInteger ) { 00917 //will be updated 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 ); //say already in counts 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 //keep 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 //will be updated 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 ); //say already in counts 00954 added--; 00955 } else { 00956 thisPivotValue = value; 00957 } 00958 } 00959 //slot in pivot 00960 elementU[put] = elementU[startColumn]; 00961 indexRowU[put] = indexRowU[startColumn]; 00962 if ( positionLargest == startColumn ) { 00963 positionLargest = put; //follow if was largest 00964 } 00965 put++; 00966 elementU[startColumn] = thisPivotValue; 00967 indexRowU[startColumn] = pivotRow; 00968 //clean up counts 00969 startColumn++; 00970 numberInColumn[iColumn] = put - startColumn; 00971 int * numberInColumnPlus = numberInColumnPlus_.array(); 00972 numberInColumnPlus[iColumn]++; 00973 startColumnU[iColumn]++; 00974 //how much space have we got 00975 int next = nextColumn[iColumn]; 00976 CoinBigIndex space; 00977 00978 space = startColumnU[next] - put - numberInColumnPlus[next]; 00979 //assume no zero elements 00980 if ( numberInPivotColumn > space ) { 00981 //getColumnSpace also moves fixed part 00982 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) { 00983 return false; 00984 } 00985 //redo starts 00986 positionLargest = positionLargest + startColumnU[iColumn] - startColumn; 00987 startColumn = startColumnU[iColumn]; 00988 put = startColumn + numberInColumn[iColumn]; 00989 } 00990 double tolerance = zeroTolerance_; 00991 00992 int *nextCount = nextCount_.array(); 00993 for ( j = 0; j < numberInPivotColumn; j++ ) { 00994 value = work[j] - thisPivotValue * multipliersL[j]; 00995 double absValue = fabs ( value ); 00996 00997 if ( absValue > tolerance ) { 00998 work[j] = 0.0; 00999 assert (put<lengthAreaU_); 01000 elementU[put] = value; 01001 indexRowU[put] = indexL[j]; 01002 if ( absValue > largest ) { 01003 largest = absValue; 01004 positionLargest = put; 01005 } 01006 put++; 01007 } else { 01008 work[j] = 0.0; 01009 added--; 01010 int word = j >> COINFACTORIZATION_SHIFT_PER_INT; 01011 int bit = j & COINFACTORIZATION_MASK_PER_INT; 01012 01013 if ( temp2[word] & ( 1 << bit ) ) { 01014 //take out of row list 01015 iRow = indexL[j]; 01016 CoinBigIndex start = startRowU[iRow]; 01017 CoinBigIndex end = start + numberInRow[iRow]; 01018 CoinBigIndex where = start; 01019 01020 while ( indexColumnU[where] != iColumn ) { 01021 where++; 01022 } /* endwhile */ 01023 #if DEBUG_COIN 01024 if ( where >= end ) { 01025 abort ( ); 01026 } 01027 #endif 01028 indexColumnU[where] = indexColumnU[end - 1]; 01029 numberInRow[iRow]--; 01030 } else { 01031 //make sure won't be added 01032 int word = j >> COINFACTORIZATION_SHIFT_PER_INT; 01033 int bit = j & COINFACTORIZATION_MASK_PER_INT; 01034 01035 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts 01036 } 01037 } 01038 } 01039 numberInColumn[iColumn] = put - startColumn; 01040 //move largest 01041 if ( positionLargest >= 0 ) { 01042 value = elementU[positionLargest]; 01043 iRow = indexRowU[positionLargest]; 01044 elementU[positionLargest] = elementU[startColumn]; 01045 indexRowU[positionLargest] = indexRowU[startColumn]; 01046 elementU[startColumn] = value; 01047 indexRowU[startColumn] = iRow; 01048 } 01049 //linked list for column 01050 if ( nextCount[iColumn + numberRows_] != -2 ) { 01051 //modify linked list 01052 deleteLink ( iColumn + numberRows_ ); 01053 addLink ( iColumn + numberRows_, numberInColumn[iColumn] ); 01054 } 01055 temp2 += increment2; 01056 } 01057 //get space for row list 01058 unsigned int *putBase = workArea2; 01059 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT; 01060 int i = 0; 01061 01062 // do linked lists and update counts 01063 while ( bigLoops ) { 01064 bigLoops--; 01065 int bit; 01066 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) { 01067 unsigned int *putThis = putBase; 01068 int iRow = indexL[i]; 01069 01070 //get space 01071 int number = 0; 01072 int jColumn; 01073 01074 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) { 01075 unsigned int test = *putThis; 01076 01077 putThis += increment2; 01078 test = 1 - ( ( test >> bit ) & 1 ); 01079 number += test; 01080 } 01081 int next = nextRow[iRow]; 01082 CoinBigIndex space; 01083 01084 space = startRowU[next] - startRowU[iRow]; 01085 number += numberInRow[iRow]; 01086 if ( space < number ) { 01087 if ( !getRowSpace ( iRow, number ) ) { 01088 return false; 01089 } 01090 } 01091 // now do 01092 putThis = putBase; 01093 next = nextRow[iRow]; 01094 number = numberInRow[iRow]; 01095 CoinBigIndex end = startRowU[iRow] + number; 01096 int saveIndex = indexColumnU[startRowU[next]]; 01097 01098 //add in 01099 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) { 01100 unsigned int test = *putThis; 01101 01102 putThis += increment2; 01103 test = 1 - ( ( test >> bit ) & 1 ); 01104 indexColumnU[end] = saveColumn[jColumn]; 01105 end += test; 01106 } 01107 //put back next one in case zapped 01108 indexColumnU[startRowU[next]] = saveIndex; 01109 markRow[iRow] = static_cast<T>(largeInteger+1); 01110 number = end - startRowU[iRow]; 01111 numberInRow[iRow] = number; 01112 deleteLink ( iRow ); 01113 addLink ( iRow, number ); 01114 } 01115 putBase++; 01116 } /* endwhile */ 01117 int bit; 01118 01119 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) { 01120 unsigned int *putThis = putBase; 01121 int iRow = indexL[i]; 01122 01123 //get space 01124 int number = 0; 01125 int jColumn; 01126 01127 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) { 01128 unsigned int test = *putThis; 01129 01130 putThis += increment2; 01131 test = 1 - ( ( test >> bit ) & 1 ); 01132 number += test; 01133 } 01134 int next = nextRow[iRow]; 01135 CoinBigIndex space; 01136 01137 space = startRowU[next] - startRowU[iRow]; 01138 number += numberInRow[iRow]; 01139 if ( space < number ) { 01140 if ( !getRowSpace ( iRow, number ) ) { 01141 return false; 01142 } 01143 } 01144 // now do 01145 putThis = putBase; 01146 next = nextRow[iRow]; 01147 number = numberInRow[iRow]; 01148 CoinBigIndex end = startRowU[iRow] + number; 01149 int saveIndex; 01150 01151 saveIndex = indexColumnU[startRowU[next]]; 01152 01153 //add in 01154 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) { 01155 unsigned int test = *putThis; 01156 01157 putThis += increment2; 01158 test = 1 - ( ( test >> bit ) & 1 ); 01159 01160 indexColumnU[end] = saveColumn[jColumn]; 01161 end += test; 01162 } 01163 indexColumnU[startRowU[next]] = saveIndex; 01164 markRow[iRow] = static_cast<T>(largeInteger+1); 01165 number = end - startRowU[iRow]; 01166 numberInRow[iRow] = number; 01167 deleteLink ( iRow ); 01168 addLink ( iRow, number ); 01169 } 01170 markRow[pivotRow] = static_cast<T>(largeInteger+1); 01171 //modify linked list for pivots 01172 deleteLink ( pivotRow ); 01173 deleteLink ( pivotColumn + numberRows_ ); 01174 totalElements_ += added; 01175 return true; 01176 } 01177 01178 /********************************* END LARGE TEMPLATE ********/ 01180 01181 protected: 01182 01185 01186 double pivotTolerance_; 01188 double zeroTolerance_; 01189 #ifndef COIN_FAST_CODE 01190 01191 double slackValue_; 01192 #else 01193 #ifndef slackValue_ 01194 #define slackValue_ -1.0 01195 #endif 01196 #endif 01197 01198 double areaFactor_; 01200 double relaxCheck_; 01202 int numberRows_; 01204 int numberRowsExtra_; 01206 int maximumRowsExtra_; 01208 int numberColumns_; 01210 int numberColumnsExtra_; 01212 int maximumColumnsExtra_; 01214 int numberGoodU_; 01216 int numberGoodL_; 01218 int maximumPivots_; 01220 int numberPivots_; 01223 CoinBigIndex totalElements_; 01225 CoinBigIndex factorElements_; 01227 CoinIntArrayWithLength pivotColumn_; 01229 CoinIntArrayWithLength permute_; 01231 CoinIntArrayWithLength permuteBack_; 01233 CoinIntArrayWithLength pivotColumnBack_; 01235 int status_; 01236 01241 //int increasingRows_; 01242 01244 int numberTrials_; 01246 CoinBigIndexArrayWithLength startRowU_; 01247 01249 CoinIntArrayWithLength numberInRow_; 01250 01252 CoinIntArrayWithLength numberInColumn_; 01253 01255 CoinIntArrayWithLength numberInColumnPlus_; 01256 01259 CoinIntArrayWithLength firstCount_; 01260 01262 CoinIntArrayWithLength nextCount_; 01263 01265 CoinIntArrayWithLength lastCount_; 01266 01268 CoinIntArrayWithLength nextColumn_; 01269 01271 CoinIntArrayWithLength lastColumn_; 01272 01274 CoinIntArrayWithLength nextRow_; 01275 01277 CoinIntArrayWithLength lastRow_; 01278 01280 CoinIntArrayWithLength saveColumn_; 01281 01283 CoinIntArrayWithLength markRow_; 01284 01286 int messageLevel_; 01287 01289 int biggerDimension_; 01290 01292 CoinIntArrayWithLength indexColumnU_; 01293 01295 CoinIntArrayWithLength pivotRowL_; 01296 01298 CoinFactorizationDoubleArrayWithLength pivotRegion_; 01299 01301 int numberSlacks_; 01302 01304 int numberU_; 01305 01307 CoinBigIndex maximumU_; 01308 01310 //int baseU_; 01311 01313 CoinBigIndex lengthU_; 01314 01316 CoinBigIndex lengthAreaU_; 01317 01319 CoinFactorizationDoubleArrayWithLength elementU_; 01320 01322 CoinIntArrayWithLength indexRowU_; 01323 01325 CoinBigIndexArrayWithLength startColumnU_; 01326 01328 CoinBigIndexArrayWithLength convertRowToColumnU_; 01329 01331 CoinBigIndex numberL_; 01332 01334 CoinBigIndex baseL_; 01335 01337 CoinBigIndex lengthL_; 01338 01340 CoinBigIndex lengthAreaL_; 01341 01343 CoinFactorizationDoubleArrayWithLength elementL_; 01344 01346 CoinIntArrayWithLength indexRowL_; 01347 01349 CoinBigIndexArrayWithLength startColumnL_; 01350 01352 bool doForrestTomlin_; 01353 01355 int numberR_; 01356 01358 CoinBigIndex lengthR_; 01359 01361 CoinBigIndex lengthAreaR_; 01362 01364 CoinFactorizationDouble *elementR_; 01365 01367 int *indexRowR_; 01368 01370 CoinBigIndexArrayWithLength startColumnR_; 01371 01373 double * denseArea_; 01374 01376 int * densePermute_; 01377 01379 int numberDense_; 01380 01382 int denseThreshold_; 01383 01385 CoinFactorizationDoubleArrayWithLength workArea_; 01386 01388 CoinUnsignedIntArrayWithLength workArea2_; 01389 01391 CoinBigIndex numberCompressions_; 01392 01394 mutable double ftranCountInput_; 01395 mutable double ftranCountAfterL_; 01396 mutable double ftranCountAfterR_; 01397 mutable double ftranCountAfterU_; 01398 mutable double btranCountInput_; 01399 mutable double btranCountAfterU_; 01400 mutable double btranCountAfterR_; 01401 mutable double btranCountAfterL_; 01402 01404 mutable int numberFtranCounts_; 01405 mutable int numberBtranCounts_; 01406 01408 double ftranAverageAfterL_; 01409 double ftranAverageAfterR_; 01410 double ftranAverageAfterU_; 01411 double btranAverageAfterU_; 01412 double btranAverageAfterR_; 01413 double btranAverageAfterL_; 01414 01416 mutable bool collectStatistics_; 01417 01419 int sparseThreshold_; 01420 01422 int sparseThreshold2_; 01423 01425 CoinBigIndexArrayWithLength startRowL_; 01426 01428 CoinIntArrayWithLength indexColumnL_; 01429 01431 CoinFactorizationDoubleArrayWithLength elementByRowL_; 01432 01434 mutable CoinIntArrayWithLength sparse_; 01438 int biasLU_; 01444 int persistenceFlag_; 01446 }; 01447 // Dense coding 01448 #ifdef COIN_HAS_LAPACK 01449 #define DENSE_CODE 1 01450 /* Type of Fortran integer translated into C */ 01451 #ifndef ipfint 01452 //typedef ipfint FORTRAN_INTEGER_TYPE ; 01453 typedef int ipfint; 01454 typedef const int cipfint; 01455 #endif 01456 #endif 01457 #endif 01458 // Extra for ugly include 01459 #ifdef UGLY_COIN_FACTOR_CODING 01460 #define FAC_UNSET (FAC_SET+1) 01461 { 01462 goodPivot=false; 01463 //store pivot columns (so can easily compress) 01464 CoinBigIndex startColumnThis = startColumn[iPivotColumn]; 01465 CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1; 01466 int put = 0; 01467 CoinBigIndex startRowThis = startRow[iPivotRow]; 01468 CoinBigIndex endRow = startRowThis + numberDoRow + 1; 01469 if ( pivotColumnPosition < 0 ) { 01470 for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) { 01471 int iColumn = indexColumn[pivotColumnPosition]; 01472 if ( iColumn != iPivotColumn ) { 01473 saveColumn[put++] = iColumn; 01474 } else { 01475 break; 01476 } 01477 } 01478 } else { 01479 for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) { 01480 saveColumn[put++] = indexColumn[i]; 01481 } 01482 } 01483 assert (pivotColumnPosition<endRow); 01484 assert (indexColumn[pivotColumnPosition]==iPivotColumn); 01485 pivotColumnPosition++; 01486 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) { 01487 saveColumn[put++] = indexColumn[pivotColumnPosition]; 01488 } 01489 //take out this bit of indexColumn 01490 int next = nextRow[iPivotRow]; 01491 int last = lastRow[iPivotRow]; 01492 01493 nextRow[last] = next; 01494 lastRow[next] = last; 01495 nextRow[iPivotRow] = numberGoodU_; //use for permute 01496 lastRow[iPivotRow] = -2; 01497 numberInRow[iPivotRow] = 0; 01498 //store column in L, compress in U and take column out 01499 CoinBigIndex l = lengthL_; 01500 // **** HORRID coding coming up but a goto seems best! 01501 { 01502 if ( l + numberDoColumn > lengthAreaL_ ) { 01503 //need more memory 01504 if ((messageLevel_&4)!=0) 01505 printf("more memory needed in middle of invert\n"); 01506 goto BAD_PIVOT; 01507 } 01508 //l+=currentAreaL_->elementByColumn-elementL; 01509 CoinBigIndex lSave = l; 01510 01511 CoinBigIndex * startColumnL = startColumnL_.array(); 01512 startColumnL[numberGoodL_] = l; //for luck and first time 01513 numberGoodL_++; 01514 startColumnL[numberGoodL_] = l + numberDoColumn; 01515 lengthL_ += numberDoColumn; 01516 if ( pivotRowPosition < 0 ) { 01517 for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) { 01518 int iRow = indexRow[pivotRowPosition]; 01519 if ( iRow != iPivotRow ) { 01520 indexRowL[l] = iRow; 01521 elementL[l] = element[pivotRowPosition]; 01522 markRow[iRow] = l - lSave; 01523 l++; 01524 //take out of row list 01525 CoinBigIndex start = startRow[iRow]; 01526 CoinBigIndex end = start + numberInRow[iRow]; 01527 CoinBigIndex where = start; 01528 01529 while ( indexColumn[where] != iPivotColumn ) { 01530 where++; 01531 } /* endwhile */ 01532 #if DEBUG_COIN 01533 if ( where >= end ) { 01534 abort ( ); 01535 } 01536 #endif 01537 indexColumn[where] = indexColumn[end - 1]; 01538 numberInRow[iRow]--; 01539 } else { 01540 break; 01541 } 01542 } 01543 } else { 01544 CoinBigIndex i; 01545 01546 for ( i = startColumnThis; i < pivotRowPosition; i++ ) { 01547 int iRow = indexRow[i]; 01548 01549 markRow[iRow] = l - lSave; 01550 indexRowL[l] = iRow; 01551 elementL[l] = element[i]; 01552 l++; 01553 //take out of row list 01554 CoinBigIndex start = startRow[iRow]; 01555 CoinBigIndex end = start + numberInRow[iRow]; 01556 CoinBigIndex where = start; 01557 01558 while ( indexColumn[where] != iPivotColumn ) { 01559 where++; 01560 } /* endwhile */ 01561 #if DEBUG_COIN 01562 if ( where >= end ) { 01563 abort ( ); 01564 } 01565 #endif 01566 indexColumn[where] = indexColumn[end - 1]; 01567 numberInRow[iRow]--; 01568 assert (numberInRow[iRow]>=0); 01569 } 01570 } 01571 assert (pivotRowPosition<endColumn); 01572 assert (indexRow[pivotRowPosition]==iPivotRow); 01573 CoinFactorizationDouble pivotElement = element[pivotRowPosition]; 01574 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement; 01575 01576 pivotRegion_.array()[numberGoodU_] = pivotMultiplier; 01577 pivotRowPosition++; 01578 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) { 01579 int iRow = indexRow[pivotRowPosition]; 01580 01581 markRow[iRow] = l - lSave; 01582 indexRowL[l] = iRow; 01583 elementL[l] = element[pivotRowPosition]; 01584 l++; 01585 //take out of row list 01586 CoinBigIndex start = startRow[iRow]; 01587 CoinBigIndex end = start + numberInRow[iRow]; 01588 CoinBigIndex where = start; 01589 01590 while ( indexColumn[where] != iPivotColumn ) { 01591 where++; 01592 } /* endwhile */ 01593 #if DEBUG_COIN 01594 if ( where >= end ) { 01595 abort ( ); 01596 } 01597 #endif 01598 indexColumn[where] = indexColumn[end - 1]; 01599 numberInRow[iRow]--; 01600 assert (numberInRow[iRow]>=0); 01601 } 01602 markRow[iPivotRow] = FAC_SET; 01603 //compress pivot column (move pivot to front including saved) 01604 numberInColumn[iPivotColumn] = 0; 01605 //use end of L for temporary space 01606 int *indexL = &indexRowL[lSave]; 01607 CoinFactorizationDouble *multipliersL = &elementL[lSave]; 01608 01609 //adjust 01610 int j; 01611 01612 for ( j = 0; j < numberDoColumn; j++ ) { 01613 multipliersL[j] *= pivotMultiplier; 01614 } 01615 //zero out fill 01616 CoinBigIndex iErase; 01617 for ( iErase = 0; iErase < increment2 * numberDoRow; 01618 iErase++ ) { 01619 workArea2[iErase] = 0; 01620 } 01621 CoinBigIndex added = numberDoRow * numberDoColumn; 01622 unsigned int *temp2 = workArea2; 01623 int * nextColumn = nextColumn_.array(); 01624 01625 //pack down and move to work 01626 int jColumn; 01627 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) { 01628 int iColumn = saveColumn[jColumn]; 01629 CoinBigIndex startColumnThis = startColumn[iColumn]; 01630 CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn]; 01631 int iRow = indexRow[startColumnThis]; 01632 CoinFactorizationDouble value = element[startColumnThis]; 01633 double largest; 01634 CoinBigIndex put = startColumnThis; 01635 CoinBigIndex positionLargest = -1; 01636 CoinFactorizationDouble thisPivotValue = 0.0; 01637 01638 //compress column and find largest not updated 01639 bool checkLargest; 01640 int mark = markRow[iRow]; 01641 01642 if ( mark == FAC_UNSET ) { 01643 largest = fabs ( value ); 01644 positionLargest = put; 01645 put++; 01646 checkLargest = false; 01647 } else { 01648 //need to find largest 01649 largest = 0.0; 01650 checkLargest = true; 01651 if ( mark != FAC_SET ) { 01652 //will be updated 01653 workArea[mark] = value; 01654 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT; 01655 int bit = mark & COINFACTORIZATION_MASK_PER_INT; 01656 01657 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts 01658 added--; 01659 } else { 01660 thisPivotValue = value; 01661 } 01662 } 01663 CoinBigIndex i; 01664 for ( i = startColumnThis + 1; i < endColumn; i++ ) { 01665 iRow = indexRow[i]; 01666 value = element[i]; 01667 int mark = markRow[iRow]; 01668 01669 if ( mark == FAC_UNSET ) { 01670 //keep 01671 indexRow[put] = iRow; 01672 element[put] = value; 01673 if ( checkLargest ) { 01674 double absValue = fabs ( value ); 01675 01676 if ( absValue > largest ) { 01677 largest = absValue; 01678 positionLargest = put; 01679 } 01680 } 01681 put++; 01682 } else if ( mark != FAC_SET ) { 01683 //will be updated 01684 workArea[mark] = value; 01685 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT; 01686 int bit = mark & COINFACTORIZATION_MASK_PER_INT; 01687 01688 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts 01689 added--; 01690 } else { 01691 thisPivotValue = value; 01692 } 01693 } 01694 //slot in pivot 01695 element[put] = element[startColumnThis]; 01696 indexRow[put] = indexRow[startColumnThis]; 01697 if ( positionLargest == startColumnThis ) { 01698 positionLargest = put; //follow if was largest 01699 } 01700 put++; 01701 element[startColumnThis] = thisPivotValue; 01702 indexRow[startColumnThis] = iPivotRow; 01703 //clean up counts 01704 startColumnThis++; 01705 numberInColumn[iColumn] = put - startColumnThis; 01706 int * numberInColumnPlus = numberInColumnPlus_.array(); 01707 numberInColumnPlus[iColumn]++; 01708 startColumn[iColumn]++; 01709 //how much space have we got 01710 int next = nextColumn[iColumn]; 01711 CoinBigIndex space; 01712 01713 space = startColumn[next] - put - numberInColumnPlus[next]; 01714 //assume no zero elements 01715 if ( numberDoColumn > space ) { 01716 //getColumnSpace also moves fixed part 01717 if ( !getColumnSpace ( iColumn, numberDoColumn ) ) { 01718 goto BAD_PIVOT; 01719 } 01720 //redo starts 01721 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis; 01722 startColumnThis = startColumn[iColumn]; 01723 put = startColumnThis + numberInColumn[iColumn]; 01724 } 01725 double tolerance = zeroTolerance_; 01726 01727 int *nextCount = nextCount_.array(); 01728 for ( j = 0; j < numberDoColumn; j++ ) { 01729 value = workArea[j] - thisPivotValue * multipliersL[j]; 01730 double absValue = fabs ( value ); 01731 01732 if ( absValue > tolerance ) { 01733 workArea[j] = 0.0; 01734 element[put] = value; 01735 indexRow[put] = indexL[j]; 01736 if ( absValue > largest ) { 01737 largest = absValue; 01738 positionLargest = put; 01739 } 01740 put++; 01741 } else { 01742 workArea[j] = 0.0; 01743 added--; 01744 int word = j >> COINFACTORIZATION_SHIFT_PER_INT; 01745 int bit = j & COINFACTORIZATION_MASK_PER_INT; 01746 01747 if ( temp2[word] & ( 1 << bit ) ) { 01748 //take out of row list 01749 iRow = indexL[j]; 01750 CoinBigIndex start = startRow[iRow]; 01751 CoinBigIndex end = start + numberInRow[iRow]; 01752 CoinBigIndex where = start; 01753 01754 while ( indexColumn[where] != iColumn ) { 01755 where++; 01756 } /* endwhile */ 01757 #if DEBUG_COIN 01758 if ( where >= end ) { 01759 abort ( ); 01760 } 01761 #endif 01762 indexColumn[where] = indexColumn[end - 1]; 01763 numberInRow[iRow]--; 01764 } else { 01765 //make sure won't be added 01766 int word = j >> COINFACTORIZATION_SHIFT_PER_INT; 01767 int bit = j & COINFACTORIZATION_MASK_PER_INT; 01768 01769 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts 01770 } 01771 } 01772 } 01773 numberInColumn[iColumn] = put - startColumnThis; 01774 //move largest 01775 if ( positionLargest >= 0 ) { 01776 value = element[positionLargest]; 01777 iRow = indexRow[positionLargest]; 01778 element[positionLargest] = element[startColumnThis]; 01779 indexRow[positionLargest] = indexRow[startColumnThis]; 01780 element[startColumnThis] = value; 01781 indexRow[startColumnThis] = iRow; 01782 } 01783 //linked list for column 01784 if ( nextCount[iColumn + numberRows_] != -2 ) { 01785 //modify linked list 01786 deleteLink ( iColumn + numberRows_ ); 01787 addLink ( iColumn + numberRows_, numberInColumn[iColumn] ); 01788 } 01789 temp2 += increment2; 01790 } 01791 //get space for row list 01792 unsigned int *putBase = workArea2; 01793 int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT; 01794 int i = 0; 01795 01796 // do linked lists and update counts 01797 while ( bigLoops ) { 01798 bigLoops--; 01799 int bit; 01800 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) { 01801 unsigned int *putThis = putBase; 01802 int iRow = indexL[i]; 01803 01804 //get space 01805 int number = 0; 01806 int jColumn; 01807 01808 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) { 01809 unsigned int test = *putThis; 01810 01811 putThis += increment2; 01812 test = 1 - ( ( test >> bit ) & 1 ); 01813 number += test; 01814 } 01815 int next = nextRow[iRow]; 01816 CoinBigIndex space; 01817 01818 space = startRow[next] - startRow[iRow]; 01819 number += numberInRow[iRow]; 01820 if ( space < number ) { 01821 if ( !getRowSpace ( iRow, number ) ) { 01822 goto BAD_PIVOT; 01823 } 01824 } 01825 // now do 01826 putThis = putBase; 01827 next = nextRow[iRow]; 01828 number = numberInRow[iRow]; 01829 CoinBigIndex end = startRow[iRow] + number; 01830 int saveIndex = indexColumn[startRow[next]]; 01831 01832 //add in 01833 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) { 01834 unsigned int test = *putThis; 01835 01836 putThis += increment2; 01837 test = 1 - ( ( test >> bit ) & 1 ); 01838 indexColumn[end] = saveColumn[jColumn]; 01839 end += test; 01840 } 01841 //put back next one in case zapped 01842 indexColumn[startRow[next]] = saveIndex; 01843 markRow[iRow] = FAC_UNSET; 01844 number = end - startRow[iRow]; 01845 numberInRow[iRow] = number; 01846 deleteLink ( iRow ); 01847 addLink ( iRow, number ); 01848 } 01849 putBase++; 01850 } /* endwhile */ 01851 int bit; 01852 01853 for ( bit = 0; i < numberDoColumn; i++, bit++ ) { 01854 unsigned int *putThis = putBase; 01855 int iRow = indexL[i]; 01856 01857 //get space 01858 int number = 0; 01859 int jColumn; 01860 01861 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) { 01862 unsigned int test = *putThis; 01863 01864 putThis += increment2; 01865 test = 1 - ( ( test >> bit ) & 1 ); 01866 number += test; 01867 } 01868 int next = nextRow[iRow]; 01869 CoinBigIndex space; 01870 01871 space = startRow[next] - startRow[iRow]; 01872 number += numberInRow[iRow]; 01873 if ( space < number ) { 01874 if ( !getRowSpace ( iRow, number ) ) { 01875 goto BAD_PIVOT; 01876 } 01877 } 01878 // now do 01879 putThis = putBase; 01880 next = nextRow[iRow]; 01881 number = numberInRow[iRow]; 01882 CoinBigIndex end = startRow[iRow] + number; 01883 int saveIndex; 01884 01885 saveIndex = indexColumn[startRow[next]]; 01886 01887 //add in 01888 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) { 01889 unsigned int test = *putThis; 01890 01891 putThis += increment2; 01892 test = 1 - ( ( test >> bit ) & 1 ); 01893 01894 indexColumn[end] = saveColumn[jColumn]; 01895 end += test; 01896 } 01897 indexColumn[startRow[next]] = saveIndex; 01898 markRow[iRow] = FAC_UNSET; 01899 number = end - startRow[iRow]; 01900 numberInRow[iRow] = number; 01901 deleteLink ( iRow ); 01902 addLink ( iRow, number ); 01903 } 01904 markRow[iPivotRow] = FAC_UNSET; 01905 //modify linked list for pivots 01906 deleteLink ( iPivotRow ); 01907 deleteLink ( iPivotColumn + numberRows_ ); 01908 totalElements_ += added; 01909 goodPivot= true; 01910 // **** UGLY UGLY UGLY 01911 } 01912 BAD_PIVOT: 01913 01914 ; 01915 } 01916 #undef FAC_UNSET 01917 #endif