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 if (positionLargest >= 0) 00987 positionLargest = positionLargest + startColumnU[iColumn] - startColumn; 00988 startColumn = startColumnU[iColumn]; 00989 put = startColumn + numberInColumn[iColumn]; 00990 } 00991 double tolerance = zeroTolerance_; 00992 00993 int *nextCount = nextCount_.array(); 00994 for ( j = 0; j < numberInPivotColumn; j++ ) { 00995 value = work[j] - thisPivotValue * multipliersL[j]; 00996 double absValue = fabs ( value ); 00997 00998 if ( absValue > tolerance ) { 00999 work[j] = 0.0; 01000 assert (put<lengthAreaU_); 01001 elementU[put] = value; 01002 indexRowU[put] = indexL[j]; 01003 if ( absValue > largest ) { 01004 largest = absValue; 01005 positionLargest = put; 01006 } 01007 put++; 01008 } else { 01009 work[j] = 0.0; 01010 added--; 01011 int word = j >> COINFACTORIZATION_SHIFT_PER_INT; 01012 int bit = j & COINFACTORIZATION_MASK_PER_INT; 01013 01014 if ( temp2[word] & ( 1 << bit ) ) { 01015 //take out of row list 01016 iRow = indexL[j]; 01017 CoinBigIndex start = startRowU[iRow]; 01018 CoinBigIndex end = start + numberInRow[iRow]; 01019 CoinBigIndex where = start; 01020 01021 while ( indexColumnU[where] != iColumn ) { 01022 where++; 01023 } /* endwhile */ 01024 #if DEBUG_COIN 01025 if ( where >= end ) { 01026 abort ( ); 01027 } 01028 #endif 01029 indexColumnU[where] = indexColumnU[end - 1]; 01030 numberInRow[iRow]--; 01031 } else { 01032 //make sure won't be added 01033 int word = j >> COINFACTORIZATION_SHIFT_PER_INT; 01034 int bit = j & COINFACTORIZATION_MASK_PER_INT; 01035 01036 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts 01037 } 01038 } 01039 } 01040 numberInColumn[iColumn] = put - startColumn; 01041 //move largest 01042 if ( positionLargest >= 0 ) { 01043 value = elementU[positionLargest]; 01044 iRow = indexRowU[positionLargest]; 01045 elementU[positionLargest] = elementU[startColumn]; 01046 indexRowU[positionLargest] = indexRowU[startColumn]; 01047 elementU[startColumn] = value; 01048 indexRowU[startColumn] = iRow; 01049 } 01050 //linked list for column 01051 if ( nextCount[iColumn + numberRows_] != -2 ) { 01052 //modify linked list 01053 deleteLink ( iColumn + numberRows_ ); 01054 addLink ( iColumn + numberRows_, numberInColumn[iColumn] ); 01055 } 01056 temp2 += increment2; 01057 } 01058 //get space for row list 01059 unsigned int *putBase = workArea2; 01060 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT; 01061 int i = 0; 01062 01063 // do linked lists and update counts 01064 while ( bigLoops ) { 01065 bigLoops--; 01066 int bit; 01067 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) { 01068 unsigned int *putThis = putBase; 01069 int iRow = indexL[i]; 01070 01071 //get space 01072 int number = 0; 01073 int jColumn; 01074 01075 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) { 01076 unsigned int test = *putThis; 01077 01078 putThis += increment2; 01079 test = 1 - ( ( test >> bit ) & 1 ); 01080 number += test; 01081 } 01082 int next = nextRow[iRow]; 01083 CoinBigIndex space; 01084 01085 space = startRowU[next] - startRowU[iRow]; 01086 number += numberInRow[iRow]; 01087 if ( space < number ) { 01088 if ( !getRowSpace ( iRow, number ) ) { 01089 return false; 01090 } 01091 } 01092 // now do 01093 putThis = putBase; 01094 next = nextRow[iRow]; 01095 number = numberInRow[iRow]; 01096 CoinBigIndex end = startRowU[iRow] + number; 01097 int saveIndex = indexColumnU[startRowU[next]]; 01098 01099 //add in 01100 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) { 01101 unsigned int test = *putThis; 01102 01103 putThis += increment2; 01104 test = 1 - ( ( test >> bit ) & 1 ); 01105 indexColumnU[end] = saveColumn[jColumn]; 01106 end += test; 01107 } 01108 //put back next one in case zapped 01109 indexColumnU[startRowU[next]] = saveIndex; 01110 markRow[iRow] = static_cast<T>(largeInteger+1); 01111 number = end - startRowU[iRow]; 01112 numberInRow[iRow] = number; 01113 deleteLink ( iRow ); 01114 addLink ( iRow, number ); 01115 } 01116 putBase++; 01117 } /* endwhile */ 01118 int bit; 01119 01120 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) { 01121 unsigned int *putThis = putBase; 01122 int iRow = indexL[i]; 01123 01124 //get space 01125 int number = 0; 01126 int jColumn; 01127 01128 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) { 01129 unsigned int test = *putThis; 01130 01131 putThis += increment2; 01132 test = 1 - ( ( test >> bit ) & 1 ); 01133 number += test; 01134 } 01135 int next = nextRow[iRow]; 01136 CoinBigIndex space; 01137 01138 space = startRowU[next] - startRowU[iRow]; 01139 number += numberInRow[iRow]; 01140 if ( space < number ) { 01141 if ( !getRowSpace ( iRow, number ) ) { 01142 return false; 01143 } 01144 } 01145 // now do 01146 putThis = putBase; 01147 next = nextRow[iRow]; 01148 number = numberInRow[iRow]; 01149 CoinBigIndex end = startRowU[iRow] + number; 01150 int saveIndex; 01151 01152 saveIndex = indexColumnU[startRowU[next]]; 01153 01154 //add in 01155 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) { 01156 unsigned int test = *putThis; 01157 01158 putThis += increment2; 01159 test = 1 - ( ( test >> bit ) & 1 ); 01160 01161 indexColumnU[end] = saveColumn[jColumn]; 01162 end += test; 01163 } 01164 indexColumnU[startRowU[next]] = saveIndex; 01165 markRow[iRow] = static_cast<T>(largeInteger+1); 01166 number = end - startRowU[iRow]; 01167 numberInRow[iRow] = number; 01168 deleteLink ( iRow ); 01169 addLink ( iRow, number ); 01170 } 01171 markRow[pivotRow] = static_cast<T>(largeInteger+1); 01172 //modify linked list for pivots 01173 deleteLink ( pivotRow ); 01174 deleteLink ( pivotColumn + numberRows_ ); 01175 totalElements_ += added; 01176 return true; 01177 } 01178 01179 /********************************* END LARGE TEMPLATE ********/ 01181 01182 protected: 01183 01186 01187 double pivotTolerance_; 01189 double zeroTolerance_; 01190 #ifndef COIN_FAST_CODE 01191 01192 double slackValue_; 01193 #else 01194 #ifndef slackValue_ 01195 #define slackValue_ -1.0 01196 #endif 01197 #endif 01198 01199 double areaFactor_; 01201 double relaxCheck_; 01203 int numberRows_; 01205 int numberRowsExtra_; 01207 int maximumRowsExtra_; 01209 int numberColumns_; 01211 int numberColumnsExtra_; 01213 int maximumColumnsExtra_; 01215 int numberGoodU_; 01217 int numberGoodL_; 01219 int maximumPivots_; 01221 int numberPivots_; 01224 CoinBigIndex totalElements_; 01226 CoinBigIndex factorElements_; 01228 CoinIntArrayWithLength pivotColumn_; 01230 CoinIntArrayWithLength permute_; 01232 CoinIntArrayWithLength permuteBack_; 01234 CoinIntArrayWithLength pivotColumnBack_; 01236 int status_; 01237 01242 //int increasingRows_; 01243 01245 int numberTrials_; 01247 CoinBigIndexArrayWithLength startRowU_; 01248 01250 CoinIntArrayWithLength numberInRow_; 01251 01253 CoinIntArrayWithLength numberInColumn_; 01254 01256 CoinIntArrayWithLength numberInColumnPlus_; 01257 01260 CoinIntArrayWithLength firstCount_; 01261 01263 CoinIntArrayWithLength nextCount_; 01264 01266 CoinIntArrayWithLength lastCount_; 01267 01269 CoinIntArrayWithLength nextColumn_; 01270 01272 CoinIntArrayWithLength lastColumn_; 01273 01275 CoinIntArrayWithLength nextRow_; 01276 01278 CoinIntArrayWithLength lastRow_; 01279 01281 CoinIntArrayWithLength saveColumn_; 01282 01284 CoinIntArrayWithLength markRow_; 01285 01287 int messageLevel_; 01288 01290 int biggerDimension_; 01291 01293 CoinIntArrayWithLength indexColumnU_; 01294 01296 CoinIntArrayWithLength pivotRowL_; 01297 01299 CoinFactorizationDoubleArrayWithLength pivotRegion_; 01300 01302 int numberSlacks_; 01303 01305 int numberU_; 01306 01308 CoinBigIndex maximumU_; 01309 01311 //int baseU_; 01312 01314 CoinBigIndex lengthU_; 01315 01317 CoinBigIndex lengthAreaU_; 01318 01320 CoinFactorizationDoubleArrayWithLength elementU_; 01321 01323 CoinIntArrayWithLength indexRowU_; 01324 01326 CoinBigIndexArrayWithLength startColumnU_; 01327 01329 CoinBigIndexArrayWithLength convertRowToColumnU_; 01330 01332 CoinBigIndex numberL_; 01333 01335 CoinBigIndex baseL_; 01336 01338 CoinBigIndex lengthL_; 01339 01341 CoinBigIndex lengthAreaL_; 01342 01344 CoinFactorizationDoubleArrayWithLength elementL_; 01345 01347 CoinIntArrayWithLength indexRowL_; 01348 01350 CoinBigIndexArrayWithLength startColumnL_; 01351 01353 bool doForrestTomlin_; 01354 01356 int numberR_; 01357 01359 CoinBigIndex lengthR_; 01360 01362 CoinBigIndex lengthAreaR_; 01363 01365 CoinFactorizationDouble *elementR_; 01366 01368 int *indexRowR_; 01369 01371 CoinBigIndexArrayWithLength startColumnR_; 01372 01374 double * denseArea_; 01375 01377 int * densePermute_; 01378 01380 int numberDense_; 01381 01383 int denseThreshold_; 01384 01386 CoinFactorizationDoubleArrayWithLength workArea_; 01387 01389 CoinUnsignedIntArrayWithLength workArea2_; 01390 01392 CoinBigIndex numberCompressions_; 01393 01395 mutable double ftranCountInput_; 01396 mutable double ftranCountAfterL_; 01397 mutable double ftranCountAfterR_; 01398 mutable double ftranCountAfterU_; 01399 mutable double btranCountInput_; 01400 mutable double btranCountAfterU_; 01401 mutable double btranCountAfterR_; 01402 mutable double btranCountAfterL_; 01403 01405 mutable int numberFtranCounts_; 01406 mutable int numberBtranCounts_; 01407 01409 double ftranAverageAfterL_; 01410 double ftranAverageAfterR_; 01411 double ftranAverageAfterU_; 01412 double btranAverageAfterU_; 01413 double btranAverageAfterR_; 01414 double btranAverageAfterL_; 01415 01417 mutable bool collectStatistics_; 01418 01420 int sparseThreshold_; 01421 01423 int sparseThreshold2_; 01424 01426 CoinBigIndexArrayWithLength startRowL_; 01427 01429 CoinIntArrayWithLength indexColumnL_; 01430 01432 CoinFactorizationDoubleArrayWithLength elementByRowL_; 01433 01435 mutable CoinIntArrayWithLength sparse_; 01439 int biasLU_; 01445 int persistenceFlag_; 01447 }; 01448 // Dense coding 01449 #ifdef COIN_HAS_LAPACK 01450 #define DENSE_CODE 1 01451 /* Type of Fortran integer translated into C */ 01452 #ifndef ipfint 01453 //typedef ipfint FORTRAN_INTEGER_TYPE ; 01454 typedef int ipfint; 01455 typedef const int cipfint; 01456 #endif 01457 #endif 01458 #endif 01459 // Extra for ugly include 01460 #ifdef UGLY_COIN_FACTOR_CODING 01461 #define FAC_UNSET (FAC_SET+1) 01462 { 01463 goodPivot=false; 01464 //store pivot columns (so can easily compress) 01465 CoinBigIndex startColumnThis = startColumn[iPivotColumn]; 01466 CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1; 01467 int put = 0; 01468 CoinBigIndex startRowThis = startRow[iPivotRow]; 01469 CoinBigIndex endRow = startRowThis + numberDoRow + 1; 01470 if ( pivotColumnPosition < 0 ) { 01471 for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) { 01472 int iColumn = indexColumn[pivotColumnPosition]; 01473 if ( iColumn != iPivotColumn ) { 01474 saveColumn[put++] = iColumn; 01475 } else { 01476 break; 01477 } 01478 } 01479 } else { 01480 for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) { 01481 saveColumn[put++] = indexColumn[i]; 01482 } 01483 } 01484 assert (pivotColumnPosition<endRow); 01485 assert (indexColumn[pivotColumnPosition]==iPivotColumn); 01486 pivotColumnPosition++; 01487 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) { 01488 saveColumn[put++] = indexColumn[pivotColumnPosition]; 01489 } 01490 //take out this bit of indexColumn 01491 int next = nextRow[iPivotRow]; 01492 int last = lastRow[iPivotRow]; 01493 01494 nextRow[last] = next; 01495 lastRow[next] = last; 01496 nextRow[iPivotRow] = numberGoodU_; //use for permute 01497 lastRow[iPivotRow] = -2; 01498 numberInRow[iPivotRow] = 0; 01499 //store column in L, compress in U and take column out 01500 CoinBigIndex l = lengthL_; 01501 // **** HORRID coding coming up but a goto seems best! 01502 { 01503 if ( l + numberDoColumn > lengthAreaL_ ) { 01504 //need more memory 01505 if ((messageLevel_&4)!=0) 01506 printf("more memory needed in middle of invert\n"); 01507 goto BAD_PIVOT; 01508 } 01509 //l+=currentAreaL_->elementByColumn-elementL; 01510 CoinBigIndex lSave = l; 01511 01512 CoinBigIndex * startColumnL = startColumnL_.array(); 01513 startColumnL[numberGoodL_] = l; //for luck and first time 01514 numberGoodL_++; 01515 startColumnL[numberGoodL_] = l + numberDoColumn; 01516 lengthL_ += numberDoColumn; 01517 if ( pivotRowPosition < 0 ) { 01518 for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) { 01519 int iRow = indexRow[pivotRowPosition]; 01520 if ( iRow != iPivotRow ) { 01521 indexRowL[l] = iRow; 01522 elementL[l] = element[pivotRowPosition]; 01523 markRow[iRow] = l - lSave; 01524 l++; 01525 //take out of row list 01526 CoinBigIndex start = startRow[iRow]; 01527 CoinBigIndex end = start + numberInRow[iRow]; 01528 CoinBigIndex where = start; 01529 01530 while ( indexColumn[where] != iPivotColumn ) { 01531 where++; 01532 } /* endwhile */ 01533 #if DEBUG_COIN 01534 if ( where >= end ) { 01535 abort ( ); 01536 } 01537 #endif 01538 indexColumn[where] = indexColumn[end - 1]; 01539 numberInRow[iRow]--; 01540 } else { 01541 break; 01542 } 01543 } 01544 } else { 01545 CoinBigIndex i; 01546 01547 for ( i = startColumnThis; i < pivotRowPosition; i++ ) { 01548 int iRow = indexRow[i]; 01549 01550 markRow[iRow] = l - lSave; 01551 indexRowL[l] = iRow; 01552 elementL[l] = element[i]; 01553 l++; 01554 //take out of row list 01555 CoinBigIndex start = startRow[iRow]; 01556 CoinBigIndex end = start + numberInRow[iRow]; 01557 CoinBigIndex where = start; 01558 01559 while ( indexColumn[where] != iPivotColumn ) { 01560 where++; 01561 } /* endwhile */ 01562 #if DEBUG_COIN 01563 if ( where >= end ) { 01564 abort ( ); 01565 } 01566 #endif 01567 indexColumn[where] = indexColumn[end - 1]; 01568 numberInRow[iRow]--; 01569 assert (numberInRow[iRow]>=0); 01570 } 01571 } 01572 assert (pivotRowPosition<endColumn); 01573 assert (indexRow[pivotRowPosition]==iPivotRow); 01574 CoinFactorizationDouble pivotElement = element[pivotRowPosition]; 01575 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement; 01576 01577 pivotRegion_.array()[numberGoodU_] = pivotMultiplier; 01578 pivotRowPosition++; 01579 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) { 01580 int iRow = indexRow[pivotRowPosition]; 01581 01582 markRow[iRow] = l - lSave; 01583 indexRowL[l] = iRow; 01584 elementL[l] = element[pivotRowPosition]; 01585 l++; 01586 //take out of row list 01587 CoinBigIndex start = startRow[iRow]; 01588 CoinBigIndex end = start + numberInRow[iRow]; 01589 CoinBigIndex where = start; 01590 01591 while ( indexColumn[where] != iPivotColumn ) { 01592 where++; 01593 } /* endwhile */ 01594 #if DEBUG_COIN 01595 if ( where >= end ) { 01596 abort ( ); 01597 } 01598 #endif 01599 indexColumn[where] = indexColumn[end - 1]; 01600 numberInRow[iRow]--; 01601 assert (numberInRow[iRow]>=0); 01602 } 01603 markRow[iPivotRow] = FAC_SET; 01604 //compress pivot column (move pivot to front including saved) 01605 numberInColumn[iPivotColumn] = 0; 01606 //use end of L for temporary space 01607 int *indexL = &indexRowL[lSave]; 01608 CoinFactorizationDouble *multipliersL = &elementL[lSave]; 01609 01610 //adjust 01611 int j; 01612 01613 for ( j = 0; j < numberDoColumn; j++ ) { 01614 multipliersL[j] *= pivotMultiplier; 01615 } 01616 //zero out fill 01617 CoinBigIndex iErase; 01618 for ( iErase = 0; iErase < increment2 * numberDoRow; 01619 iErase++ ) { 01620 workArea2[iErase] = 0; 01621 } 01622 CoinBigIndex added = numberDoRow * numberDoColumn; 01623 unsigned int *temp2 = workArea2; 01624 int * nextColumn = nextColumn_.array(); 01625 01626 //pack down and move to work 01627 int jColumn; 01628 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) { 01629 int iColumn = saveColumn[jColumn]; 01630 CoinBigIndex startColumnThis = startColumn[iColumn]; 01631 CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn]; 01632 int iRow = indexRow[startColumnThis]; 01633 CoinFactorizationDouble value = element[startColumnThis]; 01634 double largest; 01635 CoinBigIndex put = startColumnThis; 01636 CoinBigIndex positionLargest = -1; 01637 CoinFactorizationDouble thisPivotValue = 0.0; 01638 01639 //compress column and find largest not updated 01640 bool checkLargest; 01641 int mark = markRow[iRow]; 01642 01643 if ( mark == FAC_UNSET ) { 01644 largest = fabs ( value ); 01645 positionLargest = put; 01646 put++; 01647 checkLargest = false; 01648 } else { 01649 //need to find largest 01650 largest = 0.0; 01651 checkLargest = true; 01652 if ( mark != FAC_SET ) { 01653 //will be updated 01654 workArea[mark] = value; 01655 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT; 01656 int bit = mark & COINFACTORIZATION_MASK_PER_INT; 01657 01658 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts 01659 added--; 01660 } else { 01661 thisPivotValue = value; 01662 } 01663 } 01664 CoinBigIndex i; 01665 for ( i = startColumnThis + 1; i < endColumn; i++ ) { 01666 iRow = indexRow[i]; 01667 value = element[i]; 01668 int mark = markRow[iRow]; 01669 01670 if ( mark == FAC_UNSET ) { 01671 //keep 01672 indexRow[put] = iRow; 01673 element[put] = value; 01674 if ( checkLargest ) { 01675 double absValue = fabs ( value ); 01676 01677 if ( absValue > largest ) { 01678 largest = absValue; 01679 positionLargest = put; 01680 } 01681 } 01682 put++; 01683 } else if ( mark != FAC_SET ) { 01684 //will be updated 01685 workArea[mark] = value; 01686 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT; 01687 int bit = mark & COINFACTORIZATION_MASK_PER_INT; 01688 01689 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts 01690 added--; 01691 } else { 01692 thisPivotValue = value; 01693 } 01694 } 01695 //slot in pivot 01696 element[put] = element[startColumnThis]; 01697 indexRow[put] = indexRow[startColumnThis]; 01698 if ( positionLargest == startColumnThis ) { 01699 positionLargest = put; //follow if was largest 01700 } 01701 put++; 01702 element[startColumnThis] = thisPivotValue; 01703 indexRow[startColumnThis] = iPivotRow; 01704 //clean up counts 01705 startColumnThis++; 01706 numberInColumn[iColumn] = put - startColumnThis; 01707 int * numberInColumnPlus = numberInColumnPlus_.array(); 01708 numberInColumnPlus[iColumn]++; 01709 startColumn[iColumn]++; 01710 //how much space have we got 01711 int next = nextColumn[iColumn]; 01712 CoinBigIndex space; 01713 01714 space = startColumn[next] - put - numberInColumnPlus[next]; 01715 //assume no zero elements 01716 if ( numberDoColumn > space ) { 01717 //getColumnSpace also moves fixed part 01718 if ( !getColumnSpace ( iColumn, numberDoColumn ) ) { 01719 goto BAD_PIVOT; 01720 } 01721 //redo starts 01722 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis; 01723 startColumnThis = startColumn[iColumn]; 01724 put = startColumnThis + numberInColumn[iColumn]; 01725 } 01726 double tolerance = zeroTolerance_; 01727 01728 int *nextCount = nextCount_.array(); 01729 for ( j = 0; j < numberDoColumn; j++ ) { 01730 value = workArea[j] - thisPivotValue * multipliersL[j]; 01731 double absValue = fabs ( value ); 01732 01733 if ( absValue > tolerance ) { 01734 workArea[j] = 0.0; 01735 element[put] = value; 01736 indexRow[put] = indexL[j]; 01737 if ( absValue > largest ) { 01738 largest = absValue; 01739 positionLargest = put; 01740 } 01741 put++; 01742 } else { 01743 workArea[j] = 0.0; 01744 added--; 01745 int word = j >> COINFACTORIZATION_SHIFT_PER_INT; 01746 int bit = j & COINFACTORIZATION_MASK_PER_INT; 01747 01748 if ( temp2[word] & ( 1 << bit ) ) { 01749 //take out of row list 01750 iRow = indexL[j]; 01751 CoinBigIndex start = startRow[iRow]; 01752 CoinBigIndex end = start + numberInRow[iRow]; 01753 CoinBigIndex where = start; 01754 01755 while ( indexColumn[where] != iColumn ) { 01756 where++; 01757 } /* endwhile */ 01758 #if DEBUG_COIN 01759 if ( where >= end ) { 01760 abort ( ); 01761 } 01762 #endif 01763 indexColumn[where] = indexColumn[end - 1]; 01764 numberInRow[iRow]--; 01765 } else { 01766 //make sure won't be added 01767 int word = j >> COINFACTORIZATION_SHIFT_PER_INT; 01768 int bit = j & COINFACTORIZATION_MASK_PER_INT; 01769 01770 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts 01771 } 01772 } 01773 } 01774 numberInColumn[iColumn] = put - startColumnThis; 01775 //move largest 01776 if ( positionLargest >= 0 ) { 01777 value = element[positionLargest]; 01778 iRow = indexRow[positionLargest]; 01779 element[positionLargest] = element[startColumnThis]; 01780 indexRow[positionLargest] = indexRow[startColumnThis]; 01781 element[startColumnThis] = value; 01782 indexRow[startColumnThis] = iRow; 01783 } 01784 //linked list for column 01785 if ( nextCount[iColumn + numberRows_] != -2 ) { 01786 //modify linked list 01787 deleteLink ( iColumn + numberRows_ ); 01788 addLink ( iColumn + numberRows_, numberInColumn[iColumn] ); 01789 } 01790 temp2 += increment2; 01791 } 01792 //get space for row list 01793 unsigned int *putBase = workArea2; 01794 int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT; 01795 int i = 0; 01796 01797 // do linked lists and update counts 01798 while ( bigLoops ) { 01799 bigLoops--; 01800 int bit; 01801 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) { 01802 unsigned int *putThis = putBase; 01803 int iRow = indexL[i]; 01804 01805 //get space 01806 int number = 0; 01807 int jColumn; 01808 01809 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) { 01810 unsigned int test = *putThis; 01811 01812 putThis += increment2; 01813 test = 1 - ( ( test >> bit ) & 1 ); 01814 number += test; 01815 } 01816 int next = nextRow[iRow]; 01817 CoinBigIndex space; 01818 01819 space = startRow[next] - startRow[iRow]; 01820 number += numberInRow[iRow]; 01821 if ( space < number ) { 01822 if ( !getRowSpace ( iRow, number ) ) { 01823 goto BAD_PIVOT; 01824 } 01825 } 01826 // now do 01827 putThis = putBase; 01828 next = nextRow[iRow]; 01829 number = numberInRow[iRow]; 01830 CoinBigIndex end = startRow[iRow] + number; 01831 int saveIndex = indexColumn[startRow[next]]; 01832 01833 //add in 01834 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) { 01835 unsigned int test = *putThis; 01836 01837 putThis += increment2; 01838 test = 1 - ( ( test >> bit ) & 1 ); 01839 indexColumn[end] = saveColumn[jColumn]; 01840 end += test; 01841 } 01842 //put back next one in case zapped 01843 indexColumn[startRow[next]] = saveIndex; 01844 markRow[iRow] = FAC_UNSET; 01845 number = end - startRow[iRow]; 01846 numberInRow[iRow] = number; 01847 deleteLink ( iRow ); 01848 addLink ( iRow, number ); 01849 } 01850 putBase++; 01851 } /* endwhile */ 01852 int bit; 01853 01854 for ( bit = 0; i < numberDoColumn; i++, bit++ ) { 01855 unsigned int *putThis = putBase; 01856 int iRow = indexL[i]; 01857 01858 //get space 01859 int number = 0; 01860 int jColumn; 01861 01862 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) { 01863 unsigned int test = *putThis; 01864 01865 putThis += increment2; 01866 test = 1 - ( ( test >> bit ) & 1 ); 01867 number += test; 01868 } 01869 int next = nextRow[iRow]; 01870 CoinBigIndex space; 01871 01872 space = startRow[next] - startRow[iRow]; 01873 number += numberInRow[iRow]; 01874 if ( space < number ) { 01875 if ( !getRowSpace ( iRow, number ) ) { 01876 goto BAD_PIVOT; 01877 } 01878 } 01879 // now do 01880 putThis = putBase; 01881 next = nextRow[iRow]; 01882 number = numberInRow[iRow]; 01883 CoinBigIndex end = startRow[iRow] + number; 01884 int saveIndex; 01885 01886 saveIndex = indexColumn[startRow[next]]; 01887 01888 //add in 01889 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) { 01890 unsigned int test = *putThis; 01891 01892 putThis += increment2; 01893 test = 1 - ( ( test >> bit ) & 1 ); 01894 01895 indexColumn[end] = saveColumn[jColumn]; 01896 end += test; 01897 } 01898 indexColumn[startRow[next]] = saveIndex; 01899 markRow[iRow] = FAC_UNSET; 01900 number = end - startRow[iRow]; 01901 numberInRow[iRow] = number; 01902 deleteLink ( iRow ); 01903 addLink ( iRow, number ); 01904 } 01905 markRow[iPivotRow] = FAC_UNSET; 01906 //modify linked list for pivots 01907 deleteLink ( iPivotRow ); 01908 deleteLink ( iPivotColumn + numberRows_ ); 01909 totalElements_ += added; 01910 goodPivot= true; 01911 // **** UGLY UGLY UGLY 01912 } 01913 BAD_PIVOT: 01914 01915 ; 01916 } 01917 #undef FAC_UNSET 01918 #endif