- Main Page
- Modules
- Namespaces
- Classes
- Files
- Related Pages

- File List
- File Members

# /home/coin/SVN-release/Clp-1.5.0/CoinUtils/src/CoinFactorization.hpp

Go to the documentation of this file.

    00001 // Copyright (C) 2002, International Business Machines
    00002 // Corporation and others.  All Rights Reserved.
    00003
    00004 /*
    00005    Authors
    00006
    00007    John Forrest
    00008
    00009  */
    00010 #ifndef CoinFactorization_H
    00011 #define CoinFactorization_H
    00012
    00013 #include <iostream>
    00014 #include <string>
    00015 #include <cassert>
    00016 #include "CoinFinite.hpp"
    00017 #include "CoinIndexedVector.hpp"
    00018 class CoinPackedMatrix;
    00045 class CoinFactorization {
    00046    friend void CoinFactorizationUnitTest( const std::string & mpsDir );
    00047
    00048 public:
    00049
    00052
    00053     CoinFactorization (  );
    00055   CoinFactorization ( const CoinFactorization &other);
    00056
    00058    ~CoinFactorization (  );
    00060   void almostDestructor();
    00062   void show_self (  ) const;
    00064   int saveFactorization (const char * file  ) const;
    00068   int restoreFactorization (const char * file  , bool factor=false) ;
    00070   void sort (  ) const;
    00072     CoinFactorization & operator = ( const CoinFactorization & other );
    00074
    00084   int factorize ( const CoinPackedMatrix & matrix,
    00085                   int rowIsBasic[], int columnIsBasic[] ,
    00086                   double areaFactor = 0.0 );
    00097   int factorize ( int numberRows,
    00098                   int numberColumns,
    00099                   CoinBigIndex numberElements,
    00100                   CoinBigIndex maximumL,
    00101                   CoinBigIndex maximumU,
    00102                   const int indicesRow[],
    00103                   const int indicesColumn[], const double elements[] ,
    00104                   int permutation[],
    00105                   double areaFactor = 0.0);
    00110   int factorizePart1 ( int numberRows,
    00111                        int numberColumns,
    00112                        CoinBigIndex estimateNumberElements,
    00113                        int * indicesRow[],
    00114                        int * indicesColumn[],
    00115                        double * elements[],
    00116                   double areaFactor = 0.0);
    00123   int factorizePart2 (int permutation[],int exactNumberElements);
    00125
    00128
    00129   inline int status (  ) const {
    00130     return status_;
    00131   }
    00133   inline int pivots (  ) const {
    00134     return numberPivots_;
    00135   }
    00137   inline int *permute (  ) const {
    00138     return permute_.array();
    00139   }
    00141   inline int *pivotColumn (  ) const {
    00142     return pivotColumn_.array();
    00143   }
    00145   inline int *permuteBack (  ) const {
    00146     return permuteBack_.array();
    00147   }
    00149   inline int *pivotColumnBack (  ) const {
    00150     return pivotColumnBack_.array();
    00151   }
    00153   inline int numberRowsExtra (  ) const {
    00154     return numberRowsExtra_;
    00155   }
    00157   inline int numberRows (  ) const {
    00158     return numberRows_;
    00159   }
    00161   inline int maximumRowsExtra (  ) const {
    00162     return maximumRowsExtra_;
    00163   }
    00165   inline int numberColumns (  ) const {
    00166     return numberColumns_;
    00167   }
    00169   inline int numberElements (  ) const {
    00170     return totalElements_;
    00171   }
    00173   inline int numberForrestTomlin (  ) const {
    00174     return numberInColumn_.array()[numberColumnsExtra_];
    00175   }
    00177   inline int numberGoodColumns (  ) const {
    00178     return numberGoodU_;
    00179   }
    00181   inline double areaFactor (  ) const {
    00182     return areaFactor_;
    00183   }
    00184   inline void areaFactor ( double value ) {
    00185     areaFactor_=value;
    00186   }
    00188   double adjustedAreaFactor() const;
    00190   inline void relaxAccuracyCheck(double value)
    00191   { relaxCheck_ = value;}
    00192   inline double getAccuracyCheck() const
    00193   { return relaxCheck_;}
    00195   inline bool increasingRows (  ) const
    00196   { return true; }
    00202   inline void increasingRows ( int value  ) {}
    00204   inline int messageLevel (  ) const {
    00205     return messageLevel_ ;
    00206   }
    00207   void messageLevel (  int value );
    00209   inline int maximumPivots (  ) const {
    00210     return maximumPivots_ ;
    00211   }
    00212   void maximumPivots (  int value );
    00213
    00215   inline int denseThreshold() const
    00216     { return denseThreshold_;}
    00218   inline void setDenseThreshold(int value)
    00219     { denseThreshold_ = value;}
    00221   inline double pivotTolerance (  ) const {
    00222     return pivotTolerance_ ;
    00223   }
    00224   void pivotTolerance (  double value );
    00226   inline double zeroTolerance (  ) const {
    00227     return zeroTolerance_ ;
    00228   }
    00229   void zeroTolerance (  double value );
    00231   inline double slackValue (  ) const {
    00232     return slackValue_ ;
    00233   }
    00234   void slackValue (  double value );
    00236   double maximumCoefficient() const;
    00238   inline bool forrestTomlin() const
    00239   { return doForrestTomlin_;}
    00240   inline void setForrestTomlin(bool value)
    00241   { doForrestTomlin_=value;}
    00243
    00246
    00248   inline int numberDense() const
    00249   { return numberDense_;}
    00250
    00252   inline CoinBigIndex numberElementsU (  ) const {
    00253     return lengthU_;
    00254   }
    00256   inline CoinBigIndex lengthAreaU (  ) const {
    00257     return lengthAreaU_;
    00258   }
    00260   inline CoinBigIndex numberElementsL (  ) const {
    00261     return lengthL_;
    00262   }
    00264   inline CoinBigIndex lengthAreaL (  ) const {
    00265     return lengthAreaL_;
    00266   }
    00268   inline CoinBigIndex numberElementsR (  ) const {
    00269     return lengthR_;
    00270   }
    00272   inline CoinBigIndex numberCompressions() const
    00273   { return numberCompressions_;}
    00277   inline int biasLU() const
    00278   { return biasLU_;}
    00279   inline void setBiasLU(int value)
    00280   { biasLU_=value;}
    00286   inline int persistenceFlag() const
    00287   { return persistenceFlag_;}
    00288   void setPersistenceFlag(int value);
    00290
    00293
    00301   int replaceColumn ( CoinIndexedVector * regionSparse,
    00302                       int pivotRow,
    00303                       double pivotCheck ,
    00304                       bool checkBeforeModifying=false);
    00306
    00316   int updateColumnFT ( CoinIndexedVector * regionSparse,
    00317                        CoinIndexedVector * regionSparse2);
    00320   int updateColumn ( CoinIndexedVector * regionSparse,
    00321                      CoinIndexedVector * regionSparse2,
    00322                      bool noPermute=false) const;
    00327   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
    00328                               CoinIndexedVector * regionSparse2) const;
    00330   void goSparse();
    00332   int sparseThreshold ( ) const;
    00334   void sparseThreshold ( int value );
    00336
    00337
    00341
    00342   inline void clearArrays()
    00343   { gutsOfDestructor();}
    00345
    00348
    00351   int add ( CoinBigIndex numberElements,
    00352                int indicesRow[],
    00353                int indicesColumn[], double elements[] );
    00354
    00357   int addColumn ( CoinBigIndex numberElements,
    00358                      int indicesRow[], double elements[] );
    00359
    00362   int addRow ( CoinBigIndex numberElements,
    00363                   int indicesColumn[], double elements[] );
    00364
    00366   int deleteColumn ( int Row );
    00368   int deleteRow ( int Row );
    00369
    00373   int replaceRow ( int whichRow, int numberElements,
    00374                       const int indicesColumn[], const double elements[] );
    00376   void emptyRows(int numberToEmpty, const int which[]);
    00378 protected:
    00379
    00381
    00382   void getAreas ( int numberRows,
    00383                   int numberColumns,
    00384                   CoinBigIndex maximumL,
    00385                   CoinBigIndex maximumU );
    00386
    00389   void preProcess ( int state,
    00390                     int possibleDuplicates = -1 );
    00392   int factor (  );
    00395   int factorSparse (  );
    00398   int factorDense (  );
    00399
    00401   bool pivotOneOtherRow ( int pivotRow,
    00402                           int pivotColumn );
    00404   bool pivotRowSingleton ( int pivotRow,
    00405                            int pivotColumn );
    00407   bool pivotColumnSingleton ( int pivotRow,
    00408                               int pivotColumn );
    00409
    00414   bool getColumnSpace ( int iColumn,
    00415                         int extraNeeded );
    00416
    00420   bool getColumnSpaceIterateR ( int iColumn, double value,
    00421                                int iRow);
    00427   CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
    00428                                int iRow);
    00432   bool getRowSpace ( int iRow, int extraNeeded );
    00433
    00437   bool getRowSpaceIterate ( int iRow,
    00438                             int extraNeeded );
    00440   void checkConsistency (  );
    00442   inline void addLink ( int index, int count ) {
    00443     int *nextCount = nextCount_.array();
    00444     int *firstCount = firstCount_.array();
    00445     int *lastCount = lastCount_.array();
    00446     int next = firstCount[count];
    00447       lastCount[index] = -2 - count;
    00448     if ( next < 0 ) {
    00449       //first with that count
    00450       firstCount[count] = index;
    00451       nextCount[index] = -1;
    00452     } else {
    00453       firstCount[count] = index;
    00454       nextCount[index] = next;
    00455       lastCount[next] = index;
    00456   }}
    00458   inline void deleteLink ( int index ) {
    00459     int *nextCount = nextCount_.array();
    00460     int *firstCount = firstCount_.array();
    00461     int *lastCount = lastCount_.array();
    00462     int next = nextCount[index];
    00463     int last = lastCount[index];
    00464     if ( last >= 0 ) {
    00465       nextCount[last] = next;
    00466     } else {
    00467       int count = -last - 2;
    00468
    00469       firstCount[count] = next;
    00470     }
    00471     if ( next >= 0 ) {
    00472       lastCount[next] = last;
    00473     }
    00474     nextCount[index] = -2;
    00475     lastCount[index] = -2;
    00476     return;
    00477   }
    00479   void separateLinks(int count,bool rowsFirst);
    00481   void cleanup (  );
    00482
    00484   void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
    00486   void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
    00488   void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
    00490   void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
    00491
    00493   void updateColumnR ( CoinIndexedVector * region ) const;
    00496   void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
    00497
    00499   void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
    00500
    00502   void updateColumnUSparse ( CoinIndexedVector * regionSparse,
    00503                              int * indexIn) const;
    00505   void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
    00506                                int * indexIn) const;
    00508   void updateColumnUDensish ( CoinIndexedVector * regionSparse,
    00509                               int * indexIn) const;
    00511   void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
    00513   void permuteBack ( CoinIndexedVector * regionSparse,
    00514                      CoinIndexedVector * outVector) const;
    00515
    00517   void updateColumnTransposePFI ( CoinIndexedVector * region) const;
    00520   void updateColumnTransposeU ( CoinIndexedVector * region,
    00521                                 int smallestIndex) const;
    00524   void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
    00525                                         int smallestIndex) const;
    00528   void updateColumnTransposeUDensish ( CoinIndexedVector * region,
    00529                                        int smallestIndex) const;
    00532   void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
    00533
    00535   void updateColumnTransposeR ( CoinIndexedVector * region ) const;
    00537   void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
    00539   void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
    00540
    00542   void updateColumnTransposeL ( CoinIndexedVector * region ) const;
    00544   void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
    00546   void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
    00548   void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
    00550   void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
    00555   int replaceColumnPFI ( CoinIndexedVector * regionSparse,
    00556                          int pivotRow, double alpha);
    00557
    00560   int checkPivot(double saveFromU, double oldPivot) const;
    00562   void gutsOfDestructor(int type=1);
    00564   void gutsOfInitialize(int type);
    00565   void gutsOfCopy(const CoinFactorization &other);
    00566
    00568   void resetStatistics();
    00569
    00570   /********************************* START LARGE TEMPLATE ********/
    00571 #ifdef INT_IS_8
    00572 #define COINFACTORIZATION_BITS_PER_INT 64
    00573 #define COINFACTORIZATION_SHIFT_PER_INT 6
    00574 #define COINFACTORIZATION_MASK_PER_INT 0x3f
    00575 #else
    00576 #define COINFACTORIZATION_BITS_PER_INT 32
    00577 #define COINFACTORIZATION_SHIFT_PER_INT 5
    00578 #define COINFACTORIZATION_MASK_PER_INT 0x1f
    00579 #endif
    00580   template <class T>  inline bool
    00581   pivot ( int pivotRow,
    00582           int pivotColumn,
    00583           CoinBigIndex pivotRowPosition,
    00584           CoinBigIndex pivotColumnPosition,
    00585           double work[],
    00586           unsigned int workArea2[],
    00587           int increment,
    00588           int increment2,
    00589           T markRow[] ,
    00590           int largeInteger)
    00591 {
    00592   int *indexColumnU = indexColumnU_.array();
    00593   CoinBigIndex *startColumnU = startColumnU_.array();
    00594   int *numberInColumn = numberInColumn_.array();
    00595   double *elementU = elementU_.array();
    00596   int *indexRowU = indexRowU_.array();
    00597   CoinBigIndex *startRowU = startRowU_.array();
    00598   int *numberInRow = numberInRow_.array();
    00599   double *elementL = elementL_.array();
    00600   int *indexRowL = indexRowL_.array();
    00601   int *saveColumn = saveColumn_.array();
    00602   int *nextRow = nextRow_.array();
    00603   int *lastRow = lastRow_.array() ;
    00604
    00605   //store pivot columns (so can easily compress)
    00606   int numberInPivotRow = numberInRow[pivotRow] - 1;
    00607   CoinBigIndex startColumn = startColumnU[pivotColumn];
    00608   int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
    00609   CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
    00610   int put = 0;
    00611   CoinBigIndex startRow = startRowU[pivotRow];
    00612   CoinBigIndex endRow = startRow + numberInPivotRow + 1;
    00613
    00614   if ( pivotColumnPosition < 0 ) {
    00615     for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    00616       int iColumn = indexColumnU[pivotColumnPosition];
    00617       if ( iColumn != pivotColumn ) {
    00618         saveColumn[put++] = iColumn;
    00619       } else {
    00620         break;
    00621       }
    00622     }
    00623   } else {
    00624     for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
    00625       saveColumn[put++] = indexColumnU[i];
    00626     }
    00627   }
    00628   assert (pivotColumnPosition<endRow);
    00629   assert (indexColumnU[pivotColumnPosition]==pivotColumn);
    00630   pivotColumnPosition++;
    00631   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    00632     saveColumn[put++] = indexColumnU[pivotColumnPosition];
    00633   }
    00634   //take out this bit of indexColumnU
    00635   int next = nextRow[pivotRow];
    00636   int last = lastRow[pivotRow];
    00637
    00638   nextRow[last] = next;
    00639   lastRow[next] = last;
    00640   nextRow[pivotRow] = numberGoodU_;     //use for permute
    00641   lastRow[pivotRow] = -2;
    00642   numberInRow[pivotRow] = 0;
    00643   //store column in L, compress in U and take column out
    00644   CoinBigIndex l = lengthL_;
    00645
    00646   if ( l + numberInPivotColumn > lengthAreaL_ ) {
    00647     //need more memory
    00648     printf("more memory needed in middle of invert\n");
    00649     return false;
    00650   }
    00651   //l+=currentAreaL_->elementByColumn-elementL;
    00652   CoinBigIndex lSave = l;
    00653
    00654   pivotRowL_.array()[numberGoodL_] = pivotRow;
    00655   CoinBigIndex * startColumnL = startColumnL_.array();
    00656   startColumnL[numberGoodL_] = l;       //for luck and first time
    00657   numberGoodL_++;
    00658   startColumnL[numberGoodL_] = l + numberInPivotColumn;
    00659   lengthL_ += numberInPivotColumn;
    00660   if ( pivotRowPosition < 0 ) {
    00661     for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    00662       int iRow = indexRowU[pivotRowPosition];
    00663       if ( iRow != pivotRow ) {
    00664         indexRowL[l] = iRow;
    00665         elementL[l] = elementU[pivotRowPosition];
    00666         markRow[iRow] = l - lSave;
    00667         l++;
    00668         //take out of row list
    00669         CoinBigIndex start = startRowU[iRow];
    00670         CoinBigIndex end = start + numberInRow[iRow];
    00671         CoinBigIndex where = start;
    00672
    00673         while ( indexColumnU[where] != pivotColumn ) {
    00674           where++;
    00675         }                       /* endwhile */
    00676 #if DEBUG_COIN
    00677         if ( where >= end ) {
    00678           abort (  );
    00679         }
    00680 #endif
    00681         indexColumnU[where] = indexColumnU[end - 1];
    00682         numberInRow[iRow]--;
    00683       } else {
    00684         break;
    00685       }
    00686     }
    00687   } else {
    00688     CoinBigIndex i;
    00689
    00690     for ( i = startColumn; i < pivotRowPosition; i++ ) {
    00691       int iRow = indexRowU[i];
    00692
    00693       markRow[iRow] = l - lSave;
    00694       indexRowL[l] = iRow;
    00695       elementL[l] = elementU[i];
    00696       l++;
    00697       //take out of row list
    00698       CoinBigIndex start = startRowU[iRow];
    00699       CoinBigIndex end = start + numberInRow[iRow];
    00700       CoinBigIndex where = start;
    00701
    00702       while ( indexColumnU[where] != pivotColumn ) {
    00703         where++;
    00704       }                         /* endwhile */
    00705 #if DEBUG_COIN
    00706       if ( where >= end ) {
    00707         abort (  );
    00708       }
    00709 #endif
    00710       indexColumnU[where] = indexColumnU[end - 1];
    00711       numberInRow[iRow]--;
    00712       assert (numberInRow[iRow]>=0);
    00713     }
    00714   }
    00715   assert (pivotRowPosition<endColumn);
    00716   assert (indexRowU[pivotRowPosition]==pivotRow);
    00717   double pivotElement = elementU[pivotRowPosition];
    00718   double pivotMultiplier = 1.0 / pivotElement;
    00719
    00720   pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
    00721   pivotRowPosition++;
    00722   for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    00723     int iRow = indexRowU[pivotRowPosition];
    00724
    00725     markRow[iRow] = l - lSave;
    00726     indexRowL[l] = iRow;
    00727     elementL[l] = elementU[pivotRowPosition];
    00728     l++;
    00729     //take out of row list
    00730     CoinBigIndex start = startRowU[iRow];
    00731     CoinBigIndex end = start + numberInRow[iRow];
    00732     CoinBigIndex where = start;
    00733
    00734     while ( indexColumnU[where] != pivotColumn ) {
    00735       where++;
    00736     }                           /* endwhile */
    00737 #if DEBUG_COIN
    00738     if ( where >= end ) {
    00739       abort (  );
    00740     }
    00741 #endif
    00742     indexColumnU[where] = indexColumnU[end - 1];
    00743     numberInRow[iRow]--;
    00744     assert (numberInRow[iRow]>=0);
    00745   }
    00746   markRow[pivotRow] = largeInteger;
    00747   //compress pivot column (move pivot to front including saved)
    00748   numberInColumn[pivotColumn] = 0;
    00749   //use end of L for temporary space
    00750   int *indexL = &indexRowL[lSave];
    00751   double *multipliersL = &elementL[lSave];
    00752
    00753   //adjust
    00754   int j;
    00755
    00756   for ( j = 0; j < numberInPivotColumn; j++ ) {
    00757     multipliersL[j] *= pivotMultiplier;
    00758   }
    00759   //zero out fill
    00760   CoinBigIndex iErase;
    00761   for ( iErase = 0; iErase < increment2 * numberInPivotRow;
    00762         iErase++ ) {
    00763     workArea2[iErase] = 0;
    00764   }
    00765   CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
    00766   unsigned int *temp2 = workArea2;
    00767   int * nextColumn = nextColumn_.array();
    00768
    00769   //pack down and move to work
    00770   int jColumn;
    00771   for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    00772     int iColumn = saveColumn[jColumn];
    00773     CoinBigIndex startColumn = startColumnU[iColumn];
    00774     CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
    00775     int iRow = indexRowU[startColumn];
    00776     double value = elementU[startColumn];
    00777     double largest;
    00778     CoinBigIndex put = startColumn;
    00779     CoinBigIndex positionLargest = -1;
    00780     double thisPivotValue = 0.0;
    00781
    00782     //compress column and find largest not updated
    00783     bool checkLargest;
    00784     int mark = markRow[iRow];
    00785
    00786     if ( mark < 0 ) {
    00787       largest = fabs ( value );
    00788       positionLargest = put;
    00789       put++;
    00790       checkLargest = false;
    00791     } else {
    00792       //need to find largest
    00793       largest = 0.0;
    00794       checkLargest = true;
    00795       if ( mark != largeInteger ) {
    00796         //will be updated
    00797         work[mark] = value;
    00798         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    00799         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    00800
    00801         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
    00802         added--;
    00803       } else {
    00804         thisPivotValue = value;
    00805       }
    00806     }
    00807     CoinBigIndex i;
    00808     for ( i = startColumn + 1; i < endColumn; i++ ) {
    00809       iRow = indexRowU[i];
    00810       value = elementU[i];
    00811       int mark = markRow[iRow];
    00812
    00813       if ( mark < 0 ) {
    00814         //keep
    00815         indexRowU[put] = iRow;
    00816         elementU[put] = value;;
    00817         if ( checkLargest ) {
    00818           double absValue = fabs ( value );
    00819
    00820           if ( absValue > largest ) {
    00821             largest = absValue;
    00822             positionLargest = put;
    00823           }
    00824         }
    00825         put++;
    00826       } else if ( mark != largeInteger ) {
    00827         //will be updated
    00828         work[mark] = value;;
    00829         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    00830         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    00831
    00832         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
    00833         added--;
    00834       } else {
    00835         thisPivotValue = value;
    00836       }
    00837     }
    00838     //slot in pivot
    00839     elementU[put] = elementU[startColumn];
    00840     indexRowU[put] = indexRowU[startColumn];
    00841     if ( positionLargest == startColumn ) {
    00842       positionLargest = put;    //follow if was largest
    00843     }
    00844     put++;
    00845     elementU[startColumn] = thisPivotValue;
    00846     indexRowU[startColumn] = pivotRow;
    00847     //clean up counts
    00848     startColumn++;
    00849     numberInColumn[iColumn] = put - startColumn;
    00850     int * numberInColumnPlus = numberInColumnPlus_.array();
    00851     numberInColumnPlus[iColumn]++;
    00852     startColumnU[iColumn]++;
    00853     //how much space have we got
    00854     int next = nextColumn[iColumn];
    00855     CoinBigIndex space;
    00856
    00857     space = startColumnU[next] - put - numberInColumnPlus[next];
    00858     //assume no zero elements
    00859     if ( numberInPivotColumn > space ) {
    00860       //getColumnSpace also moves fixed part
    00861       if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
    00862         return false;
    00863       }
    00864       //redo starts
    00865       positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
    00866       startColumn = startColumnU[iColumn];
    00867       put = startColumn + numberInColumn[iColumn];
    00868     }
    00869     double tolerance = zeroTolerance_;
    00870
    00871     int *nextCount = nextCount_.array();
    00872     for ( j = 0; j < numberInPivotColumn; j++ ) {
    00873       value = work[j] - thisPivotValue * multipliersL[j];
    00874       double absValue = fabs ( value );
    00875
    00876       if ( absValue > tolerance ) {
    00877         work[j] = 0.0;
    00878         elementU[put] = value;
    00879         indexRowU[put] = indexL[j];
    00880         if ( absValue > largest ) {
    00881           largest = absValue;
    00882           positionLargest = put;
    00883         }
    00884         put++;
    00885       } else {
    00886         work[j] = 0.0;
    00887         added--;
    00888         int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    00889         int bit = j & COINFACTORIZATION_MASK_PER_INT;
    00890
    00891         if ( temp2[word] & ( 1 << bit ) ) {
    00892           //take out of row list
    00893           iRow = indexL[j];
    00894           CoinBigIndex start = startRowU[iRow];
    00895           CoinBigIndex end = start + numberInRow[iRow];
    00896           CoinBigIndex where = start;
    00897
    00898           while ( indexColumnU[where] != iColumn ) {
    00899             where++;
    00900           }                     /* endwhile */
    00901 #if DEBUG_COIN
    00902           if ( where >= end ) {
    00903             abort (  );
    00904           }
    00905 #endif
    00906           indexColumnU[where] = indexColumnU[end - 1];
    00907           numberInRow[iRow]--;
    00908         } else {
    00909           //make sure won't be added
    00910           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    00911           int bit = j & COINFACTORIZATION_MASK_PER_INT;
    00912
    00913           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    00914         }
    00915       }
    00916     }
    00917     numberInColumn[iColumn] = put - startColumn;
    00918     //move largest
    00919     if ( positionLargest >= 0 ) {
    00920       value = elementU[positionLargest];
    00921       iRow = indexRowU[positionLargest];
    00922       elementU[positionLargest] = elementU[startColumn];
    00923       indexRowU[positionLargest] = indexRowU[startColumn];
    00924       elementU[startColumn] = value;
    00925       indexRowU[startColumn] = iRow;
    00926     }
    00927     //linked list for column
    00928     if ( nextCount[iColumn + numberRows_] != -2 ) {
    00929       //modify linked list
    00930       deleteLink ( iColumn + numberRows_ );
    00931       addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
    00932     }
    00933     temp2 += increment2;
    00934   }
    00935   //get space for row list
    00936   unsigned int *putBase = workArea2;
    00937   int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
    00938   int i = 0;
    00939
    00940   // do linked lists and update counts
    00941   while ( bigLoops ) {
    00942     bigLoops--;
    00943     int bit;
    00944     for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
    00945       unsigned int *putThis = putBase;
    00946       int iRow = indexL[i];
    00947
    00948       //get space
    00949       int number = 0;
    00950       int jColumn;
    00951
    00952       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    00953         unsigned int test = *putThis;
    00954
    00955         putThis += increment2;
    00956         test = 1 - ( ( test >> bit ) & 1 );
    00957         number += test;
    00958       }
    00959       int next = nextRow[iRow];
    00960       CoinBigIndex space;
    00961
    00962       space = startRowU[next] - startRowU[iRow];
    00963       number += numberInRow[iRow];
    00964       if ( space < number ) {
    00965         if ( !getRowSpace ( iRow, number ) ) {
    00966           return false;
    00967         }
    00968       }
    00969       // now do
    00970       putThis = putBase;
    00971       next = nextRow[iRow];
    00972       number = numberInRow[iRow];
    00973       CoinBigIndex end = startRowU[iRow] + number;
    00974       int saveIndex = indexColumnU[startRowU[next]];
    00975
    00976       //add in
    00977       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    00978         unsigned int test = *putThis;
    00979
    00980         putThis += increment2;
    00981         test = 1 - ( ( test >> bit ) & 1 );
    00982         indexColumnU[end] = saveColumn[jColumn];
    00983         end += test;
    00984       }
    00985       //put back next one in case zapped
    00986       indexColumnU[startRowU[next]] = saveIndex;
    00987       markRow[iRow] = -1;
    00988       number = end - startRowU[iRow];
    00989       numberInRow[iRow] = number;
    00990       deleteLink ( iRow );
    00991       addLink ( iRow, number );
    00992     }
    00993     putBase++;
    00994   }                             /* endwhile */
    00995   int bit;
    00996
    00997   for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
    00998     unsigned int *putThis = putBase;
    00999     int iRow = indexL[i];
    01000
    01001     //get space
    01002     int number = 0;
    01003     int jColumn;
    01004
    01005     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01006       unsigned int test = *putThis;
    01007
    01008       putThis += increment2;
    01009       test = 1 - ( ( test >> bit ) & 1 );
    01010       number += test;
    01011     }
    01012     int next = nextRow[iRow];
    01013     CoinBigIndex space;
    01014
    01015     space = startRowU[next] - startRowU[iRow];
    01016     number += numberInRow[iRow];
    01017     if ( space < number ) {
    01018       if ( !getRowSpace ( iRow, number ) ) {
    01019         return false;
    01020       }
    01021     }
    01022     // now do
    01023     putThis = putBase;
    01024     next = nextRow[iRow];
    01025     number = numberInRow[iRow];
    01026     CoinBigIndex end = startRowU[iRow] + number;
    01027     int saveIndex;
    01028
    01029     saveIndex = indexColumnU[startRowU[next]];
    01030
    01031     //add in
    01032     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01033       unsigned int test = *putThis;
    01034
    01035       putThis += increment2;
    01036       test = 1 - ( ( test >> bit ) & 1 );
    01037
    01038       indexColumnU[end] = saveColumn[jColumn];
    01039       end += test;
    01040     }
    01041     indexColumnU[startRowU[next]] = saveIndex;
    01042     markRow[iRow] = -1;
    01043     number = end - startRowU[iRow];
    01044     numberInRow[iRow] = number;
    01045     deleteLink ( iRow );
    01046     addLink ( iRow, number );
    01047   }
    01048   markRow[pivotRow] = -1;
    01049   //modify linked list for pivots
    01050   deleteLink ( pivotRow );
    01051   deleteLink ( pivotColumn + numberRows_ );
    01052   totalElements_ += added;
    01053   return true;
    01054 }
    01055
    01056   /********************************* END LARGE TEMPLATE ********/
    01058
    01059 protected:
    01060
    01063
    01064   double pivotTolerance_;
    01066   double zeroTolerance_;
    01068   double slackValue_;
    01070   double areaFactor_;
    01072   double relaxCheck_;
    01074   int numberRows_;
    01076   int numberRowsExtra_;
    01078   int maximumRowsExtra_;
    01080   int numberColumns_;
    01082   int numberColumnsExtra_;
    01084   int maximumColumnsExtra_;
    01086   int numberGoodU_;
    01088   int numberGoodL_;
    01090   int maximumPivots_;
    01092   int numberPivots_;
    01095   CoinBigIndex totalElements_;
    01097   CoinBigIndex factorElements_;
    01099   CoinIntArrayWithLength pivotColumn_;
    01101   CoinIntArrayWithLength permute_;
    01103   CoinIntArrayWithLength permuteBack_;
    01105   CoinIntArrayWithLength pivotColumnBack_;
    01107   int status_;
    01108
    01113   //int increasingRows_;
    01114
    01116   int numberTrials_;
    01118   CoinBigIndexArrayWithLength startRowU_;
    01119
    01121   CoinIntArrayWithLength numberInRow_;
    01122
    01124   CoinIntArrayWithLength numberInColumn_;
    01125
    01127   CoinIntArrayWithLength numberInColumnPlus_;
    01128
    01131   CoinIntArrayWithLength firstCount_;
    01132
    01134   CoinIntArrayWithLength nextCount_;
    01135
    01137   CoinIntArrayWithLength lastCount_;
    01138
    01140   CoinIntArrayWithLength nextColumn_;
    01141
    01143   CoinIntArrayWithLength lastColumn_;
    01144
    01146   CoinIntArrayWithLength nextRow_;
    01147
    01149   CoinIntArrayWithLength lastRow_;
    01150
    01152   CoinIntArrayWithLength saveColumn_;
    01153
    01155   CoinIntArrayWithLength markRow_;
    01156
    01158   int messageLevel_;
    01159
    01161   int biggerDimension_;
    01162
    01164   CoinIntArrayWithLength indexColumnU_;
    01165
    01167   CoinIntArrayWithLength pivotRowL_;
    01168
    01170   CoinDoubleArrayWithLength pivotRegion_;
    01171
    01173   int numberSlacks_;
    01174
    01176   int numberU_;
    01177
    01179   CoinBigIndex maximumU_;
    01180
    01182   //int baseU_;
    01183
    01185   CoinBigIndex lengthU_;
    01186
    01188   CoinBigIndex lengthAreaU_;
    01189
    01191   CoinDoubleArrayWithLength elementU_;
    01192
    01194   CoinIntArrayWithLength indexRowU_;
    01195
    01197   CoinBigIndexArrayWithLength startColumnU_;
    01198
    01200   CoinBigIndexArrayWithLength convertRowToColumnU_;
    01201
    01203   int numberL_;
    01204
    01206   int baseL_;
    01207
    01209   CoinBigIndex lengthL_;
    01210
    01212   CoinBigIndex lengthAreaL_;
    01213
    01215   CoinDoubleArrayWithLength elementL_;
    01216
    01218   CoinIntArrayWithLength indexRowL_;
    01219
    01221   CoinBigIndexArrayWithLength startColumnL_;
    01222
    01224   bool doForrestTomlin_;
    01225
    01227   int numberR_;
    01228
    01230   CoinBigIndex lengthR_;
    01231
    01233   CoinBigIndex lengthAreaR_;
    01234
    01236   double *elementR_;
    01237
    01239   int *indexRowR_;
    01240
    01242   CoinBigIndexArrayWithLength startColumnR_;
    01243
    01245   double  * denseArea_;
    01246
    01248   int * densePermute_;
    01249
    01251   int numberDense_;
    01252
    01254   int denseThreshold_;
    01255
    01257   CoinDoubleArrayWithLength workArea_;
    01258
    01260   CoinUnsignedIntArrayWithLength workArea2_;
    01261
    01263   CoinBigIndex numberCompressions_;
    01264
    01266   mutable double ftranCountInput_;
    01267   mutable double ftranCountAfterL_;
    01268   mutable double ftranCountAfterR_;
    01269   mutable double ftranCountAfterU_;
    01270   mutable double btranCountInput_;
    01271   mutable double btranCountAfterU_;
    01272   mutable double btranCountAfterR_;
    01273   mutable double btranCountAfterL_;
    01274
    01276   mutable int numberFtranCounts_;
    01277   mutable int numberBtranCounts_;
    01278
    01280   double ftranAverageAfterL_;
    01281   double ftranAverageAfterR_;
    01282   double ftranAverageAfterU_;
    01283   double btranAverageAfterU_;
    01284   double btranAverageAfterR_;
    01285   double btranAverageAfterL_;
    01286
    01288   mutable bool collectStatistics_;
    01289
    01291   int sparseThreshold_;
    01292
    01294   int sparseThreshold2_;
    01295
    01297   CoinBigIndexArrayWithLength startRowL_;
    01298
    01300   CoinIntArrayWithLength indexColumnL_;
    01301
    01303   CoinDoubleArrayWithLength elementByRowL_;
    01304
    01306   mutable CoinIntArrayWithLength sparse_;
    01310   int biasLU_;
    01316   int persistenceFlag_;
    01318 };
    01319 // Dense coding
    01320 #ifdef COIN_HAS_LAPACK
    01321 #define DENSE_CODE 1
    01322 /* Type of Fortran integer translated into C */
    01323 #ifndef ipfint
    01324 //typedef ipfint FORTRAN_INTEGER_TYPE ;
    01325 typedef int ipfint;
    01326 typedef const int cipfint;
    01327 #endif
    01328 #endif
    01329 #endif