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

- File List
- File Members

# /home/coin/SVN-release/Alps-1.2.1/CoinUtils/src/CoinFactorization.hpp

Go to the documentation of this file.

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