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

- File List
- File Members

# /home/coin/SVN-release/Cgl-0.55.0/CoinUtils/src/CoinFactorization.hpp

Go to the documentation of this file.

    00001 /* $Id: CoinFactorization.hpp 1215 2009-11-05 11:03:04Z forrest $ */
    00002 // Copyright (C) 2002, International Business Machines
    00003 // Corporation and others.  All Rights Reserved.
    00004
    00005 /*
    00006    Authors
    00007
    00008    John Forrest
    00009
    00010  */
    00011 #ifndef CoinFactorization_H
    00012 #define CoinFactorization_H
    00013 //#define COIN_ONE_ETA_COPY 100
    00014
    00015 #include <iostream>
    00016 #include <string>
    00017 #include <cassert>
    00018 #include <cstdio>
    00019 #include "CoinFinite.hpp"
    00020 #include "CoinIndexedVector.hpp"
    00021 class CoinPackedMatrix;
    00047 class CoinFactorization {
    00048    friend void CoinFactorizationUnitTest( const std::string & mpsDir );
    00049
    00050 public:
    00051
    00054
    00055     CoinFactorization (  );
    00057   CoinFactorization ( const CoinFactorization &other);
    00058
    00060    ~CoinFactorization (  );
    00062   void almostDestructor();
    00064   void show_self (  ) const;
    00066   int saveFactorization (const char * file  ) const;
    00070   int restoreFactorization (const char * file  , bool factor=false) ;
    00072   void sort (  ) const;
    00074     CoinFactorization & operator = ( const CoinFactorization & other );
    00076
    00086   int factorize ( const CoinPackedMatrix & matrix,
    00087                   int rowIsBasic[], int columnIsBasic[] ,
    00088                   double areaFactor = 0.0 );
    00099   int factorize ( int numberRows,
    00100                   int numberColumns,
    00101                   CoinBigIndex numberElements,
    00102                   CoinBigIndex maximumL,
    00103                   CoinBigIndex maximumU,
    00104                   const int indicesRow[],
    00105                   const int indicesColumn[], const double elements[] ,
    00106                   int permutation[],
    00107                   double areaFactor = 0.0);
    00112   int factorizePart1 ( int numberRows,
    00113                        int numberColumns,
    00114                        CoinBigIndex estimateNumberElements,
    00115                        int * indicesRow[],
    00116                        int * indicesColumn[],
    00117                        CoinFactorizationDouble * elements[],
    00118                   double areaFactor = 0.0);
    00125   int factorizePart2 (int permutation[],int exactNumberElements);
    00127   double conditionNumber() const;
    00128
    00130
    00133
    00134   inline int status (  ) const {
    00135     return status_;
    00136   }
    00138   inline void setStatus (  int value)
    00139   {  status_=value;  }
    00141   inline int pivots (  ) const {
    00142     return numberPivots_;
    00143   }
    00145   inline void setPivots (  int value )
    00146   { numberPivots_=value; }
    00148   inline int *permute (  ) const {
    00149     return permute_.array();
    00150   }
    00152   inline int *pivotColumn (  ) const {
    00153     return pivotColumn_.array();
    00154   }
    00156   inline CoinFactorizationDouble *pivotRegion (  ) const {
    00157     return pivotRegion_.array();
    00158   }
    00160   inline int *permuteBack (  ) const {
    00161     return permuteBack_.array();
    00162   }
    00165   inline int *pivotColumnBack (  ) const {
    00166     //return firstCount_.array();
    00167     return pivotColumnBack_.array();
    00168   }
    00170   inline CoinBigIndex * startRowL() const
    00171   { return startRowL_.array();}
    00172
    00174   inline CoinBigIndex * startColumnL() const
    00175   { return startColumnL_.array();}
    00176
    00178   inline int * indexColumnL() const
    00179   { return indexColumnL_.array();}
    00180
    00182   inline int * indexRowL() const
    00183   { return indexRowL_.array();}
    00184
    00186   inline CoinFactorizationDouble * elementByRowL() const
    00187   { return elementByRowL_.array();}
    00188
    00190   inline int numberRowsExtra (  ) const {
    00191     return numberRowsExtra_;
    00192   }
    00194   inline void setNumberRows(int value)
    00195   { numberRows_ = value; }
    00197   inline int numberRows (  ) const {
    00198     return numberRows_;
    00199   }
    00201   inline CoinBigIndex numberL() const
    00202   { return numberL_;}
    00203
    00205   inline CoinBigIndex baseL() const
    00206   { return baseL_;}
    00208   inline int maximumRowsExtra (  ) const {
    00209     return maximumRowsExtra_;
    00210   }
    00212   inline int numberColumns (  ) const {
    00213     return numberColumns_;
    00214   }
    00216   inline int numberElements (  ) const {
    00217     return totalElements_;
    00218   }
    00220   inline int numberForrestTomlin (  ) const {
    00221     return numberInColumn_.array()[numberColumnsExtra_];
    00222   }
    00224   inline int numberGoodColumns (  ) const {
    00225     return numberGoodU_;
    00226   }
    00228   inline double areaFactor (  ) const {
    00229     return areaFactor_;
    00230   }
    00231   inline void areaFactor ( double value ) {
    00232     areaFactor_=value;
    00233   }
    00235   double adjustedAreaFactor() const;
    00237   inline void relaxAccuracyCheck(double value)
    00238   { relaxCheck_ = value;}
    00239   inline double getAccuracyCheck() const
    00240   { return relaxCheck_;}
    00242   inline int messageLevel (  ) const {
    00243     return messageLevel_ ;
    00244   }
    00245   void messageLevel (  int value );
    00247   inline int maximumPivots (  ) const {
    00248     return maximumPivots_ ;
    00249   }
    00250   void maximumPivots (  int value );
    00251
    00253   inline int denseThreshold() const
    00254     { return denseThreshold_;}
    00256   inline void setDenseThreshold(int value)
    00257     { denseThreshold_ = value;}
    00259   inline double pivotTolerance (  ) const {
    00260     return pivotTolerance_ ;
    00261   }
    00262   void pivotTolerance (  double value );
    00264   inline double zeroTolerance (  ) const {
    00265     return zeroTolerance_ ;
    00266   }
    00267   void zeroTolerance (  double value );
    00268 #ifndef COIN_FAST_CODE
    00270   inline double slackValue (  ) const {
    00271     return slackValue_ ;
    00272   }
    00273   void slackValue (  double value );
    00274 #endif
    00276   double maximumCoefficient() const;
    00278   inline bool forrestTomlin() const
    00279   { return doForrestTomlin_;}
    00280   inline void setForrestTomlin(bool value)
    00281   { doForrestTomlin_=value;}
    00283   inline bool spaceForForrestTomlin() const
    00284   {
    00285     CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
    00286     CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
    00287     return (space>=0)&&doForrestTomlin_;
    00288   }
    00290
    00293
    00295   inline int numberDense() const
    00296   { return numberDense_;}
    00297
    00299   inline CoinBigIndex numberElementsU (  ) const {
    00300     return lengthU_;
    00301   }
    00303   inline void setNumberElementsU(CoinBigIndex value)
    00304   { lengthU_ = value; }
    00306   inline CoinBigIndex lengthAreaU (  ) const {
    00307     return lengthAreaU_;
    00308   }
    00310   inline CoinBigIndex numberElementsL (  ) const {
    00311     return lengthL_;
    00312   }
    00314   inline CoinBigIndex lengthAreaL (  ) const {
    00315     return lengthAreaL_;
    00316   }
    00318   inline CoinBigIndex numberElementsR (  ) const {
    00319     return lengthR_;
    00320   }
    00322   inline CoinBigIndex numberCompressions() const
    00323   { return numberCompressions_;}
    00325   inline int * numberInRow() const
    00326   { return numberInRow_.array();}
    00328   inline int * numberInColumn() const
    00329   { return numberInColumn_.array();}
    00331   inline CoinFactorizationDouble * elementU() const
    00332   { return elementU_.array();}
    00334   inline int * indexRowU() const
    00335   { return indexRowU_.array();}
    00337   inline CoinBigIndex * startColumnU() const
    00338   { return startColumnU_.array();}
    00340   inline int maximumColumnsExtra()
    00341   { return maximumColumnsExtra_;}
    00345   inline int biasLU() const
    00346   { return biasLU_;}
    00347   inline void setBiasLU(int value)
    00348   { biasLU_=value;}
    00354   inline int persistenceFlag() const
    00355   { return persistenceFlag_;}
    00356   void setPersistenceFlag(int value);
    00358
    00361
    00369   int replaceColumn ( CoinIndexedVector * regionSparse,
    00370                       int pivotRow,
    00371                       double pivotCheck ,
    00372                       bool checkBeforeModifying=false,
    00373                       double acceptablePivot=1.0e-8);
    00378   void replaceColumnU ( CoinIndexedVector * regionSparse,
    00379                         CoinBigIndex * deleted,
    00380                         int internalPivotRow);
    00382
    00392   int updateColumnFT ( CoinIndexedVector * regionSparse,
    00393                        CoinIndexedVector * regionSparse2);
    00396   int updateColumn ( CoinIndexedVector * regionSparse,
    00397                      CoinIndexedVector * regionSparse2,
    00398                      bool noPermute=false) const;
    00404   int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
    00405                            CoinIndexedVector * regionSparse2,
    00406                            CoinIndexedVector * regionSparse3,
    00407                            bool noPermuteRegion3=false) ;
    00412   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
    00413                               CoinIndexedVector * regionSparse2) const;
    00415   void goSparse();
    00417   inline int sparseThreshold ( ) const
    00418   { return sparseThreshold_;}
    00420   void sparseThreshold ( int value );
    00422
    00423
    00427
    00428   inline void clearArrays()
    00429   { gutsOfDestructor();}
    00431
    00434
    00437   int add ( CoinBigIndex numberElements,
    00438                int indicesRow[],
    00439                int indicesColumn[], double elements[] );
    00440
    00443   int addColumn ( CoinBigIndex numberElements,
    00444                      int indicesRow[], double elements[] );
    00445
    00448   int addRow ( CoinBigIndex numberElements,
    00449                   int indicesColumn[], double elements[] );
    00450
    00452   int deleteColumn ( int Row );
    00454   int deleteRow ( int Row );
    00455
    00459   int replaceRow ( int whichRow, int numberElements,
    00460                       const int indicesColumn[], const double elements[] );
    00462   void emptyRows(int numberToEmpty, const int which[]);
    00464
    00465
    00466   void checkSparse();
    00468   inline bool collectStatistics() const
    00469   { return collectStatistics_;}
    00471   inline void setCollectStatistics(bool onOff) const
    00472   { collectStatistics_ = onOff;}
    00474   void gutsOfDestructor(int type=1);
    00476   void gutsOfInitialize(int type);
    00477   void gutsOfCopy(const CoinFactorization &other);
    00478
    00480   void resetStatistics();
    00481
    00482
    00484
    00486
    00487   void getAreas ( int numberRows,
    00488                   int numberColumns,
    00489                   CoinBigIndex maximumL,
    00490                   CoinBigIndex maximumU );
    00491
    00494   void preProcess ( int state,
    00495                     int possibleDuplicates = -1 );
    00497   int factor (  );
    00498 protected:
    00501   int factorSparse (  );
    00504   int factorSparseSmall (  );
    00507   int factorSparseLarge (  );
    00510   int factorDense (  );
    00511
    00513   bool pivotOneOtherRow ( int pivotRow,
    00514                           int pivotColumn );
    00516   bool pivotRowSingleton ( int pivotRow,
    00517                            int pivotColumn );
    00519   bool pivotColumnSingleton ( int pivotRow,
    00520                               int pivotColumn );
    00521
    00526   bool getColumnSpace ( int iColumn,
    00527                         int extraNeeded );
    00528
    00531   bool reorderU();
    00535   bool getColumnSpaceIterateR ( int iColumn, double value,
    00536                                int iRow);
    00542   CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
    00543                                int iRow);
    00547   bool getRowSpace ( int iRow, int extraNeeded );
    00548
    00552   bool getRowSpaceIterate ( int iRow,
    00553                             int extraNeeded );
    00555   void checkConsistency (  );
    00557   inline void addLink ( int index, int count ) {
    00558     int *nextCount = nextCount_.array();
    00559     int *firstCount = firstCount_.array();
    00560     int *lastCount = lastCount_.array();
    00561     int next = firstCount[count];
    00562       lastCount[index] = -2 - count;
    00563     if ( next < 0 ) {
    00564       //first with that count
    00565       firstCount[count] = index;
    00566       nextCount[index] = -1;
    00567     } else {
    00568       firstCount[count] = index;
    00569       nextCount[index] = next;
    00570       lastCount[next] = index;
    00571   }}
    00573   inline void deleteLink ( int index ) {
    00574     int *nextCount = nextCount_.array();
    00575     int *firstCount = firstCount_.array();
    00576     int *lastCount = lastCount_.array();
    00577     int next = nextCount[index];
    00578     int last = lastCount[index];
    00579     if ( last >= 0 ) {
    00580       nextCount[last] = next;
    00581     } else {
    00582       int count = -last - 2;
    00583
    00584       firstCount[count] = next;
    00585     }
    00586     if ( next >= 0 ) {
    00587       lastCount[next] = last;
    00588     }
    00589     nextCount[index] = -2;
    00590     lastCount[index] = -2;
    00591     return;
    00592   }
    00594   void separateLinks(int count,bool rowsFirst);
    00596   void cleanup (  );
    00597
    00599   void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
    00601   void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
    00603   void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
    00605   void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
    00606
    00608   void updateColumnR ( CoinIndexedVector * region ) const;
    00611   void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
    00612
    00614   void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
    00615
    00617   void updateColumnUSparse ( CoinIndexedVector * regionSparse,
    00618                              int * indexIn) const;
    00620   void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
    00621                                int * indexIn) const;
    00623   int updateColumnUDensish ( double * COIN_RESTRICT region,
    00624                              int * COIN_RESTRICT regionIndex) const;
    00626   void updateTwoColumnsUDensish (
    00627                                  int & numberNonZero1,
    00628                                  double * COIN_RESTRICT region1,
    00629                                  int * COIN_RESTRICT index1,
    00630                                  int & numberNonZero2,
    00631                                  double * COIN_RESTRICT region2,
    00632                                  int * COIN_RESTRICT index2) const;
    00634   void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
    00636   void permuteBack ( CoinIndexedVector * regionSparse,
    00637                      CoinIndexedVector * outVector) const;
    00638
    00640   void updateColumnTransposePFI ( CoinIndexedVector * region) const;
    00643   void updateColumnTransposeU ( CoinIndexedVector * region,
    00644                                 int smallestIndex) const;
    00647   void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
    00648                                         int smallestIndex) const;
    00651   void updateColumnTransposeUDensish ( CoinIndexedVector * region,
    00652                                        int smallestIndex) const;
    00655   void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
    00658   void updateColumnTransposeUByColumn ( CoinIndexedVector * region,
    00659                                         int smallestIndex) const;
    00660
    00662   void updateColumnTransposeR ( CoinIndexedVector * region ) const;
    00664   void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
    00666   void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
    00667
    00669   void updateColumnTransposeL ( CoinIndexedVector * region ) const;
    00671   void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
    00673   void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
    00675   void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
    00677   void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
    00678 public:
    00683   int replaceColumnPFI ( CoinIndexedVector * regionSparse,
    00684                          int pivotRow, double alpha);
    00685 protected:
    00688   int checkPivot(double saveFromU, double oldPivot) const;
    00689   /********************************* START LARGE TEMPLATE ********/
    00690 #ifdef INT_IS_8
    00691 #define COINFACTORIZATION_BITS_PER_INT 64
    00692 #define COINFACTORIZATION_SHIFT_PER_INT 6
    00693 #define COINFACTORIZATION_MASK_PER_INT 0x3f
    00694 #else
    00695 #define COINFACTORIZATION_BITS_PER_INT 32
    00696 #define COINFACTORIZATION_SHIFT_PER_INT 5
    00697 #define COINFACTORIZATION_MASK_PER_INT 0x1f
    00698 #endif
    00699   template <class T>  inline bool
    00700   pivot ( int pivotRow,
    00701           int pivotColumn,
    00702           CoinBigIndex pivotRowPosition,
    00703           CoinBigIndex pivotColumnPosition,
    00704           CoinFactorizationDouble work[],
    00705           unsigned int workArea2[],
    00706           int increment2,
    00707           T markRow[] ,
    00708           int largeInteger)
    00709 {
    00710   int *indexColumnU = indexColumnU_.array();
    00711   CoinBigIndex *startColumnU = startColumnU_.array();
    00712   int *numberInColumn = numberInColumn_.array();
    00713   CoinFactorizationDouble *elementU = elementU_.array();
    00714   int *indexRowU = indexRowU_.array();
    00715   CoinBigIndex *startRowU = startRowU_.array();
    00716   int *numberInRow = numberInRow_.array();
    00717   CoinFactorizationDouble *elementL = elementL_.array();
    00718   int *indexRowL = indexRowL_.array();
    00719   int *saveColumn = saveColumn_.array();
    00720   int *nextRow = nextRow_.array();
    00721   int *lastRow = lastRow_.array() ;
    00722
    00723   //store pivot columns (so can easily compress)
    00724   int numberInPivotRow = numberInRow[pivotRow] - 1;
    00725   CoinBigIndex startColumn = startColumnU[pivotColumn];
    00726   int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
    00727   CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
    00728   int put = 0;
    00729   CoinBigIndex startRow = startRowU[pivotRow];
    00730   CoinBigIndex endRow = startRow + numberInPivotRow + 1;
    00731
    00732   if ( pivotColumnPosition < 0 ) {
    00733     for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    00734       int iColumn = indexColumnU[pivotColumnPosition];
    00735       if ( iColumn != pivotColumn ) {
    00736         saveColumn[put++] = iColumn;
    00737       } else {
    00738         break;
    00739       }
    00740     }
    00741   } else {
    00742     for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
    00743       saveColumn[put++] = indexColumnU[i];
    00744     }
    00745   }
    00746   assert (pivotColumnPosition<endRow);
    00747   assert (indexColumnU[pivotColumnPosition]==pivotColumn);
    00748   pivotColumnPosition++;
    00749   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    00750     saveColumn[put++] = indexColumnU[pivotColumnPosition];
    00751   }
    00752   //take out this bit of indexColumnU
    00753   int next = nextRow[pivotRow];
    00754   int last = lastRow[pivotRow];
    00755
    00756   nextRow[last] = next;
    00757   lastRow[next] = last;
    00758   nextRow[pivotRow] = numberGoodU_;     //use for permute
    00759   lastRow[pivotRow] = -2;
    00760   numberInRow[pivotRow] = 0;
    00761   //store column in L, compress in U and take column out
    00762   CoinBigIndex l = lengthL_;
    00763
    00764   if ( l + numberInPivotColumn > lengthAreaL_ ) {
    00765     //need more memory
    00766     if ((messageLevel_&4)!=0)
    00767       printf("more memory needed in middle of invert\n");
    00768     return false;
    00769   }
    00770   //l+=currentAreaL_->elementByColumn-elementL;
    00771   CoinBigIndex lSave = l;
    00772
    00773   CoinBigIndex * startColumnL = startColumnL_.array();
    00774   startColumnL[numberGoodL_] = l;       //for luck and first time
    00775   numberGoodL_++;
    00776   startColumnL[numberGoodL_] = l + numberInPivotColumn;
    00777   lengthL_ += numberInPivotColumn;
    00778   if ( pivotRowPosition < 0 ) {
    00779     for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    00780       int iRow = indexRowU[pivotRowPosition];
    00781       if ( iRow != pivotRow ) {
    00782         indexRowL[l] = iRow;
    00783         elementL[l] = elementU[pivotRowPosition];
    00784         markRow[iRow] = static_cast<T>(l - lSave);
    00785         l++;
    00786         //take out of row list
    00787         CoinBigIndex start = startRowU[iRow];
    00788         CoinBigIndex end = start + numberInRow[iRow];
    00789         CoinBigIndex where = start;
    00790
    00791         while ( indexColumnU[where] != pivotColumn ) {
    00792           where++;
    00793         }                       /* endwhile */
    00794 #if DEBUG_COIN
    00795         if ( where >= end ) {
    00796           abort (  );
    00797         }
    00798 #endif
    00799         indexColumnU[where] = indexColumnU[end - 1];
    00800         numberInRow[iRow]--;
    00801       } else {
    00802         break;
    00803       }
    00804     }
    00805   } else {
    00806     CoinBigIndex i;
    00807
    00808     for ( i = startColumn; i < pivotRowPosition; i++ ) {
    00809       int iRow = indexRowU[i];
    00810
    00811       markRow[iRow] = static_cast<T>(l - lSave);
    00812       indexRowL[l] = iRow;
    00813       elementL[l] = elementU[i];
    00814       l++;
    00815       //take out of row list
    00816       CoinBigIndex start = startRowU[iRow];
    00817       CoinBigIndex end = start + numberInRow[iRow];
    00818       CoinBigIndex where = start;
    00819
    00820       while ( indexColumnU[where] != pivotColumn ) {
    00821         where++;
    00822       }                         /* endwhile */
    00823 #if DEBUG_COIN
    00824       if ( where >= end ) {
    00825         abort (  );
    00826       }
    00827 #endif
    00828       indexColumnU[where] = indexColumnU[end - 1];
    00829       numberInRow[iRow]--;
    00830       assert (numberInRow[iRow]>=0);
    00831     }
    00832   }
    00833   assert (pivotRowPosition<endColumn);
    00834   assert (indexRowU[pivotRowPosition]==pivotRow);
    00835   CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
    00836   CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
    00837
    00838   pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
    00839   pivotRowPosition++;
    00840   for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    00841     int iRow = indexRowU[pivotRowPosition];
    00842
    00843     markRow[iRow] = static_cast<T>(l - lSave);
    00844     indexRowL[l] = iRow;
    00845     elementL[l] = elementU[pivotRowPosition];
    00846     l++;
    00847     //take out of row list
    00848     CoinBigIndex start = startRowU[iRow];
    00849     CoinBigIndex end = start + numberInRow[iRow];
    00850     CoinBigIndex where = start;
    00851
    00852     while ( indexColumnU[where] != pivotColumn ) {
    00853       where++;
    00854     }                           /* endwhile */
    00855 #if DEBUG_COIN
    00856     if ( where >= end ) {
    00857       abort (  );
    00858     }
    00859 #endif
    00860     indexColumnU[where] = indexColumnU[end - 1];
    00861     numberInRow[iRow]--;
    00862     assert (numberInRow[iRow]>=0);
    00863   }
    00864   markRow[pivotRow] = static_cast<T>(largeInteger);
    00865   //compress pivot column (move pivot to front including saved)
    00866   numberInColumn[pivotColumn] = 0;
    00867   //use end of L for temporary space
    00868   int *indexL = &indexRowL[lSave];
    00869   CoinFactorizationDouble *multipliersL = &elementL[lSave];
    00870
    00871   //adjust
    00872   int j;
    00873
    00874   for ( j = 0; j < numberInPivotColumn; j++ ) {
    00875     multipliersL[j] *= pivotMultiplier;
    00876   }
    00877   //zero out fill
    00878   CoinBigIndex iErase;
    00879   for ( iErase = 0; iErase < increment2 * numberInPivotRow;
    00880         iErase++ ) {
    00881     workArea2[iErase] = 0;
    00882   }
    00883   CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
    00884   unsigned int *temp2 = workArea2;
    00885   int * nextColumn = nextColumn_.array();
    00886
    00887   //pack down and move to work
    00888   int jColumn;
    00889   for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    00890     int iColumn = saveColumn[jColumn];
    00891     CoinBigIndex startColumn = startColumnU[iColumn];
    00892     CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
    00893     int iRow = indexRowU[startColumn];
    00894     CoinFactorizationDouble value = elementU[startColumn];
    00895     double largest;
    00896     CoinBigIndex put = startColumn;
    00897     CoinBigIndex positionLargest = -1;
    00898     CoinFactorizationDouble thisPivotValue = 0.0;
    00899
    00900     //compress column and find largest not updated
    00901     bool checkLargest;
    00902     int mark = markRow[iRow];
    00903
    00904     if ( mark == largeInteger+1 ) {
    00905       largest = fabs ( value );
    00906       positionLargest = put;
    00907       put++;
    00908       checkLargest = false;
    00909     } else {
    00910       //need to find largest
    00911       largest = 0.0;
    00912       checkLargest = true;
    00913       if ( mark != largeInteger ) {
    00914         //will be updated
    00915         work[mark] = value;
    00916         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    00917         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    00918
    00919         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
    00920         added--;
    00921       } else {
    00922         thisPivotValue = value;
    00923       }
    00924     }
    00925     CoinBigIndex i;
    00926     for ( i = startColumn + 1; i < endColumn; i++ ) {
    00927       iRow = indexRowU[i];
    00928       value = elementU[i];
    00929       int mark = markRow[iRow];
    00930
    00931       if ( mark == largeInteger+1 ) {
    00932         //keep
    00933         indexRowU[put] = iRow;
    00934         elementU[put] = value;
    00935         if ( checkLargest ) {
    00936           double absValue = fabs ( value );
    00937
    00938           if ( absValue > largest ) {
    00939             largest = absValue;
    00940             positionLargest = put;
    00941           }
    00942         }
    00943         put++;
    00944       } else if ( mark != largeInteger ) {
    00945         //will be updated
    00946         work[mark] = value;
    00947         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    00948         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    00949
    00950         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
    00951         added--;
    00952       } else {
    00953         thisPivotValue = value;
    00954       }
    00955     }
    00956     //slot in pivot
    00957     elementU[put] = elementU[startColumn];
    00958     indexRowU[put] = indexRowU[startColumn];
    00959     if ( positionLargest == startColumn ) {
    00960       positionLargest = put;    //follow if was largest
    00961     }
    00962     put++;
    00963     elementU[startColumn] = thisPivotValue;
    00964     indexRowU[startColumn] = pivotRow;
    00965     //clean up counts
    00966     startColumn++;
    00967     numberInColumn[iColumn] = put - startColumn;
    00968     int * numberInColumnPlus = numberInColumnPlus_.array();
    00969     numberInColumnPlus[iColumn]++;
    00970     startColumnU[iColumn]++;
    00971     //how much space have we got
    00972     int next = nextColumn[iColumn];
    00973     CoinBigIndex space;
    00974
    00975     space = startColumnU[next] - put - numberInColumnPlus[next];
    00976     //assume no zero elements
    00977     if ( numberInPivotColumn > space ) {
    00978       //getColumnSpace also moves fixed part
    00979       if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
    00980         return false;
    00981       }
    00982       //redo starts
    00983       positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
    00984       startColumn = startColumnU[iColumn];
    00985       put = startColumn + numberInColumn[iColumn];
    00986     }
    00987     double tolerance = zeroTolerance_;
    00988
    00989     int *nextCount = nextCount_.array();
    00990     for ( j = 0; j < numberInPivotColumn; j++ ) {
    00991       value = work[j] - thisPivotValue * multipliersL[j];
    00992       double absValue = fabs ( value );
    00993
    00994       if ( absValue > tolerance ) {
    00995         work[j] = 0.0;
    00996         elementU[put] = value;
    00997         indexRowU[put] = indexL[j];
    00998         if ( absValue > largest ) {
    00999           largest = absValue;
    01000           positionLargest = put;
    01001         }
    01002         put++;
    01003       } else {
    01004         work[j] = 0.0;
    01005         added--;
    01006         int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01007         int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01008
    01009         if ( temp2[word] & ( 1 << bit ) ) {
    01010           //take out of row list
    01011           iRow = indexL[j];
    01012           CoinBigIndex start = startRowU[iRow];
    01013           CoinBigIndex end = start + numberInRow[iRow];
    01014           CoinBigIndex where = start;
    01015
    01016           while ( indexColumnU[where] != iColumn ) {
    01017             where++;
    01018           }                     /* endwhile */
    01019 #if DEBUG_COIN
    01020           if ( where >= end ) {
    01021             abort (  );
    01022           }
    01023 #endif
    01024           indexColumnU[where] = indexColumnU[end - 1];
    01025           numberInRow[iRow]--;
    01026         } else {
    01027           //make sure won't be added
    01028           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01029           int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01030
    01031           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    01032         }
    01033       }
    01034     }
    01035     numberInColumn[iColumn] = put - startColumn;
    01036     //move largest
    01037     if ( positionLargest >= 0 ) {
    01038       value = elementU[positionLargest];
    01039       iRow = indexRowU[positionLargest];
    01040       elementU[positionLargest] = elementU[startColumn];
    01041       indexRowU[positionLargest] = indexRowU[startColumn];
    01042       elementU[startColumn] = value;
    01043       indexRowU[startColumn] = iRow;
    01044     }
    01045     //linked list for column
    01046     if ( nextCount[iColumn + numberRows_] != -2 ) {
    01047       //modify linked list
    01048       deleteLink ( iColumn + numberRows_ );
    01049       addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
    01050     }
    01051     temp2 += increment2;
    01052   }
    01053   //get space for row list
    01054   unsigned int *putBase = workArea2;
    01055   int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
    01056   int i = 0;
    01057
    01058   // do linked lists and update counts
    01059   while ( bigLoops ) {
    01060     bigLoops--;
    01061     int bit;
    01062     for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
    01063       unsigned int *putThis = putBase;
    01064       int iRow = indexL[i];
    01065
    01066       //get space
    01067       int number = 0;
    01068       int jColumn;
    01069
    01070       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01071         unsigned int test = *putThis;
    01072
    01073         putThis += increment2;
    01074         test = 1 - ( ( test >> bit ) & 1 );
    01075         number += test;
    01076       }
    01077       int next = nextRow[iRow];
    01078       CoinBigIndex space;
    01079
    01080       space = startRowU[next] - startRowU[iRow];
    01081       number += numberInRow[iRow];
    01082       if ( space < number ) {
    01083         if ( !getRowSpace ( iRow, number ) ) {
    01084           return false;
    01085         }
    01086       }
    01087       // now do
    01088       putThis = putBase;
    01089       next = nextRow[iRow];
    01090       number = numberInRow[iRow];
    01091       CoinBigIndex end = startRowU[iRow] + number;
    01092       int saveIndex = indexColumnU[startRowU[next]];
    01093
    01094       //add in
    01095       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01096         unsigned int test = *putThis;
    01097
    01098         putThis += increment2;
    01099         test = 1 - ( ( test >> bit ) & 1 );
    01100         indexColumnU[end] = saveColumn[jColumn];
    01101         end += test;
    01102       }
    01103       //put back next one in case zapped
    01104       indexColumnU[startRowU[next]] = saveIndex;
    01105       markRow[iRow] = static_cast<T>(largeInteger+1);
    01106       number = end - startRowU[iRow];
    01107       numberInRow[iRow] = number;
    01108       deleteLink ( iRow );
    01109       addLink ( iRow, number );
    01110     }
    01111     putBase++;
    01112   }                             /* endwhile */
    01113   int bit;
    01114
    01115   for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
    01116     unsigned int *putThis = putBase;
    01117     int iRow = indexL[i];
    01118
    01119     //get space
    01120     int number = 0;
    01121     int jColumn;
    01122
    01123     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01124       unsigned int test = *putThis;
    01125
    01126       putThis += increment2;
    01127       test = 1 - ( ( test >> bit ) & 1 );
    01128       number += test;
    01129     }
    01130     int next = nextRow[iRow];
    01131     CoinBigIndex space;
    01132
    01133     space = startRowU[next] - startRowU[iRow];
    01134     number += numberInRow[iRow];
    01135     if ( space < number ) {
    01136       if ( !getRowSpace ( iRow, number ) ) {
    01137         return false;
    01138       }
    01139     }
    01140     // now do
    01141     putThis = putBase;
    01142     next = nextRow[iRow];
    01143     number = numberInRow[iRow];
    01144     CoinBigIndex end = startRowU[iRow] + number;
    01145     int saveIndex;
    01146
    01147     saveIndex = indexColumnU[startRowU[next]];
    01148
    01149     //add in
    01150     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01151       unsigned int test = *putThis;
    01152
    01153       putThis += increment2;
    01154       test = 1 - ( ( test >> bit ) & 1 );
    01155
    01156       indexColumnU[end] = saveColumn[jColumn];
    01157       end += test;
    01158     }
    01159     indexColumnU[startRowU[next]] = saveIndex;
    01160     markRow[iRow] = static_cast<T>(largeInteger+1);
    01161     number = end - startRowU[iRow];
    01162     numberInRow[iRow] = number;
    01163     deleteLink ( iRow );
    01164     addLink ( iRow, number );
    01165   }
    01166   markRow[pivotRow] = static_cast<T>(largeInteger+1);
    01167   //modify linked list for pivots
    01168   deleteLink ( pivotRow );
    01169   deleteLink ( pivotColumn + numberRows_ );
    01170   totalElements_ += added;
    01171   return true;
    01172 }
    01173
    01174   /********************************* END LARGE TEMPLATE ********/
    01176
    01177 protected:
    01178
    01181
    01182   double pivotTolerance_;
    01184   double zeroTolerance_;
    01185 #ifndef COIN_FAST_CODE
    01187   double slackValue_;
    01188 #else
    01189 #ifndef slackValue_
    01190 #define slackValue_ -1.0
    01191 #endif
    01192 #endif
    01194   double areaFactor_;
    01196   double relaxCheck_;
    01198   int numberRows_;
    01200   int numberRowsExtra_;
    01202   int maximumRowsExtra_;
    01204   int numberColumns_;
    01206   int numberColumnsExtra_;
    01208   int maximumColumnsExtra_;
    01210   int numberGoodU_;
    01212   int numberGoodL_;
    01214   int maximumPivots_;
    01216   int numberPivots_;
    01219   CoinBigIndex totalElements_;
    01221   CoinBigIndex factorElements_;
    01223   CoinIntArrayWithLength pivotColumn_;
    01225   CoinIntArrayWithLength permute_;
    01227   CoinIntArrayWithLength permuteBack_;
    01229   CoinIntArrayWithLength pivotColumnBack_;
    01231   int status_;
    01232
    01237   //int increasingRows_;
    01238
    01240   int numberTrials_;
    01242   CoinBigIndexArrayWithLength startRowU_;
    01243
    01245   CoinIntArrayWithLength numberInRow_;
    01246
    01248   CoinIntArrayWithLength numberInColumn_;
    01249
    01251   CoinIntArrayWithLength numberInColumnPlus_;
    01252
    01255   CoinIntArrayWithLength firstCount_;
    01256
    01258   CoinIntArrayWithLength nextCount_;
    01259
    01261   CoinIntArrayWithLength lastCount_;
    01262
    01264   CoinIntArrayWithLength nextColumn_;
    01265
    01267   CoinIntArrayWithLength lastColumn_;
    01268
    01270   CoinIntArrayWithLength nextRow_;
    01271
    01273   CoinIntArrayWithLength lastRow_;
    01274
    01276   CoinIntArrayWithLength saveColumn_;
    01277
    01279   CoinIntArrayWithLength markRow_;
    01280
    01282   int messageLevel_;
    01283
    01285   int biggerDimension_;
    01286
    01288   CoinIntArrayWithLength indexColumnU_;
    01289
    01291   CoinIntArrayWithLength pivotRowL_;
    01292
    01294   CoinFactorizationDoubleArrayWithLength pivotRegion_;
    01295
    01297   int numberSlacks_;
    01298
    01300   int numberU_;
    01301
    01303   CoinBigIndex maximumU_;
    01304
    01306   //int baseU_;
    01307
    01309   CoinBigIndex lengthU_;
    01310
    01312   CoinBigIndex lengthAreaU_;
    01313
    01315   CoinFactorizationDoubleArrayWithLength elementU_;
    01316
    01318   CoinIntArrayWithLength indexRowU_;
    01319
    01321   CoinBigIndexArrayWithLength startColumnU_;
    01322
    01324   CoinBigIndexArrayWithLength convertRowToColumnU_;
    01325
    01327   CoinBigIndex numberL_;
    01328
    01330   CoinBigIndex baseL_;
    01331
    01333   CoinBigIndex lengthL_;
    01334
    01336   CoinBigIndex lengthAreaL_;
    01337
    01339   CoinFactorizationDoubleArrayWithLength elementL_;
    01340
    01342   CoinIntArrayWithLength indexRowL_;
    01343
    01345   CoinBigIndexArrayWithLength startColumnL_;
    01346
    01348   bool doForrestTomlin_;
    01349
    01351   int numberR_;
    01352
    01354   CoinBigIndex lengthR_;
    01355
    01357   CoinBigIndex lengthAreaR_;
    01358
    01360   CoinFactorizationDouble *elementR_;
    01361
    01363   int *indexRowR_;
    01364
    01366   CoinBigIndexArrayWithLength startColumnR_;
    01367
    01369   double  * denseArea_;
    01370
    01372   int * densePermute_;
    01373
    01375   int numberDense_;
    01376
    01378   int denseThreshold_;
    01379
    01381   CoinFactorizationDoubleArrayWithLength workArea_;
    01382
    01384   CoinUnsignedIntArrayWithLength workArea2_;
    01385
    01387   CoinBigIndex numberCompressions_;
    01388
    01390   mutable double ftranCountInput_;
    01391   mutable double ftranCountAfterL_;
    01392   mutable double ftranCountAfterR_;
    01393   mutable double ftranCountAfterU_;
    01394   mutable double btranCountInput_;
    01395   mutable double btranCountAfterU_;
    01396   mutable double btranCountAfterR_;
    01397   mutable double btranCountAfterL_;
    01398
    01400   mutable int numberFtranCounts_;
    01401   mutable int numberBtranCounts_;
    01402
    01404   double ftranAverageAfterL_;
    01405   double ftranAverageAfterR_;
    01406   double ftranAverageAfterU_;
    01407   double btranAverageAfterU_;
    01408   double btranAverageAfterR_;
    01409   double btranAverageAfterL_;
    01410
    01412   mutable bool collectStatistics_;
    01413
    01415   int sparseThreshold_;
    01416
    01418   int sparseThreshold2_;
    01419
    01421   CoinBigIndexArrayWithLength startRowL_;
    01422
    01424   CoinIntArrayWithLength indexColumnL_;
    01425
    01427   CoinFactorizationDoubleArrayWithLength elementByRowL_;
    01428
    01430   mutable CoinIntArrayWithLength sparse_;
    01434   int biasLU_;
    01440   int persistenceFlag_;
    01442 };
    01443 // Dense coding
    01444 #ifdef COIN_HAS_LAPACK
    01445 #define DENSE_CODE 1
    01446 /* Type of Fortran integer translated into C */
    01447 #ifndef ipfint
    01448 //typedef ipfint FORTRAN_INTEGER_TYPE ;
    01449 typedef int ipfint;
    01450 typedef const int cipfint;
    01451 #endif
    01452 #endif
    01453 #endif
    01454 // Extra for ugly include
    01455 #ifdef UGLY_COIN_FACTOR_CODING
    01456 #define FAC_UNSET (FAC_SET+1)
    01457 {
    01458   goodPivot=false;
    01459   //store pivot columns (so can easily compress)
    01460   CoinBigIndex startColumnThis = startColumn[iPivotColumn];
    01461   CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
    01462   int put = 0;
    01463   CoinBigIndex startRowThis = startRow[iPivotRow];
    01464   CoinBigIndex endRow = startRowThis + numberDoRow + 1;
    01465   if ( pivotColumnPosition < 0 ) {
    01466     for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    01467       int iColumn = indexColumn[pivotColumnPosition];
    01468       if ( iColumn != iPivotColumn ) {
    01469         saveColumn[put++] = iColumn;
    01470       } else {
    01471         break;
    01472       }
    01473     }
    01474   } else {
    01475     for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
    01476       saveColumn[put++] = indexColumn[i];
    01477     }
    01478   }
    01479   assert (pivotColumnPosition<endRow);
    01480   assert (indexColumn[pivotColumnPosition]==iPivotColumn);
    01481   pivotColumnPosition++;
    01482   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    01483     saveColumn[put++] = indexColumn[pivotColumnPosition];
    01484   }
    01485   //take out this bit of indexColumn
    01486   int next = nextRow[iPivotRow];
    01487   int last = lastRow[iPivotRow];
    01488
    01489   nextRow[last] = next;
    01490   lastRow[next] = last;
    01491   nextRow[iPivotRow] = numberGoodU_;    //use for permute
    01492   lastRow[iPivotRow] = -2;
    01493   numberInRow[iPivotRow] = 0;
    01494   //store column in L, compress in U and take column out
    01495   CoinBigIndex l = lengthL_;
    01496   // **** HORRID coding coming up but a goto seems best!
    01497   {
    01498     if ( l + numberDoColumn > lengthAreaL_ ) {
    01499       //need more memory
    01500       if ((messageLevel_&4)!=0)
    01501         printf("more memory needed in middle of invert\n");
    01502       goto BAD_PIVOT;
    01503     }
    01504     //l+=currentAreaL_->elementByColumn-elementL;
    01505     CoinBigIndex lSave = l;
    01506
    01507     CoinBigIndex * startColumnL = startColumnL_.array();
    01508     startColumnL[numberGoodL_] = l;     //for luck and first time
    01509     numberGoodL_++;
    01510     startColumnL[numberGoodL_] = l + numberDoColumn;
    01511     lengthL_ += numberDoColumn;
    01512     if ( pivotRowPosition < 0 ) {
    01513       for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    01514         int iRow = indexRow[pivotRowPosition];
    01515         if ( iRow != iPivotRow ) {
    01516           indexRowL[l] = iRow;
    01517           elementL[l] = element[pivotRowPosition];
    01518           markRow[iRow] = l - lSave;
    01519           l++;
    01520           //take out of row list
    01521           CoinBigIndex start = startRow[iRow];
    01522           CoinBigIndex end = start + numberInRow[iRow];
    01523           CoinBigIndex where = start;
    01524
    01525           while ( indexColumn[where] != iPivotColumn ) {
    01526             where++;
    01527           }                     /* endwhile */
    01528 #if DEBUG_COIN
    01529           if ( where >= end ) {
    01530             abort (  );
    01531           }
    01532 #endif
    01533           indexColumn[where] = indexColumn[end - 1];
    01534           numberInRow[iRow]--;
    01535         } else {
    01536           break;
    01537         }
    01538       }
    01539     } else {
    01540       CoinBigIndex i;
    01541
    01542       for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
    01543         int iRow = indexRow[i];
    01544
    01545         markRow[iRow] = l - lSave;
    01546         indexRowL[l] = iRow;
    01547         elementL[l] = element[i];
    01548         l++;
    01549         //take out of row list
    01550         CoinBigIndex start = startRow[iRow];
    01551         CoinBigIndex end = start + numberInRow[iRow];
    01552         CoinBigIndex where = start;
    01553
    01554         while ( indexColumn[where] != iPivotColumn ) {
    01555           where++;
    01556         }                               /* endwhile */
    01557 #if DEBUG_COIN
    01558         if ( where >= end ) {
    01559           abort (  );
    01560         }
    01561 #endif
    01562         indexColumn[where] = indexColumn[end - 1];
    01563         numberInRow[iRow]--;
    01564         assert (numberInRow[iRow]>=0);
    01565       }
    01566     }
    01567     assert (pivotRowPosition<endColumn);
    01568     assert (indexRow[pivotRowPosition]==iPivotRow);
    01569     CoinFactorizationDouble pivotElement = element[pivotRowPosition];
    01570     CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
    01571
    01572     pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
    01573     pivotRowPosition++;
    01574     for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    01575       int iRow = indexRow[pivotRowPosition];
    01576
    01577       markRow[iRow] = l - lSave;
    01578       indexRowL[l] = iRow;
    01579       elementL[l] = element[pivotRowPosition];
    01580       l++;
    01581       //take out of row list
    01582       CoinBigIndex start = startRow[iRow];
    01583       CoinBigIndex end = start + numberInRow[iRow];
    01584       CoinBigIndex where = start;
    01585
    01586       while ( indexColumn[where] != iPivotColumn ) {
    01587         where++;
    01588       }                         /* endwhile */
    01589 #if DEBUG_COIN
    01590       if ( where >= end ) {
    01591         abort (  );
    01592       }
    01593 #endif
    01594       indexColumn[where] = indexColumn[end - 1];
    01595       numberInRow[iRow]--;
    01596       assert (numberInRow[iRow]>=0);
    01597     }
    01598     markRow[iPivotRow] = FAC_SET;
    01599     //compress pivot column (move pivot to front including saved)
    01600     numberInColumn[iPivotColumn] = 0;
    01601     //use end of L for temporary space
    01602     int *indexL = &indexRowL[lSave];
    01603     CoinFactorizationDouble *multipliersL = &elementL[lSave];
    01604
    01605     //adjust
    01606     int j;
    01607
    01608     for ( j = 0; j < numberDoColumn; j++ ) {
    01609       multipliersL[j] *= pivotMultiplier;
    01610     }
    01611     //zero out fill
    01612     CoinBigIndex iErase;
    01613     for ( iErase = 0; iErase < increment2 * numberDoRow;
    01614           iErase++ ) {
    01615       workArea2[iErase] = 0;
    01616     }
    01617     CoinBigIndex added = numberDoRow * numberDoColumn;
    01618     unsigned int *temp2 = workArea2;
    01619     int * nextColumn = nextColumn_.array();
    01620
    01621     //pack down and move to work
    01622     int jColumn;
    01623     for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01624       int iColumn = saveColumn[jColumn];
    01625       CoinBigIndex startColumnThis = startColumn[iColumn];
    01626       CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
    01627       int iRow = indexRow[startColumnThis];
    01628       CoinFactorizationDouble value = element[startColumnThis];
    01629       double largest;
    01630       CoinBigIndex put = startColumnThis;
    01631       CoinBigIndex positionLargest = -1;
    01632       CoinFactorizationDouble thisPivotValue = 0.0;
    01633
    01634       //compress column and find largest not updated
    01635       bool checkLargest;
    01636       int mark = markRow[iRow];
    01637
    01638       if ( mark == FAC_UNSET ) {
    01639         largest = fabs ( value );
    01640         positionLargest = put;
    01641         put++;
    01642         checkLargest = false;
    01643       } else {
    01644         //need to find largest
    01645         largest = 0.0;
    01646         checkLargest = true;
    01647         if ( mark != FAC_SET ) {
    01648           //will be updated
    01649           workArea[mark] = value;
    01650           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    01651           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    01652
    01653           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    01654           added--;
    01655         } else {
    01656           thisPivotValue = value;
    01657         }
    01658       }
    01659       CoinBigIndex i;
    01660       for ( i = startColumnThis + 1; i < endColumn; i++ ) {
    01661         iRow = indexRow[i];
    01662         value = element[i];
    01663         int mark = markRow[iRow];
    01664
    01665         if ( mark == FAC_UNSET ) {
    01666           //keep
    01667           indexRow[put] = iRow;
    01668           element[put] = value;
    01669           if ( checkLargest ) {
    01670             double absValue = fabs ( value );
    01671
    01672             if ( absValue > largest ) {
    01673               largest = absValue;
    01674               positionLargest = put;
    01675             }
    01676           }
    01677           put++;
    01678         } else if ( mark != FAC_SET ) {
    01679           //will be updated
    01680           workArea[mark] = value;
    01681           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    01682           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    01683
    01684           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    01685           added--;
    01686         } else {
    01687           thisPivotValue = value;
    01688         }
    01689       }
    01690       //slot in pivot
    01691       element[put] = element[startColumnThis];
    01692       indexRow[put] = indexRow[startColumnThis];
    01693       if ( positionLargest == startColumnThis ) {
    01694         positionLargest = put;  //follow if was largest
    01695       }
    01696       put++;
    01697       element[startColumnThis] = thisPivotValue;
    01698       indexRow[startColumnThis] = iPivotRow;
    01699       //clean up counts
    01700       startColumnThis++;
    01701       numberInColumn[iColumn] = put - startColumnThis;
    01702       int * numberInColumnPlus = numberInColumnPlus_.array();
    01703       numberInColumnPlus[iColumn]++;
    01704       startColumn[iColumn]++;
    01705       //how much space have we got
    01706       int next = nextColumn[iColumn];
    01707       CoinBigIndex space;
    01708
    01709       space = startColumn[next] - put - numberInColumnPlus[next];
    01710       //assume no zero elements
    01711       if ( numberDoColumn > space ) {
    01712         //getColumnSpace also moves fixed part
    01713         if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
    01714           goto BAD_PIVOT;
    01715         }
    01716         //redo starts
    01717         positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
    01718         startColumnThis = startColumn[iColumn];
    01719         put = startColumnThis + numberInColumn[iColumn];
    01720       }
    01721       double tolerance = zeroTolerance_;
    01722
    01723       int *nextCount = nextCount_.array();
    01724       for ( j = 0; j < numberDoColumn; j++ ) {
    01725         value = workArea[j] - thisPivotValue * multipliersL[j];
    01726         double absValue = fabs ( value );
    01727
    01728         if ( absValue > tolerance ) {
    01729           workArea[j] = 0.0;
    01730           element[put] = value;
    01731           indexRow[put] = indexL[j];
    01732           if ( absValue > largest ) {
    01733             largest = absValue;
    01734             positionLargest = put;
    01735           }
    01736           put++;
    01737         } else {
    01738           workArea[j] = 0.0;
    01739           added--;
    01740           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01741           int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01742
    01743           if ( temp2[word] & ( 1 << bit ) ) {
    01744             //take out of row list
    01745             iRow = indexL[j];
    01746             CoinBigIndex start = startRow[iRow];
    01747             CoinBigIndex end = start + numberInRow[iRow];
    01748             CoinBigIndex where = start;
    01749
    01750             while ( indexColumn[where] != iColumn ) {
    01751               where++;
    01752             }                   /* endwhile */
    01753 #if DEBUG_COIN
    01754             if ( where >= end ) {
    01755               abort (  );
    01756             }
    01757 #endif
    01758             indexColumn[where] = indexColumn[end - 1];
    01759             numberInRow[iRow]--;
    01760           } else {
    01761             //make sure won't be added
    01762             int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01763             int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01764
    01765             temp2[word] = temp2[word] | ( 1 << bit );   //say already in counts
    01766           }
    01767         }
    01768       }
    01769       numberInColumn[iColumn] = put - startColumnThis;
    01770       //move largest
    01771       if ( positionLargest >= 0 ) {
    01772         value = element[positionLargest];
    01773         iRow = indexRow[positionLargest];
    01774         element[positionLargest] = element[startColumnThis];
    01775         indexRow[positionLargest] = indexRow[startColumnThis];
    01776         element[startColumnThis] = value;
    01777         indexRow[startColumnThis] = iRow;
    01778       }
    01779       //linked list for column
    01780       if ( nextCount[iColumn + numberRows_] != -2 ) {
    01781         //modify linked list
    01782         deleteLink ( iColumn + numberRows_ );
    01783         addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
    01784       }
    01785       temp2 += increment2;
    01786     }
    01787     //get space for row list
    01788     unsigned int *putBase = workArea2;
    01789     int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
    01790     int i = 0;
    01791
    01792     // do linked lists and update counts
    01793     while ( bigLoops ) {
    01794       bigLoops--;
    01795       int bit;
    01796       for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
    01797         unsigned int *putThis = putBase;
    01798         int iRow = indexL[i];
    01799
    01800         //get space
    01801         int number = 0;
    01802         int jColumn;
    01803
    01804         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01805           unsigned int test = *putThis;
    01806
    01807           putThis += increment2;
    01808           test = 1 - ( ( test >> bit ) & 1 );
    01809           number += test;
    01810         }
    01811         int next = nextRow[iRow];
    01812         CoinBigIndex space;
    01813
    01814         space = startRow[next] - startRow[iRow];
    01815         number += numberInRow[iRow];
    01816         if ( space < number ) {
    01817           if ( !getRowSpace ( iRow, number ) ) {
    01818             goto BAD_PIVOT;
    01819           }
    01820         }
    01821         // now do
    01822         putThis = putBase;
    01823         next = nextRow[iRow];
    01824         number = numberInRow[iRow];
    01825         CoinBigIndex end = startRow[iRow] + number;
    01826         int saveIndex = indexColumn[startRow[next]];
    01827
    01828         //add in
    01829         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01830           unsigned int test = *putThis;
    01831
    01832           putThis += increment2;
    01833           test = 1 - ( ( test >> bit ) & 1 );
    01834           indexColumn[end] = saveColumn[jColumn];
    01835           end += test;
    01836         }
    01837         //put back next one in case zapped
    01838         indexColumn[startRow[next]] = saveIndex;
    01839         markRow[iRow] = FAC_UNSET;
    01840         number = end - startRow[iRow];
    01841         numberInRow[iRow] = number;
    01842         deleteLink ( iRow );
    01843         addLink ( iRow, number );
    01844       }
    01845       putBase++;
    01846     }                           /* endwhile */
    01847     int bit;
    01848
    01849     for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
    01850       unsigned int *putThis = putBase;
    01851       int iRow = indexL[i];
    01852
    01853       //get space
    01854       int number = 0;
    01855       int jColumn;
    01856
    01857       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01858         unsigned int test = *putThis;
    01859
    01860         putThis += increment2;
    01861         test = 1 - ( ( test >> bit ) & 1 );
    01862         number += test;
    01863       }
    01864       int next = nextRow[iRow];
    01865       CoinBigIndex space;
    01866
    01867       space = startRow[next] - startRow[iRow];
    01868       number += numberInRow[iRow];
    01869       if ( space < number ) {
    01870         if ( !getRowSpace ( iRow, number ) ) {
    01871           goto BAD_PIVOT;
    01872         }
    01873       }
    01874       // now do
    01875       putThis = putBase;
    01876       next = nextRow[iRow];
    01877       number = numberInRow[iRow];
    01878       CoinBigIndex end = startRow[iRow] + number;
    01879       int saveIndex;
    01880
    01881       saveIndex = indexColumn[startRow[next]];
    01882
    01883       //add in
    01884       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01885         unsigned int test = *putThis;
    01886
    01887         putThis += increment2;
    01888         test = 1 - ( ( test >> bit ) & 1 );
    01889
    01890         indexColumn[end] = saveColumn[jColumn];
    01891         end += test;
    01892       }
    01893       indexColumn[startRow[next]] = saveIndex;
    01894       markRow[iRow] = FAC_UNSET;
    01895       number = end - startRow[iRow];
    01896       numberInRow[iRow] = number;
    01897       deleteLink ( iRow );
    01898       addLink ( iRow, number );
    01899     }
    01900     markRow[iPivotRow] = FAC_UNSET;
    01901     //modify linked list for pivots
    01902     deleteLink ( iPivotRow );
    01903     deleteLink ( iPivotColumn + numberRows_ );
    01904     totalElements_ += added;
    01905     goodPivot= true;
    01906     // **** UGLY UGLY UGLY
    01907   }
    01908  BAD_PIVOT:
    01909
    01910   ;
    01911 }
    01912 #undef FAC_UNSET
    01913 #endif