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

- File List
- File Members

# /home/coin/SVN-release/Cgl-0.52.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   double conditionNumber() const;
    00126
    00128
    00131
    00132   inline int status (  ) const {
    00133     return status_;
    00134   }
    00136   inline void setStatus (  int value)
    00137   {  status_=value;  }
    00139   inline int pivots (  ) const {
    00140     return numberPivots_;
    00141   }
    00143   inline void setPivots (  int value )
    00144   { numberPivots_=value; }
    00146   inline int *permute (  ) const {
    00147     return permute_.array();
    00148   }
    00150   inline int *pivotColumn (  ) const {
    00151     return pivotColumn_.array();
    00152   }
    00154   inline double *pivotRegion (  ) const {
    00155     return pivotRegion_.array();
    00156   }
    00158   inline int *permuteBack (  ) const {
    00159     return permuteBack_.array();
    00160   }
    00162   inline int *pivotColumnBack (  ) const {
    00163     return pivotColumnBack_.array();
    00164   }
    00166   inline CoinBigIndex * startRowL() const
    00167   { return startRowL_.array();}
    00168
    00170   inline CoinBigIndex * startColumnL() const
    00171   { return startColumnL_.array();}
    00172
    00174   inline int * indexColumnL() const
    00175   { return indexColumnL_.array();}
    00176
    00178   inline int * indexRowL() const
    00179   { return indexRowL_.array();}
    00180
    00182   inline double * elementByRowL() const
    00183   { return elementByRowL_.array();}
    00184
    00186   inline int numberRowsExtra (  ) const {
    00187     return numberRowsExtra_;
    00188   }
    00190   inline void setNumberRows(int value)
    00191   { numberRows_ = value; }
    00193   inline int numberRows (  ) const {
    00194     return numberRows_;
    00195   }
    00197   inline CoinBigIndex numberL() const
    00198   { return numberL_;}
    00199
    00201   inline CoinBigIndex baseL() const
    00202   { return baseL_;}
    00204   inline int maximumRowsExtra (  ) const {
    00205     return maximumRowsExtra_;
    00206   }
    00208   inline int numberColumns (  ) const {
    00209     return numberColumns_;
    00210   }
    00212   inline int numberElements (  ) const {
    00213     return totalElements_;
    00214   }
    00216   inline int numberForrestTomlin (  ) const {
    00217     return numberInColumn_.array()[numberColumnsExtra_];
    00218   }
    00220   inline int numberGoodColumns (  ) const {
    00221     return numberGoodU_;
    00222   }
    00224   inline double areaFactor (  ) const {
    00225     return areaFactor_;
    00226   }
    00227   inline void areaFactor ( double value ) {
    00228     areaFactor_=value;
    00229   }
    00231   double adjustedAreaFactor() const;
    00233   inline void relaxAccuracyCheck(double value)
    00234   { relaxCheck_ = value;}
    00235   inline double getAccuracyCheck() const
    00236   { return relaxCheck_;}
    00238   inline bool increasingRows (  ) const
    00239   { return true; }
    00245   inline void increasingRows ( int value  ) {}
    00247   inline int messageLevel (  ) const {
    00248     return messageLevel_ ;
    00249   }
    00250   void messageLevel (  int value );
    00252   inline int maximumPivots (  ) const {
    00253     return maximumPivots_ ;
    00254   }
    00255   void maximumPivots (  int value );
    00256
    00258   inline int denseThreshold() const
    00259     { return denseThreshold_;}
    00261   inline void setDenseThreshold(int value)
    00262     { denseThreshold_ = value;}
    00264   inline double pivotTolerance (  ) const {
    00265     return pivotTolerance_ ;
    00266   }
    00267   void pivotTolerance (  double value );
    00269   inline double zeroTolerance (  ) const {
    00270     return zeroTolerance_ ;
    00271   }
    00272   void zeroTolerance (  double value );
    00274   inline double slackValue (  ) const {
    00275     return slackValue_ ;
    00276   }
    00277   void slackValue (  double value );
    00279   double maximumCoefficient() const;
    00281   inline bool forrestTomlin() const
    00282   { return doForrestTomlin_;}
    00283   inline void setForrestTomlin(bool value)
    00284   { doForrestTomlin_=value;}
    00286
    00289
    00291   inline int numberDense() const
    00292   { return numberDense_;}
    00293
    00295   inline CoinBigIndex numberElementsU (  ) const {
    00296     return lengthU_;
    00297   }
    00299   inline void setNumberElementsU(CoinBigIndex value)
    00300   { lengthU_ = value; }
    00302   inline CoinBigIndex lengthAreaU (  ) const {
    00303     return lengthAreaU_;
    00304   }
    00306   inline CoinBigIndex numberElementsL (  ) const {
    00307     return lengthL_;
    00308   }
    00310   inline CoinBigIndex lengthAreaL (  ) const {
    00311     return lengthAreaL_;
    00312   }
    00314   inline CoinBigIndex numberElementsR (  ) const {
    00315     return lengthR_;
    00316   }
    00318   inline CoinBigIndex numberCompressions() const
    00319   { return numberCompressions_;}
    00321   inline int * numberInRow() const
    00322   { return numberInRow_.array();}
    00324   inline int * numberInColumn() const
    00325   { return numberInColumn_.array();}
    00327   inline double * elementU() const
    00328   { return elementU_.array();}
    00330   inline int * indexRowU() const
    00331   { return indexRowU_.array();}
    00333   inline CoinBigIndex * startColumnU() const
    00334   { return startColumnU_.array();}
    00336   inline int maximumColumnsExtra()
    00337   { return maximumColumnsExtra_;}
    00341   inline int biasLU() const
    00342   { return biasLU_;}
    00343   inline void setBiasLU(int value)
    00344   { biasLU_=value;}
    00350   inline int persistenceFlag() const
    00351   { return persistenceFlag_;}
    00352   void setPersistenceFlag(int value);
    00354
    00357
    00365   int replaceColumn ( CoinIndexedVector * regionSparse,
    00366                       int pivotRow,
    00367                       double pivotCheck ,
    00368                       bool checkBeforeModifying=false);
    00370
    00380   int updateColumnFT ( CoinIndexedVector * regionSparse,
    00381                        CoinIndexedVector * regionSparse2);
    00384   int updateColumn ( CoinIndexedVector * regionSparse,
    00385                      CoinIndexedVector * regionSparse2,
    00386                      bool noPermute=false) const;
    00391   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
    00392                               CoinIndexedVector * regionSparse2) const;
    00394   void goSparse();
    00396   inline int sparseThreshold ( ) const
    00397   { return sparseThreshold_;}
    00399   void sparseThreshold ( int value );
    00401
    00402
    00406
    00407   inline void clearArrays()
    00408   { gutsOfDestructor();}
    00410
    00413
    00416   int add ( CoinBigIndex numberElements,
    00417                int indicesRow[],
    00418                int indicesColumn[], double elements[] );
    00419
    00422   int addColumn ( CoinBigIndex numberElements,
    00423                      int indicesRow[], double elements[] );
    00424
    00427   int addRow ( CoinBigIndex numberElements,
    00428                   int indicesColumn[], double elements[] );
    00429
    00431   int deleteColumn ( int Row );
    00433   int deleteRow ( int Row );
    00434
    00438   int replaceRow ( int whichRow, int numberElements,
    00439                       const int indicesColumn[], const double elements[] );
    00441   void emptyRows(int numberToEmpty, const int which[]);
    00443
    00444
    00445   void checkSparse();
    00447   inline bool collectStatistics() const
    00448   { return collectStatistics_;}
    00450   inline void setCollectStatistics(bool onOff) const
    00451   { collectStatistics_ = onOff;}
    00453   void gutsOfDestructor(int type=1);
    00455   void gutsOfInitialize(int type);
    00456   void gutsOfCopy(const CoinFactorization &other);
    00457
    00459   void resetStatistics();
    00460
    00461
    00463
    00465
    00466   void getAreas ( int numberRows,
    00467                   int numberColumns,
    00468                   CoinBigIndex maximumL,
    00469                   CoinBigIndex maximumU );
    00470
    00473   void preProcess ( int state,
    00474                     int possibleDuplicates = -1 );
    00476   int factor (  );
    00477 protected:
    00480   int factorSparse (  );
    00483   int factorDense (  );
    00484
    00486   bool pivotOneOtherRow ( int pivotRow,
    00487                           int pivotColumn );
    00489   bool pivotRowSingleton ( int pivotRow,
    00490                            int pivotColumn );
    00492   bool pivotColumnSingleton ( int pivotRow,
    00493                               int pivotColumn );
    00494
    00499   bool getColumnSpace ( int iColumn,
    00500                         int extraNeeded );
    00501
    00505   bool getColumnSpaceIterateR ( int iColumn, double value,
    00506                                int iRow);
    00512   CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
    00513                                int iRow);
    00517   bool getRowSpace ( int iRow, int extraNeeded );
    00518
    00522   bool getRowSpaceIterate ( int iRow,
    00523                             int extraNeeded );
    00525   void checkConsistency (  );
    00527   inline void addLink ( int index, int count ) {
    00528     int *nextCount = nextCount_.array();
    00529     int *firstCount = firstCount_.array();
    00530     int *lastCount = lastCount_.array();
    00531     int next = firstCount[count];
    00532       lastCount[index] = -2 - count;
    00533     if ( next < 0 ) {
    00534       //first with that count
    00535       firstCount[count] = index;
    00536       nextCount[index] = -1;
    00537     } else {
    00538       firstCount[count] = index;
    00539       nextCount[index] = next;
    00540       lastCount[next] = index;
    00541   }}
    00543   inline void deleteLink ( int index ) {
    00544     int *nextCount = nextCount_.array();
    00545     int *firstCount = firstCount_.array();
    00546     int *lastCount = lastCount_.array();
    00547     int next = nextCount[index];
    00548     int last = lastCount[index];
    00549     if ( last >= 0 ) {
    00550       nextCount[last] = next;
    00551     } else {
    00552       int count = -last - 2;
    00553
    00554       firstCount[count] = next;
    00555     }
    00556     if ( next >= 0 ) {
    00557       lastCount[next] = last;
    00558     }
    00559     nextCount[index] = -2;
    00560     lastCount[index] = -2;
    00561     return;
    00562   }
    00564   void separateLinks(int count,bool rowsFirst);
    00566   void cleanup (  );
    00567
    00569   void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
    00571   void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
    00573   void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
    00575   void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
    00576
    00578   void updateColumnR ( CoinIndexedVector * region ) const;
    00581   void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
    00582
    00584   void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
    00585
    00587   void updateColumnUSparse ( CoinIndexedVector * regionSparse,
    00588                              int * indexIn) const;
    00590   void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
    00591                                int * indexIn) const;
    00593   void updateColumnUDensish ( CoinIndexedVector * regionSparse,
    00594                               int * indexIn) const;
    00596   void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
    00598   void permuteBack ( CoinIndexedVector * regionSparse,
    00599                      CoinIndexedVector * outVector) const;
    00600
    00602   void updateColumnTransposePFI ( CoinIndexedVector * region) const;
    00605   void updateColumnTransposeU ( CoinIndexedVector * region,
    00606                                 int smallestIndex) const;
    00609   void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
    00610                                         int smallestIndex) const;
    00613   void updateColumnTransposeUDensish ( CoinIndexedVector * region,
    00614                                        int smallestIndex) const;
    00617   void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
    00618
    00620   void updateColumnTransposeR ( CoinIndexedVector * region ) const;
    00622   void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
    00624   void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
    00625
    00627   void updateColumnTransposeL ( CoinIndexedVector * region ) const;
    00629   void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
    00631   void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
    00633   void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
    00635   void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
    00636 public:
    00641   int replaceColumnPFI ( CoinIndexedVector * regionSparse,
    00642                          int pivotRow, double alpha);
    00643 protected:
    00646   int checkPivot(double saveFromU, double oldPivot) const;
    00647   /********************************* START LARGE TEMPLATE ********/
    00648 #ifdef INT_IS_8
    00649 #define COINFACTORIZATION_BITS_PER_INT 64
    00650 #define COINFACTORIZATION_SHIFT_PER_INT 6
    00651 #define COINFACTORIZATION_MASK_PER_INT 0x3f
    00652 #else
    00653 #define COINFACTORIZATION_BITS_PER_INT 32
    00654 #define COINFACTORIZATION_SHIFT_PER_INT 5
    00655 #define COINFACTORIZATION_MASK_PER_INT 0x1f
    00656 #endif
    00657   template <class T>  inline bool
    00658   pivot ( int pivotRow,
    00659           int pivotColumn,
    00660           CoinBigIndex pivotRowPosition,
    00661           CoinBigIndex pivotColumnPosition,
    00662           double work[],
    00663           unsigned int workArea2[],
    00664           int increment,
    00665           int increment2,
    00666           T markRow[] ,
    00667           int largeInteger)
    00668 {
    00669   int *indexColumnU = indexColumnU_.array();
    00670   CoinBigIndex *startColumnU = startColumnU_.array();
    00671   int *numberInColumn = numberInColumn_.array();
    00672   double *elementU = elementU_.array();
    00673   int *indexRowU = indexRowU_.array();
    00674   CoinBigIndex *startRowU = startRowU_.array();
    00675   int *numberInRow = numberInRow_.array();
    00676   double *elementL = elementL_.array();
    00677   int *indexRowL = indexRowL_.array();
    00678   int *saveColumn = saveColumn_.array();
    00679   int *nextRow = nextRow_.array();
    00680   int *lastRow = lastRow_.array() ;
    00681
    00682   //store pivot columns (so can easily compress)
    00683   int numberInPivotRow = numberInRow[pivotRow] - 1;
    00684   CoinBigIndex startColumn = startColumnU[pivotColumn];
    00685   int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
    00686   CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
    00687   int put = 0;
    00688   CoinBigIndex startRow = startRowU[pivotRow];
    00689   CoinBigIndex endRow = startRow + numberInPivotRow + 1;
    00690
    00691   if ( pivotColumnPosition < 0 ) {
    00692     for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    00693       int iColumn = indexColumnU[pivotColumnPosition];
    00694       if ( iColumn != pivotColumn ) {
    00695         saveColumn[put++] = iColumn;
    00696       } else {
    00697         break;
    00698       }
    00699     }
    00700   } else {
    00701     for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
    00702       saveColumn[put++] = indexColumnU[i];
    00703     }
    00704   }
    00705   assert (pivotColumnPosition<endRow);
    00706   assert (indexColumnU[pivotColumnPosition]==pivotColumn);
    00707   pivotColumnPosition++;
    00708   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    00709     saveColumn[put++] = indexColumnU[pivotColumnPosition];
    00710   }
    00711   //take out this bit of indexColumnU
    00712   int next = nextRow[pivotRow];
    00713   int last = lastRow[pivotRow];
    00714
    00715   nextRow[last] = next;
    00716   lastRow[next] = last;
    00717   nextRow[pivotRow] = numberGoodU_;     //use for permute
    00718   lastRow[pivotRow] = -2;
    00719   numberInRow[pivotRow] = 0;
    00720   //store column in L, compress in U and take column out
    00721   CoinBigIndex l = lengthL_;
    00722
    00723   if ( l + numberInPivotColumn > lengthAreaL_ ) {
    00724     //need more memory
    00725     printf("more memory needed in middle of invert\n");
    00726     return false;
    00727   }
    00728   //l+=currentAreaL_->elementByColumn-elementL;
    00729   CoinBigIndex lSave = l;
    00730
    00731   pivotRowL_.array()[numberGoodL_] = pivotRow;
    00732   CoinBigIndex * startColumnL = startColumnL_.array();
    00733   startColumnL[numberGoodL_] = l;       //for luck and first time
    00734   numberGoodL_++;
    00735   startColumnL[numberGoodL_] = l + numberInPivotColumn;
    00736   lengthL_ += numberInPivotColumn;
    00737   if ( pivotRowPosition < 0 ) {
    00738     for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    00739       int iRow = indexRowU[pivotRowPosition];
    00740       if ( iRow != pivotRow ) {
    00741         indexRowL[l] = iRow;
    00742         elementL[l] = elementU[pivotRowPosition];
    00743         markRow[iRow] = l - lSave;
    00744         l++;
    00745         //take out of row list
    00746         CoinBigIndex start = startRowU[iRow];
    00747         CoinBigIndex end = start + numberInRow[iRow];
    00748         CoinBigIndex where = start;
    00749
    00750         while ( indexColumnU[where] != pivotColumn ) {
    00751           where++;
    00752         }                       /* endwhile */
    00753 #if DEBUG_COIN
    00754         if ( where >= end ) {
    00755           abort (  );
    00756         }
    00757 #endif
    00758         indexColumnU[where] = indexColumnU[end - 1];
    00759         numberInRow[iRow]--;
    00760       } else {
    00761         break;
    00762       }
    00763     }
    00764   } else {
    00765     CoinBigIndex i;
    00766
    00767     for ( i = startColumn; i < pivotRowPosition; i++ ) {
    00768       int iRow = indexRowU[i];
    00769
    00770       markRow[iRow] = l - lSave;
    00771       indexRowL[l] = iRow;
    00772       elementL[l] = elementU[i];
    00773       l++;
    00774       //take out of row list
    00775       CoinBigIndex start = startRowU[iRow];
    00776       CoinBigIndex end = start + numberInRow[iRow];
    00777       CoinBigIndex where = start;
    00778
    00779       while ( indexColumnU[where] != pivotColumn ) {
    00780         where++;
    00781       }                         /* endwhile */
    00782 #if DEBUG_COIN
    00783       if ( where >= end ) {
    00784         abort (  );
    00785       }
    00786 #endif
    00787       indexColumnU[where] = indexColumnU[end - 1];
    00788       numberInRow[iRow]--;
    00789       assert (numberInRow[iRow]>=0);
    00790     }
    00791   }
    00792   assert (pivotRowPosition<endColumn);
    00793   assert (indexRowU[pivotRowPosition]==pivotRow);
    00794   double pivotElement = elementU[pivotRowPosition];
    00795   double pivotMultiplier = 1.0 / pivotElement;
    00796
    00797   pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
    00798   pivotRowPosition++;
    00799   for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    00800     int iRow = indexRowU[pivotRowPosition];
    00801
    00802     markRow[iRow] = l - lSave;
    00803     indexRowL[l] = iRow;
    00804     elementL[l] = elementU[pivotRowPosition];
    00805     l++;
    00806     //take out of row list
    00807     CoinBigIndex start = startRowU[iRow];
    00808     CoinBigIndex end = start + numberInRow[iRow];
    00809     CoinBigIndex where = start;
    00810
    00811     while ( indexColumnU[where] != pivotColumn ) {
    00812       where++;
    00813     }                           /* endwhile */
    00814 #if DEBUG_COIN
    00815     if ( where >= end ) {
    00816       abort (  );
    00817     }
    00818 #endif
    00819     indexColumnU[where] = indexColumnU[end - 1];
    00820     numberInRow[iRow]--;
    00821     assert (numberInRow[iRow]>=0);
    00822   }
    00823   markRow[pivotRow] = largeInteger;
    00824   //compress pivot column (move pivot to front including saved)
    00825   numberInColumn[pivotColumn] = 0;
    00826   //use end of L for temporary space
    00827   int *indexL = &indexRowL[lSave];
    00828   double *multipliersL = &elementL[lSave];
    00829
    00830   //adjust
    00831   int j;
    00832
    00833   for ( j = 0; j < numberInPivotColumn; j++ ) {
    00834     multipliersL[j] *= pivotMultiplier;
    00835   }
    00836   //zero out fill
    00837   CoinBigIndex iErase;
    00838   for ( iErase = 0; iErase < increment2 * numberInPivotRow;
    00839         iErase++ ) {
    00840     workArea2[iErase] = 0;
    00841   }
    00842   CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
    00843   unsigned int *temp2 = workArea2;
    00844   int * nextColumn = nextColumn_.array();
    00845
    00846   //pack down and move to work
    00847   int jColumn;
    00848   for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    00849     int iColumn = saveColumn[jColumn];
    00850     CoinBigIndex startColumn = startColumnU[iColumn];
    00851     CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
    00852     int iRow = indexRowU[startColumn];
    00853     double value = elementU[startColumn];
    00854     double largest;
    00855     CoinBigIndex put = startColumn;
    00856     CoinBigIndex positionLargest = -1;
    00857     double thisPivotValue = 0.0;
    00858
    00859     //compress column and find largest not updated
    00860     bool checkLargest;
    00861     int mark = markRow[iRow];
    00862
    00863     if ( mark < 0 ) {
    00864       largest = fabs ( value );
    00865       positionLargest = put;
    00866       put++;
    00867       checkLargest = false;
    00868     } else {
    00869       //need to find largest
    00870       largest = 0.0;
    00871       checkLargest = true;
    00872       if ( mark != largeInteger ) {
    00873         //will be updated
    00874         work[mark] = value;
    00875         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    00876         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    00877
    00878         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
    00879         added--;
    00880       } else {
    00881         thisPivotValue = value;
    00882       }
    00883     }
    00884     CoinBigIndex i;
    00885     for ( i = startColumn + 1; i < endColumn; i++ ) {
    00886       iRow = indexRowU[i];
    00887       value = elementU[i];
    00888       int mark = markRow[iRow];
    00889
    00890       if ( mark < 0 ) {
    00891         //keep
    00892         indexRowU[put] = iRow;
    00893         elementU[put] = value;
    00894         if ( checkLargest ) {
    00895           double absValue = fabs ( value );
    00896
    00897           if ( absValue > largest ) {
    00898             largest = absValue;
    00899             positionLargest = put;
    00900           }
    00901         }
    00902         put++;
    00903       } else if ( mark != largeInteger ) {
    00904         //will be updated
    00905         work[mark] = value;
    00906         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    00907         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    00908
    00909         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
    00910         added--;
    00911       } else {
    00912         thisPivotValue = value;
    00913       }
    00914     }
    00915     //slot in pivot
    00916     elementU[put] = elementU[startColumn];
    00917     indexRowU[put] = indexRowU[startColumn];
    00918     if ( positionLargest == startColumn ) {
    00919       positionLargest = put;    //follow if was largest
    00920     }
    00921     put++;
    00922     elementU[startColumn] = thisPivotValue;
    00923     indexRowU[startColumn] = pivotRow;
    00924     //clean up counts
    00925     startColumn++;
    00926     numberInColumn[iColumn] = put - startColumn;
    00927     int * numberInColumnPlus = numberInColumnPlus_.array();
    00928     numberInColumnPlus[iColumn]++;
    00929     startColumnU[iColumn]++;
    00930     //how much space have we got
    00931     int next = nextColumn[iColumn];
    00932     CoinBigIndex space;
    00933
    00934     space = startColumnU[next] - put - numberInColumnPlus[next];
    00935     //assume no zero elements
    00936     if ( numberInPivotColumn > space ) {
    00937       //getColumnSpace also moves fixed part
    00938       if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
    00939         return false;
    00940       }
    00941       //redo starts
    00942       positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
    00943       startColumn = startColumnU[iColumn];
    00944       put = startColumn + numberInColumn[iColumn];
    00945     }
    00946     double tolerance = zeroTolerance_;
    00947
    00948     int *nextCount = nextCount_.array();
    00949     for ( j = 0; j < numberInPivotColumn; j++ ) {
    00950       value = work[j] - thisPivotValue * multipliersL[j];
    00951       double absValue = fabs ( value );
    00952
    00953       if ( absValue > tolerance ) {
    00954         work[j] = 0.0;
    00955         elementU[put] = value;
    00956         indexRowU[put] = indexL[j];
    00957         if ( absValue > largest ) {
    00958           largest = absValue;
    00959           positionLargest = put;
    00960         }
    00961         put++;
    00962       } else {
    00963         work[j] = 0.0;
    00964         added--;
    00965         int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    00966         int bit = j & COINFACTORIZATION_MASK_PER_INT;
    00967
    00968         if ( temp2[word] & ( 1 << bit ) ) {
    00969           //take out of row list
    00970           iRow = indexL[j];
    00971           CoinBigIndex start = startRowU[iRow];
    00972           CoinBigIndex end = start + numberInRow[iRow];
    00973           CoinBigIndex where = start;
    00974
    00975           while ( indexColumnU[where] != iColumn ) {
    00976             where++;
    00977           }                     /* endwhile */
    00978 #if DEBUG_COIN
    00979           if ( where >= end ) {
    00980             abort (  );
    00981           }
    00982 #endif
    00983           indexColumnU[where] = indexColumnU[end - 1];
    00984           numberInRow[iRow]--;
    00985         } else {
    00986           //make sure won't be added
    00987           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    00988           int bit = j & COINFACTORIZATION_MASK_PER_INT;
    00989
    00990           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    00991         }
    00992       }
    00993     }
    00994     numberInColumn[iColumn] = put - startColumn;
    00995     //move largest
    00996     if ( positionLargest >= 0 ) {
    00997       value = elementU[positionLargest];
    00998       iRow = indexRowU[positionLargest];
    00999       elementU[positionLargest] = elementU[startColumn];
    01000       indexRowU[positionLargest] = indexRowU[startColumn];
    01001       elementU[startColumn] = value;
    01002       indexRowU[startColumn] = iRow;
    01003     }
    01004     //linked list for column
    01005     if ( nextCount[iColumn + numberRows_] != -2 ) {
    01006       //modify linked list
    01007       deleteLink ( iColumn + numberRows_ );
    01008       addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
    01009     }
    01010     temp2 += increment2;
    01011   }
    01012   //get space for row list
    01013   unsigned int *putBase = workArea2;
    01014   int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
    01015   int i = 0;
    01016
    01017   // do linked lists and update counts
    01018   while ( bigLoops ) {
    01019     bigLoops--;
    01020     int bit;
    01021     for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
    01022       unsigned int *putThis = putBase;
    01023       int iRow = indexL[i];
    01024
    01025       //get space
    01026       int number = 0;
    01027       int jColumn;
    01028
    01029       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01030         unsigned int test = *putThis;
    01031
    01032         putThis += increment2;
    01033         test = 1 - ( ( test >> bit ) & 1 );
    01034         number += test;
    01035       }
    01036       int next = nextRow[iRow];
    01037       CoinBigIndex space;
    01038
    01039       space = startRowU[next] - startRowU[iRow];
    01040       number += numberInRow[iRow];
    01041       if ( space < number ) {
    01042         if ( !getRowSpace ( iRow, number ) ) {
    01043           return false;
    01044         }
    01045       }
    01046       // now do
    01047       putThis = putBase;
    01048       next = nextRow[iRow];
    01049       number = numberInRow[iRow];
    01050       CoinBigIndex end = startRowU[iRow] + number;
    01051       int saveIndex = indexColumnU[startRowU[next]];
    01052
    01053       //add in
    01054       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01055         unsigned int test = *putThis;
    01056
    01057         putThis += increment2;
    01058         test = 1 - ( ( test >> bit ) & 1 );
    01059         indexColumnU[end] = saveColumn[jColumn];
    01060         end += test;
    01061       }
    01062       //put back next one in case zapped
    01063       indexColumnU[startRowU[next]] = saveIndex;
    01064       markRow[iRow] = -1;
    01065       number = end - startRowU[iRow];
    01066       numberInRow[iRow] = number;
    01067       deleteLink ( iRow );
    01068       addLink ( iRow, number );
    01069     }
    01070     putBase++;
    01071   }                             /* endwhile */
    01072   int bit;
    01073
    01074   for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
    01075     unsigned int *putThis = putBase;
    01076     int iRow = indexL[i];
    01077
    01078     //get space
    01079     int number = 0;
    01080     int jColumn;
    01081
    01082     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01083       unsigned int test = *putThis;
    01084
    01085       putThis += increment2;
    01086       test = 1 - ( ( test >> bit ) & 1 );
    01087       number += test;
    01088     }
    01089     int next = nextRow[iRow];
    01090     CoinBigIndex space;
    01091
    01092     space = startRowU[next] - startRowU[iRow];
    01093     number += numberInRow[iRow];
    01094     if ( space < number ) {
    01095       if ( !getRowSpace ( iRow, number ) ) {
    01096         return false;
    01097       }
    01098     }
    01099     // now do
    01100     putThis = putBase;
    01101     next = nextRow[iRow];
    01102     number = numberInRow[iRow];
    01103     CoinBigIndex end = startRowU[iRow] + number;
    01104     int saveIndex;
    01105
    01106     saveIndex = indexColumnU[startRowU[next]];
    01107
    01108     //add in
    01109     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01110       unsigned int test = *putThis;
    01111
    01112       putThis += increment2;
    01113       test = 1 - ( ( test >> bit ) & 1 );
    01114
    01115       indexColumnU[end] = saveColumn[jColumn];
    01116       end += test;
    01117     }
    01118     indexColumnU[startRowU[next]] = saveIndex;
    01119     markRow[iRow] = -1;
    01120     number = end - startRowU[iRow];
    01121     numberInRow[iRow] = number;
    01122     deleteLink ( iRow );
    01123     addLink ( iRow, number );
    01124   }
    01125   markRow[pivotRow] = -1;
    01126   //modify linked list for pivots
    01127   deleteLink ( pivotRow );
    01128   deleteLink ( pivotColumn + numberRows_ );
    01129   totalElements_ += added;
    01130   return true;
    01131 }
    01132
    01133   /********************************* END LARGE TEMPLATE ********/
    01135
    01136 protected:
    01137
    01140
    01141   double pivotTolerance_;
    01143   double zeroTolerance_;
    01145   double slackValue_;
    01147   double areaFactor_;
    01149   double relaxCheck_;
    01151   int numberRows_;
    01153   int numberRowsExtra_;
    01155   int maximumRowsExtra_;
    01157   int numberColumns_;
    01159   int numberColumnsExtra_;
    01161   int maximumColumnsExtra_;
    01163   int numberGoodU_;
    01165   int numberGoodL_;
    01167   int maximumPivots_;
    01169   int numberPivots_;
    01172   CoinBigIndex totalElements_;
    01174   CoinBigIndex factorElements_;
    01176   CoinIntArrayWithLength pivotColumn_;
    01178   CoinIntArrayWithLength permute_;
    01180   CoinIntArrayWithLength permuteBack_;
    01182   CoinIntArrayWithLength pivotColumnBack_;
    01184   int status_;
    01185
    01190   //int increasingRows_;
    01191
    01193   int numberTrials_;
    01195   CoinBigIndexArrayWithLength startRowU_;
    01196
    01198   CoinIntArrayWithLength numberInRow_;
    01199
    01201   CoinIntArrayWithLength numberInColumn_;
    01202
    01204   CoinIntArrayWithLength numberInColumnPlus_;
    01205
    01208   CoinIntArrayWithLength firstCount_;
    01209
    01211   CoinIntArrayWithLength nextCount_;
    01212
    01214   CoinIntArrayWithLength lastCount_;
    01215
    01217   CoinIntArrayWithLength nextColumn_;
    01218
    01220   CoinIntArrayWithLength lastColumn_;
    01221
    01223   CoinIntArrayWithLength nextRow_;
    01224
    01226   CoinIntArrayWithLength lastRow_;
    01227
    01229   CoinIntArrayWithLength saveColumn_;
    01230
    01232   CoinIntArrayWithLength markRow_;
    01233
    01235   int messageLevel_;
    01236
    01238   int biggerDimension_;
    01239
    01241   CoinIntArrayWithLength indexColumnU_;
    01242
    01244   CoinIntArrayWithLength pivotRowL_;
    01245
    01247   CoinDoubleArrayWithLength pivotRegion_;
    01248
    01250   int numberSlacks_;
    01251
    01253   int numberU_;
    01254
    01256   CoinBigIndex maximumU_;
    01257
    01259   //int baseU_;
    01260
    01262   CoinBigIndex lengthU_;
    01263
    01265   CoinBigIndex lengthAreaU_;
    01266
    01268   CoinDoubleArrayWithLength elementU_;
    01269
    01271   CoinIntArrayWithLength indexRowU_;
    01272
    01274   CoinBigIndexArrayWithLength startColumnU_;
    01275
    01277   CoinBigIndexArrayWithLength convertRowToColumnU_;
    01278
    01280   CoinBigIndex numberL_;
    01281
    01283   CoinBigIndex baseL_;
    01284
    01286   CoinBigIndex lengthL_;
    01287
    01289   CoinBigIndex lengthAreaL_;
    01290
    01292   CoinDoubleArrayWithLength elementL_;
    01293
    01295   CoinIntArrayWithLength indexRowL_;
    01296
    01298   CoinBigIndexArrayWithLength startColumnL_;
    01299
    01301   bool doForrestTomlin_;
    01302
    01304   int numberR_;
    01305
    01307   CoinBigIndex lengthR_;
    01308
    01310   CoinBigIndex lengthAreaR_;
    01311
    01313   double *elementR_;
    01314
    01316   int *indexRowR_;
    01317
    01319   CoinBigIndexArrayWithLength startColumnR_;
    01320
    01322   double  * denseArea_;
    01323
    01325   int * densePermute_;
    01326
    01328   int numberDense_;
    01329
    01331   int denseThreshold_;
    01332
    01334   CoinDoubleArrayWithLength workArea_;
    01335
    01337   CoinUnsignedIntArrayWithLength workArea2_;
    01338
    01340   CoinBigIndex numberCompressions_;
    01341
    01343   mutable double ftranCountInput_;
    01344   mutable double ftranCountAfterL_;
    01345   mutable double ftranCountAfterR_;
    01346   mutable double ftranCountAfterU_;
    01347   mutable double btranCountInput_;
    01348   mutable double btranCountAfterU_;
    01349   mutable double btranCountAfterR_;
    01350   mutable double btranCountAfterL_;
    01351
    01353   mutable int numberFtranCounts_;
    01354   mutable int numberBtranCounts_;
    01355
    01357   double ftranAverageAfterL_;
    01358   double ftranAverageAfterR_;
    01359   double ftranAverageAfterU_;
    01360   double btranAverageAfterU_;
    01361   double btranAverageAfterR_;
    01362   double btranAverageAfterL_;
    01363
    01365   mutable bool collectStatistics_;
    01366
    01368   int sparseThreshold_;
    01369
    01371   int sparseThreshold2_;
    01372
    01374   CoinBigIndexArrayWithLength startRowL_;
    01375
    01377   CoinIntArrayWithLength indexColumnL_;
    01378
    01380   CoinDoubleArrayWithLength elementByRowL_;
    01381
    01383   mutable CoinIntArrayWithLength sparse_;
    01387   int biasLU_;
    01393   int persistenceFlag_;
    01395 };
    01396 // Dense coding
    01397 #ifdef COIN_HAS_LAPACK
    01398 #define DENSE_CODE 1
    01399 /* Type of Fortran integer translated into C */
    01400 #ifndef ipfint
    01401 //typedef ipfint FORTRAN_INTEGER_TYPE ;
    01402 typedef int ipfint;
    01403 typedef const int cipfint;
    01404 #endif
    01405 #endif
    01406 #endif