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

- File List
- File Members

# /home/coin/SVN-release/Bcps-0.91.0/CoinUtils/src/CoinFactorization.hpp

Go to the documentation of this file.

    00001 // Copyright (C) 2002, International Business Machines
    00002 // Corporation and others.  All Rights Reserved.
    00003
    00004 /*
    00005    Authors
    00006
    00007    John Forrest
    00008
    00009  */
    00010 #ifndef CoinFactorization_H
    00011 #define CoinFactorization_H
    00012 //#define COIN_ONE_ETA_COPY 100
    00013
    00014 #include <iostream>
    00015 #include <string>
    00016 #include <cassert>
    00017 #include <cstdio>
    00018 #include "CoinFinite.hpp"
    00019 #include "CoinIndexedVector.hpp"
    00020 class CoinPackedMatrix;
    00046 class CoinFactorization {
    00047    friend void CoinFactorizationUnitTest( const std::string & mpsDir );
    00048
    00049 public:
    00050
    00053
    00054     CoinFactorization (  );
    00056   CoinFactorization ( const CoinFactorization &other);
    00057
    00059    ~CoinFactorization (  );
    00061   void almostDestructor();
    00063   void show_self (  ) const;
    00065   int saveFactorization (const char * file  ) const;
    00069   int restoreFactorization (const char * file  , bool factor=false) ;
    00071   void sort (  ) const;
    00073     CoinFactorization & operator = ( const CoinFactorization & other );
    00075
    00085   int factorize ( const CoinPackedMatrix & matrix,
    00086                   int rowIsBasic[], int columnIsBasic[] ,
    00087                   double areaFactor = 0.0 );
    00098   int factorize ( int numberRows,
    00099                   int numberColumns,
    00100                   CoinBigIndex numberElements,
    00101                   CoinBigIndex maximumL,
    00102                   CoinBigIndex maximumU,
    00103                   const int indicesRow[],
    00104                   const int indicesColumn[], const double elements[] ,
    00105                   int permutation[],
    00106                   double areaFactor = 0.0);
    00111   int factorizePart1 ( int numberRows,
    00112                        int numberColumns,
    00113                        CoinBigIndex estimateNumberElements,
    00114                        int * indicesRow[],
    00115                        int * indicesColumn[],
    00116                        CoinFactorizationDouble * elements[],
    00117                   double areaFactor = 0.0);
    00124   int factorizePart2 (int permutation[],int exactNumberElements);
    00126   double conditionNumber() const;
    00127
    00129
    00132
    00133   inline int status (  ) const {
    00134     return status_;
    00135   }
    00137   inline void setStatus (  int value)
    00138   {  status_=value;  }
    00140   inline int pivots (  ) const {
    00141     return numberPivots_;
    00142   }
    00144   inline void setPivots (  int value )
    00145   { numberPivots_=value; }
    00147   inline int *permute (  ) const {
    00148     return permute_.array();
    00149   }
    00151   inline int *pivotColumn (  ) const {
    00152     return pivotColumn_.array();
    00153   }
    00155   inline CoinFactorizationDouble *pivotRegion (  ) const {
    00156     return pivotRegion_.array();
    00157   }
    00159   inline int *permuteBack (  ) const {
    00160     return permuteBack_.array();
    00161   }
    00164   inline int *pivotColumnBack (  ) const {
    00165     //return firstCount_.array();
    00166     return pivotColumnBack_.array();
    00167   }
    00169   inline CoinBigIndex * startRowL() const
    00170   { return startRowL_.array();}
    00171
    00173   inline CoinBigIndex * startColumnL() const
    00174   { return startColumnL_.array();}
    00175
    00177   inline int * indexColumnL() const
    00178   { return indexColumnL_.array();}
    00179
    00181   inline int * indexRowL() const
    00182   { return indexRowL_.array();}
    00183
    00185   inline CoinFactorizationDouble * elementByRowL() const
    00186   { return elementByRowL_.array();}
    00187
    00189   inline int numberRowsExtra (  ) const {
    00190     return numberRowsExtra_;
    00191   }
    00193   inline void setNumberRows(int value)
    00194   { numberRows_ = value; }
    00196   inline int numberRows (  ) const {
    00197     return numberRows_;
    00198   }
    00200   inline CoinBigIndex numberL() const
    00201   { return numberL_;}
    00202
    00204   inline CoinBigIndex baseL() const
    00205   { return baseL_;}
    00207   inline int maximumRowsExtra (  ) const {
    00208     return maximumRowsExtra_;
    00209   }
    00211   inline int numberColumns (  ) const {
    00212     return numberColumns_;
    00213   }
    00215   inline int numberElements (  ) const {
    00216     return totalElements_;
    00217   }
    00219   inline int numberForrestTomlin (  ) const {
    00220     return numberInColumn_.array()[numberColumnsExtra_];
    00221   }
    00223   inline int numberGoodColumns (  ) const {
    00224     return numberGoodU_;
    00225   }
    00227   inline double areaFactor (  ) const {
    00228     return areaFactor_;
    00229   }
    00230   inline void areaFactor ( double value ) {
    00231     areaFactor_=value;
    00232   }
    00234   double adjustedAreaFactor() const;
    00236   inline void relaxAccuracyCheck(double value)
    00237   { relaxCheck_ = value;}
    00238   inline double getAccuracyCheck() const
    00239   { return relaxCheck_;}
    00241   inline int messageLevel (  ) const {
    00242     return messageLevel_ ;
    00243   }
    00244   void messageLevel (  int value );
    00246   inline int maximumPivots (  ) const {
    00247     return maximumPivots_ ;
    00248   }
    00249   void maximumPivots (  int value );
    00250
    00252   inline int denseThreshold() const
    00253     { return denseThreshold_;}
    00255   inline void setDenseThreshold(int value)
    00256     { denseThreshold_ = value;}
    00258   inline double pivotTolerance (  ) const {
    00259     return pivotTolerance_ ;
    00260   }
    00261   void pivotTolerance (  double value );
    00263   inline double zeroTolerance (  ) const {
    00264     return zeroTolerance_ ;
    00265   }
    00266   void zeroTolerance (  double value );
    00267 #ifndef COIN_FAST_CODE
    00269   inline double slackValue (  ) const {
    00270     return slackValue_ ;
    00271   }
    00272   void slackValue (  double value );
    00273 #endif
    00275   double maximumCoefficient() const;
    00277   inline bool forrestTomlin() const
    00278   { return doForrestTomlin_;}
    00279   inline void setForrestTomlin(bool value)
    00280   { doForrestTomlin_=value;}
    00282   inline bool spaceForForrestTomlin() const
    00283   {
    00284     CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
    00285     CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
    00286     return (space>=0)&&doForrestTomlin_;
    00287   }
    00289
    00292
    00294   inline int numberDense() const
    00295   { return numberDense_;}
    00296
    00298   inline CoinBigIndex numberElementsU (  ) const {
    00299     return lengthU_;
    00300   }
    00302   inline void setNumberElementsU(CoinBigIndex value)
    00303   { lengthU_ = value; }
    00305   inline CoinBigIndex lengthAreaU (  ) const {
    00306     return lengthAreaU_;
    00307   }
    00309   inline CoinBigIndex numberElementsL (  ) const {
    00310     return lengthL_;
    00311   }
    00313   inline CoinBigIndex lengthAreaL (  ) const {
    00314     return lengthAreaL_;
    00315   }
    00317   inline CoinBigIndex numberElementsR (  ) const {
    00318     return lengthR_;
    00319   }
    00321   inline CoinBigIndex numberCompressions() const
    00322   { return numberCompressions_;}
    00324   inline int * numberInRow() const
    00325   { return numberInRow_.array();}
    00327   inline int * numberInColumn() const
    00328   { return numberInColumn_.array();}
    00330   inline CoinFactorizationDouble * elementU() const
    00331   { return elementU_.array();}
    00333   inline int * indexRowU() const
    00334   { return indexRowU_.array();}
    00336   inline CoinBigIndex * startColumnU() const
    00337   { return startColumnU_.array();}
    00339   inline int maximumColumnsExtra()
    00340   { return maximumColumnsExtra_;}
    00344   inline int biasLU() const
    00345   { return biasLU_;}
    00346   inline void setBiasLU(int value)
    00347   { biasLU_=value;}
    00353   inline int persistenceFlag() const
    00354   { return persistenceFlag_;}
    00355   void setPersistenceFlag(int value);
    00357
    00360
    00368   int replaceColumn ( CoinIndexedVector * regionSparse,
    00369                       int pivotRow,
    00370                       double pivotCheck ,
    00371                       bool checkBeforeModifying=false);
    00376   void replaceColumnU ( CoinIndexedVector * regionSparse,
    00377                         CoinBigIndex * deleted,
    00378                         int internalPivotRow);
    00380
    00390   int updateColumnFT ( CoinIndexedVector * regionSparse,
    00391                        CoinIndexedVector * regionSparse2);
    00394   int updateColumn ( CoinIndexedVector * regionSparse,
    00395                      CoinIndexedVector * regionSparse2,
    00396                      bool noPermute=false) const;
    00402   int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
    00403                            CoinIndexedVector * regionSparse2,
    00404                            CoinIndexedVector * regionSparse3,
    00405                            bool noPermuteRegion3=false) ;
    00410   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
    00411                               CoinIndexedVector * regionSparse2) const;
    00413   void goSparse();
    00415   inline int sparseThreshold ( ) const
    00416   { return sparseThreshold_;}
    00418   void sparseThreshold ( int value );
    00420
    00421
    00425
    00426   inline void clearArrays()
    00427   { gutsOfDestructor();}
    00429
    00432
    00435   int add ( CoinBigIndex numberElements,
    00436                int indicesRow[],
    00437                int indicesColumn[], double elements[] );
    00438
    00441   int addColumn ( CoinBigIndex numberElements,
    00442                      int indicesRow[], double elements[] );
    00443
    00446   int addRow ( CoinBigIndex numberElements,
    00447                   int indicesColumn[], double elements[] );
    00448
    00450   int deleteColumn ( int Row );
    00452   int deleteRow ( int Row );
    00453
    00457   int replaceRow ( int whichRow, int numberElements,
    00458                       const int indicesColumn[], const double elements[] );
    00460   void emptyRows(int numberToEmpty, const int which[]);
    00462
    00463
    00464   void checkSparse();
    00466   inline bool collectStatistics() const
    00467   { return collectStatistics_;}
    00469   inline void setCollectStatistics(bool onOff) const
    00470   { collectStatistics_ = onOff;}
    00472   void gutsOfDestructor(int type=1);
    00474   void gutsOfInitialize(int type);
    00475   void gutsOfCopy(const CoinFactorization &other);
    00476
    00478   void resetStatistics();
    00479
    00480
    00482
    00484
    00485   void getAreas ( int numberRows,
    00486                   int numberColumns,
    00487                   CoinBigIndex maximumL,
    00488                   CoinBigIndex maximumU );
    00489
    00492   void preProcess ( int state,
    00493                     int possibleDuplicates = -1 );
    00495   int factor (  );
    00496 protected:
    00499   int factorSparse (  );
    00502   int factorSparseSmall (  );
    00505   int factorSparseLarge (  );
    00508   int factorDense (  );
    00509
    00511   bool pivotOneOtherRow ( int pivotRow,
    00512                           int pivotColumn );
    00514   bool pivotRowSingleton ( int pivotRow,
    00515                            int pivotColumn );
    00517   bool pivotColumnSingleton ( int pivotRow,
    00518                               int pivotColumn );
    00519
    00524   bool getColumnSpace ( int iColumn,
    00525                         int extraNeeded );
    00526
    00529   bool reorderU();
    00533   bool getColumnSpaceIterateR ( int iColumn, double value,
    00534                                int iRow);
    00540   CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
    00541                                int iRow);
    00545   bool getRowSpace ( int iRow, int extraNeeded );
    00546
    00550   bool getRowSpaceIterate ( int iRow,
    00551                             int extraNeeded );
    00553   void checkConsistency (  );
    00555   inline void addLink ( int index, int count ) {
    00556     int *nextCount = nextCount_.array();
    00557     int *firstCount = firstCount_.array();
    00558     int *lastCount = lastCount_.array();
    00559     int next = firstCount[count];
    00560       lastCount[index] = -2 - count;
    00561     if ( next < 0 ) {
    00562       //first with that count
    00563       firstCount[count] = index;
    00564       nextCount[index] = -1;
    00565     } else {
    00566       firstCount[count] = index;
    00567       nextCount[index] = next;
    00568       lastCount[next] = index;
    00569   }}
    00571   inline void deleteLink ( int index ) {
    00572     int *nextCount = nextCount_.array();
    00573     int *firstCount = firstCount_.array();
    00574     int *lastCount = lastCount_.array();
    00575     int next = nextCount[index];
    00576     int last = lastCount[index];
    00577     if ( last >= 0 ) {
    00578       nextCount[last] = next;
    00579     } else {
    00580       int count = -last - 2;
    00581
    00582       firstCount[count] = next;
    00583     }
    00584     if ( next >= 0 ) {
    00585       lastCount[next] = last;
    00586     }
    00587     nextCount[index] = -2;
    00588     lastCount[index] = -2;
    00589     return;
    00590   }
    00592   void separateLinks(int count,bool rowsFirst);
    00594   void cleanup (  );
    00595
    00597   void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
    00599   void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
    00601   void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
    00603   void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
    00604
    00606   void updateColumnR ( CoinIndexedVector * region ) const;
    00609   void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
    00610
    00612   void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
    00613
    00615   void updateColumnUSparse ( CoinIndexedVector * regionSparse,
    00616                              int * indexIn) const;
    00618   void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
    00619                                int * indexIn) const;
    00621   int updateColumnUDensish ( double * COIN_RESTRICT region,
    00622                              int * COIN_RESTRICT regionIndex) const;
    00624   void updateTwoColumnsUDensish (
    00625                                  int & numberNonZero1,
    00626                                  double * COIN_RESTRICT region1,
    00627                                  int * COIN_RESTRICT index1,
    00628                                  int & numberNonZero2,
    00629                                  double * COIN_RESTRICT region2,
    00630                                  int * COIN_RESTRICT index2) const;
    00632   void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
    00634   void permuteBack ( CoinIndexedVector * regionSparse,
    00635                      CoinIndexedVector * outVector) const;
    00636
    00638   void updateColumnTransposePFI ( CoinIndexedVector * region) const;
    00641   void updateColumnTransposeU ( CoinIndexedVector * region,
    00642                                 int smallestIndex) const;
    00645   void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
    00646                                         int smallestIndex) const;
    00649   void updateColumnTransposeUDensish ( CoinIndexedVector * region,
    00650                                        int smallestIndex) const;
    00653   void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
    00656   void updateColumnTransposeUByColumn ( CoinIndexedVector * region,
    00657                                         int smallestIndex) const;
    00658
    00660   void updateColumnTransposeR ( CoinIndexedVector * region ) const;
    00662   void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
    00664   void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
    00665
    00667   void updateColumnTransposeL ( CoinIndexedVector * region ) const;
    00669   void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
    00671   void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
    00673   void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
    00675   void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
    00676 public:
    00681   int replaceColumnPFI ( CoinIndexedVector * regionSparse,
    00682                          int pivotRow, double alpha);
    00683 protected:
    00686   int checkPivot(double saveFromU, double oldPivot) const;
    00687   /********************************* START LARGE TEMPLATE ********/
    00688 #ifdef INT_IS_8
    00689 #define COINFACTORIZATION_BITS_PER_INT 64
    00690 #define COINFACTORIZATION_SHIFT_PER_INT 6
    00691 #define COINFACTORIZATION_MASK_PER_INT 0x3f
    00692 #else
    00693 #define COINFACTORIZATION_BITS_PER_INT 32
    00694 #define COINFACTORIZATION_SHIFT_PER_INT 5
    00695 #define COINFACTORIZATION_MASK_PER_INT 0x1f
    00696 #endif
    00697   template <class T>  inline bool
    00698   pivot ( int pivotRow,
    00699           int pivotColumn,
    00700           CoinBigIndex pivotRowPosition,
    00701           CoinBigIndex pivotColumnPosition,
    00702           CoinFactorizationDouble work[],
    00703           unsigned int workArea2[],
    00704           int increment,
    00705           int increment2,
    00706           T markRow[] ,
    00707           int largeInteger)
    00708 {
    00709   int *indexColumnU = indexColumnU_.array();
    00710   CoinBigIndex *startColumnU = startColumnU_.array();
    00711   int *numberInColumn = numberInColumn_.array();
    00712   CoinFactorizationDouble *elementU = elementU_.array();
    00713   int *indexRowU = indexRowU_.array();
    00714   CoinBigIndex *startRowU = startRowU_.array();
    00715   int *numberInRow = numberInRow_.array();
    00716   CoinFactorizationDouble *elementL = elementL_.array();
    00717   int *indexRowL = indexRowL_.array();
    00718   int *saveColumn = saveColumn_.array();
    00719   int *nextRow = nextRow_.array();
    00720   int *lastRow = lastRow_.array() ;
    00721
    00722   //store pivot columns (so can easily compress)
    00723   int numberInPivotRow = numberInRow[pivotRow] - 1;
    00724   CoinBigIndex startColumn = startColumnU[pivotColumn];
    00725   int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
    00726   CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
    00727   int put = 0;
    00728   CoinBigIndex startRow = startRowU[pivotRow];
    00729   CoinBigIndex endRow = startRow + numberInPivotRow + 1;
    00730
    00731   if ( pivotColumnPosition < 0 ) {
    00732     for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    00733       int iColumn = indexColumnU[pivotColumnPosition];
    00734       if ( iColumn != pivotColumn ) {
    00735         saveColumn[put++] = iColumn;
    00736       } else {
    00737         break;
    00738       }
    00739     }
    00740   } else {
    00741     for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
    00742       saveColumn[put++] = indexColumnU[i];
    00743     }
    00744   }
    00745   assert (pivotColumnPosition<endRow);
    00746   assert (indexColumnU[pivotColumnPosition]==pivotColumn);
    00747   pivotColumnPosition++;
    00748   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    00749     saveColumn[put++] = indexColumnU[pivotColumnPosition];
    00750   }
    00751   //take out this bit of indexColumnU
    00752   int next = nextRow[pivotRow];
    00753   int last = lastRow[pivotRow];
    00754
    00755   nextRow[last] = next;
    00756   lastRow[next] = last;
    00757   nextRow[pivotRow] = numberGoodU_;     //use for permute
    00758   lastRow[pivotRow] = -2;
    00759   numberInRow[pivotRow] = 0;
    00760   //store column in L, compress in U and take column out
    00761   CoinBigIndex l = lengthL_;
    00762
    00763   if ( l + numberInPivotColumn > lengthAreaL_ ) {
    00764     //need more memory
    00765     if ((messageLevel_&4)!=0)
    00766       printf("more memory needed in middle of invert\n");
    00767     return false;
    00768   }
    00769   //l+=currentAreaL_->elementByColumn-elementL;
    00770   CoinBigIndex lSave = l;
    00771
    00772   CoinBigIndex * startColumnL = startColumnL_.array();
    00773   startColumnL[numberGoodL_] = l;       //for luck and first time
    00774   numberGoodL_++;
    00775   startColumnL[numberGoodL_] = l + numberInPivotColumn;
    00776   lengthL_ += numberInPivotColumn;
    00777   if ( pivotRowPosition < 0 ) {
    00778     for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    00779       int iRow = indexRowU[pivotRowPosition];
    00780       if ( iRow != pivotRow ) {
    00781         indexRowL[l] = iRow;
    00782         elementL[l] = elementU[pivotRowPosition];
    00783         markRow[iRow] = static_cast<T>(l - lSave);
    00784         l++;
    00785         //take out of row list
    00786         CoinBigIndex start = startRowU[iRow];
    00787         CoinBigIndex end = start + numberInRow[iRow];
    00788         CoinBigIndex where = start;
    00789
    00790         while ( indexColumnU[where] != pivotColumn ) {
    00791           where++;
    00792         }                       /* endwhile */
    00793 #if DEBUG_COIN
    00794         if ( where >= end ) {
    00795           abort (  );
    00796         }
    00797 #endif
    00798         indexColumnU[where] = indexColumnU[end - 1];
    00799         numberInRow[iRow]--;
    00800       } else {
    00801         break;
    00802       }
    00803     }
    00804   } else {
    00805     CoinBigIndex i;
    00806
    00807     for ( i = startColumn; i < pivotRowPosition; i++ ) {
    00808       int iRow = indexRowU[i];
    00809
    00810       markRow[iRow] = static_cast<T>(l - lSave);
    00811       indexRowL[l] = iRow;
    00812       elementL[l] = elementU[i];
    00813       l++;
    00814       //take out of row list
    00815       CoinBigIndex start = startRowU[iRow];
    00816       CoinBigIndex end = start + numberInRow[iRow];
    00817       CoinBigIndex where = start;
    00818
    00819       while ( indexColumnU[where] != pivotColumn ) {
    00820         where++;
    00821       }                         /* endwhile */
    00822 #if DEBUG_COIN
    00823       if ( where >= end ) {
    00824         abort (  );
    00825       }
    00826 #endif
    00827       indexColumnU[where] = indexColumnU[end - 1];
    00828       numberInRow[iRow]--;
    00829       assert (numberInRow[iRow]>=0);
    00830     }
    00831   }
    00832   assert (pivotRowPosition<endColumn);
    00833   assert (indexRowU[pivotRowPosition]==pivotRow);
    00834   CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
    00835   CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
    00836
    00837   pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
    00838   pivotRowPosition++;
    00839   for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    00840     int iRow = indexRowU[pivotRowPosition];
    00841
    00842     markRow[iRow] = static_cast<T>(l - lSave);
    00843     indexRowL[l] = iRow;
    00844     elementL[l] = elementU[pivotRowPosition];
    00845     l++;
    00846     //take out of row list
    00847     CoinBigIndex start = startRowU[iRow];
    00848     CoinBigIndex end = start + numberInRow[iRow];
    00849     CoinBigIndex where = start;
    00850
    00851     while ( indexColumnU[where] != pivotColumn ) {
    00852       where++;
    00853     }                           /* endwhile */
    00854 #if DEBUG_COIN
    00855     if ( where >= end ) {
    00856       abort (  );
    00857     }
    00858 #endif
    00859     indexColumnU[where] = indexColumnU[end - 1];
    00860     numberInRow[iRow]--;
    00861     assert (numberInRow[iRow]>=0);
    00862   }
    00863   markRow[pivotRow] = static_cast<T>(largeInteger);
    00864   //compress pivot column (move pivot to front including saved)
    00865   numberInColumn[pivotColumn] = 0;
    00866   //use end of L for temporary space
    00867   int *indexL = &indexRowL[lSave];
    00868   CoinFactorizationDouble *multipliersL = &elementL[lSave];
    00869
    00870   //adjust
    00871   int j;
    00872
    00873   for ( j = 0; j < numberInPivotColumn; j++ ) {
    00874     multipliersL[j] *= pivotMultiplier;
    00875   }
    00876   //zero out fill
    00877   CoinBigIndex iErase;
    00878   for ( iErase = 0; iErase < increment2 * numberInPivotRow;
    00879         iErase++ ) {
    00880     workArea2[iErase] = 0;
    00881   }
    00882   CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
    00883   unsigned int *temp2 = workArea2;
    00884   int * nextColumn = nextColumn_.array();
    00885
    00886   //pack down and move to work
    00887   int jColumn;
    00888   for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    00889     int iColumn = saveColumn[jColumn];
    00890     CoinBigIndex startColumn = startColumnU[iColumn];
    00891     CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
    00892     int iRow = indexRowU[startColumn];
    00893     CoinFactorizationDouble value = elementU[startColumn];
    00894     double largest;
    00895     CoinBigIndex put = startColumn;
    00896     CoinBigIndex positionLargest = -1;
    00897     CoinFactorizationDouble thisPivotValue = 0.0;
    00898
    00899     //compress column and find largest not updated
    00900     bool checkLargest;
    00901     int mark = markRow[iRow];
    00902
    00903     if ( mark == largeInteger+1 ) {
    00904       largest = fabs ( value );
    00905       positionLargest = put;
    00906       put++;
    00907       checkLargest = false;
    00908     } else {
    00909       //need to find largest
    00910       largest = 0.0;
    00911       checkLargest = true;
    00912       if ( mark != largeInteger ) {
    00913         //will be updated
    00914         work[mark] = value;
    00915         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    00916         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    00917
    00918         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
    00919         added--;
    00920       } else {
    00921         thisPivotValue = value;
    00922       }
    00923     }
    00924     CoinBigIndex i;
    00925     for ( i = startColumn + 1; i < endColumn; i++ ) {
    00926       iRow = indexRowU[i];
    00927       value = elementU[i];
    00928       int mark = markRow[iRow];
    00929
    00930       if ( mark == largeInteger+1 ) {
    00931         //keep
    00932         indexRowU[put] = iRow;
    00933         elementU[put] = value;
    00934         if ( checkLargest ) {
    00935           double absValue = fabs ( value );
    00936
    00937           if ( absValue > largest ) {
    00938             largest = absValue;
    00939             positionLargest = put;
    00940           }
    00941         }
    00942         put++;
    00943       } else if ( mark != largeInteger ) {
    00944         //will be updated
    00945         work[mark] = value;
    00946         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    00947         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    00948
    00949         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
    00950         added--;
    00951       } else {
    00952         thisPivotValue = value;
    00953       }
    00954     }
    00955     //slot in pivot
    00956     elementU[put] = elementU[startColumn];
    00957     indexRowU[put] = indexRowU[startColumn];
    00958     if ( positionLargest == startColumn ) {
    00959       positionLargest = put;    //follow if was largest
    00960     }
    00961     put++;
    00962     elementU[startColumn] = thisPivotValue;
    00963     indexRowU[startColumn] = pivotRow;
    00964     //clean up counts
    00965     startColumn++;
    00966     numberInColumn[iColumn] = put - startColumn;
    00967     int * numberInColumnPlus = numberInColumnPlus_.array();
    00968     numberInColumnPlus[iColumn]++;
    00969     startColumnU[iColumn]++;
    00970     //how much space have we got
    00971     int next = nextColumn[iColumn];
    00972     CoinBigIndex space;
    00973
    00974     space = startColumnU[next] - put - numberInColumnPlus[next];
    00975     //assume no zero elements
    00976     if ( numberInPivotColumn > space ) {
    00977       //getColumnSpace also moves fixed part
    00978       if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
    00979         return false;
    00980       }
    00981       //redo starts
    00982       positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
    00983       startColumn = startColumnU[iColumn];
    00984       put = startColumn + numberInColumn[iColumn];
    00985     }
    00986     double tolerance = zeroTolerance_;
    00987
    00988     int *nextCount = nextCount_.array();
    00989     for ( j = 0; j < numberInPivotColumn; j++ ) {
    00990       value = work[j] - thisPivotValue * multipliersL[j];
    00991       double absValue = fabs ( value );
    00992
    00993       if ( absValue > tolerance ) {
    00994         work[j] = 0.0;
    00995         elementU[put] = value;
    00996         indexRowU[put] = indexL[j];
    00997         if ( absValue > largest ) {
    00998           largest = absValue;
    00999           positionLargest = put;
    01000         }
    01001         put++;
    01002       } else {
    01003         work[j] = 0.0;
    01004         added--;
    01005         int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01006         int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01007
    01008         if ( temp2[word] & ( 1 << bit ) ) {
    01009           //take out of row list
    01010           iRow = indexL[j];
    01011           CoinBigIndex start = startRowU[iRow];
    01012           CoinBigIndex end = start + numberInRow[iRow];
    01013           CoinBigIndex where = start;
    01014
    01015           while ( indexColumnU[where] != iColumn ) {
    01016             where++;
    01017           }                     /* endwhile */
    01018 #if DEBUG_COIN
    01019           if ( where >= end ) {
    01020             abort (  );
    01021           }
    01022 #endif
    01023           indexColumnU[where] = indexColumnU[end - 1];
    01024           numberInRow[iRow]--;
    01025         } else {
    01026           //make sure won't be added
    01027           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01028           int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01029
    01030           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    01031         }
    01032       }
    01033     }
    01034     numberInColumn[iColumn] = put - startColumn;
    01035     //move largest
    01036     if ( positionLargest >= 0 ) {
    01037       value = elementU[positionLargest];
    01038       iRow = indexRowU[positionLargest];
    01039       elementU[positionLargest] = elementU[startColumn];
    01040       indexRowU[positionLargest] = indexRowU[startColumn];
    01041       elementU[startColumn] = value;
    01042       indexRowU[startColumn] = iRow;
    01043     }
    01044     //linked list for column
    01045     if ( nextCount[iColumn + numberRows_] != -2 ) {
    01046       //modify linked list
    01047       deleteLink ( iColumn + numberRows_ );
    01048       addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
    01049     }
    01050     temp2 += increment2;
    01051   }
    01052   //get space for row list
    01053   unsigned int *putBase = workArea2;
    01054   int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
    01055   int i = 0;
    01056
    01057   // do linked lists and update counts
    01058   while ( bigLoops ) {
    01059     bigLoops--;
    01060     int bit;
    01061     for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
    01062       unsigned int *putThis = putBase;
    01063       int iRow = indexL[i];
    01064
    01065       //get space
    01066       int number = 0;
    01067       int jColumn;
    01068
    01069       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01070         unsigned int test = *putThis;
    01071
    01072         putThis += increment2;
    01073         test = 1 - ( ( test >> bit ) & 1 );
    01074         number += test;
    01075       }
    01076       int next = nextRow[iRow];
    01077       CoinBigIndex space;
    01078
    01079       space = startRowU[next] - startRowU[iRow];
    01080       number += numberInRow[iRow];
    01081       if ( space < number ) {
    01082         if ( !getRowSpace ( iRow, number ) ) {
    01083           return false;
    01084         }
    01085       }
    01086       // now do
    01087       putThis = putBase;
    01088       next = nextRow[iRow];
    01089       number = numberInRow[iRow];
    01090       CoinBigIndex end = startRowU[iRow] + number;
    01091       int saveIndex = indexColumnU[startRowU[next]];
    01092
    01093       //add in
    01094       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01095         unsigned int test = *putThis;
    01096
    01097         putThis += increment2;
    01098         test = 1 - ( ( test >> bit ) & 1 );
    01099         indexColumnU[end] = saveColumn[jColumn];
    01100         end += test;
    01101       }
    01102       //put back next one in case zapped
    01103       indexColumnU[startRowU[next]] = saveIndex;
    01104       markRow[iRow] = static_cast<T>(largeInteger+1);
    01105       number = end - startRowU[iRow];
    01106       numberInRow[iRow] = number;
    01107       deleteLink ( iRow );
    01108       addLink ( iRow, number );
    01109     }
    01110     putBase++;
    01111   }                             /* endwhile */
    01112   int bit;
    01113
    01114   for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
    01115     unsigned int *putThis = putBase;
    01116     int iRow = indexL[i];
    01117
    01118     //get space
    01119     int number = 0;
    01120     int jColumn;
    01121
    01122     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01123       unsigned int test = *putThis;
    01124
    01125       putThis += increment2;
    01126       test = 1 - ( ( test >> bit ) & 1 );
    01127       number += test;
    01128     }
    01129     int next = nextRow[iRow];
    01130     CoinBigIndex space;
    01131
    01132     space = startRowU[next] - startRowU[iRow];
    01133     number += numberInRow[iRow];
    01134     if ( space < number ) {
    01135       if ( !getRowSpace ( iRow, number ) ) {
    01136         return false;
    01137       }
    01138     }
    01139     // now do
    01140     putThis = putBase;
    01141     next = nextRow[iRow];
    01142     number = numberInRow[iRow];
    01143     CoinBigIndex end = startRowU[iRow] + number;
    01144     int saveIndex;
    01145
    01146     saveIndex = indexColumnU[startRowU[next]];
    01147
    01148     //add in
    01149     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01150       unsigned int test = *putThis;
    01151
    01152       putThis += increment2;
    01153       test = 1 - ( ( test >> bit ) & 1 );
    01154
    01155       indexColumnU[end] = saveColumn[jColumn];
    01156       end += test;
    01157     }
    01158     indexColumnU[startRowU[next]] = saveIndex;
    01159     markRow[iRow] = static_cast<T>(largeInteger+1);
    01160     number = end - startRowU[iRow];
    01161     numberInRow[iRow] = number;
    01162     deleteLink ( iRow );
    01163     addLink ( iRow, number );
    01164   }
    01165   markRow[pivotRow] = static_cast<T>(largeInteger+1);
    01166   //modify linked list for pivots
    01167   deleteLink ( pivotRow );
    01168   deleteLink ( pivotColumn + numberRows_ );
    01169   totalElements_ += added;
    01170   return true;
    01171 }
    01172
    01173   /********************************* END LARGE TEMPLATE ********/
    01175
    01176 protected:
    01177
    01180
    01181   double pivotTolerance_;
    01183   double zeroTolerance_;
    01184 #ifndef COIN_FAST_CODE
    01186   double slackValue_;
    01187 #else
    01188 #ifndef slackValue_
    01189 #define slackValue_ -1.0
    01190 #endif
    01191 #endif
    01193   double areaFactor_;
    01195   double relaxCheck_;
    01197   int numberRows_;
    01199   int numberRowsExtra_;
    01201   int maximumRowsExtra_;
    01203   int numberColumns_;
    01205   int numberColumnsExtra_;
    01207   int maximumColumnsExtra_;
    01209   int numberGoodU_;
    01211   int numberGoodL_;
    01213   int maximumPivots_;
    01215   int numberPivots_;
    01218   CoinBigIndex totalElements_;
    01220   CoinBigIndex factorElements_;
    01222   CoinIntArrayWithLength pivotColumn_;
    01224   CoinIntArrayWithLength permute_;
    01226   CoinIntArrayWithLength permuteBack_;
    01228   CoinIntArrayWithLength pivotColumnBack_;
    01230   int status_;
    01231
    01236   //int increasingRows_;
    01237
    01239   int numberTrials_;
    01241   CoinBigIndexArrayWithLength startRowU_;
    01242
    01244   CoinIntArrayWithLength numberInRow_;
    01245
    01247   CoinIntArrayWithLength numberInColumn_;
    01248
    01250   CoinIntArrayWithLength numberInColumnPlus_;
    01251
    01254   CoinIntArrayWithLength firstCount_;
    01255
    01257   CoinIntArrayWithLength nextCount_;
    01258
    01260   CoinIntArrayWithLength lastCount_;
    01261
    01263   CoinIntArrayWithLength nextColumn_;
    01264
    01266   CoinIntArrayWithLength lastColumn_;
    01267
    01269   CoinIntArrayWithLength nextRow_;
    01270
    01272   CoinIntArrayWithLength lastRow_;
    01273
    01275   CoinIntArrayWithLength saveColumn_;
    01276
    01278   CoinIntArrayWithLength markRow_;
    01279
    01281   int messageLevel_;
    01282
    01284   int biggerDimension_;
    01285
    01287   CoinIntArrayWithLength indexColumnU_;
    01288
    01290   CoinIntArrayWithLength pivotRowL_;
    01291
    01293   CoinFactorizationDoubleArrayWithLength pivotRegion_;
    01294
    01296   int numberSlacks_;
    01297
    01299   int numberU_;
    01300
    01302   CoinBigIndex maximumU_;
    01303
    01305   //int baseU_;
    01306
    01308   CoinBigIndex lengthU_;
    01309
    01311   CoinBigIndex lengthAreaU_;
    01312
    01314   CoinFactorizationDoubleArrayWithLength elementU_;
    01315
    01317   CoinIntArrayWithLength indexRowU_;
    01318
    01320   CoinBigIndexArrayWithLength startColumnU_;
    01321
    01323   CoinBigIndexArrayWithLength convertRowToColumnU_;
    01324
    01326   CoinBigIndex numberL_;
    01327
    01329   CoinBigIndex baseL_;
    01330
    01332   CoinBigIndex lengthL_;
    01333
    01335   CoinBigIndex lengthAreaL_;
    01336
    01338   CoinFactorizationDoubleArrayWithLength elementL_;
    01339
    01341   CoinIntArrayWithLength indexRowL_;
    01342
    01344   CoinBigIndexArrayWithLength startColumnL_;
    01345
    01347   bool doForrestTomlin_;
    01348
    01350   int numberR_;
    01351
    01353   CoinBigIndex lengthR_;
    01354
    01356   CoinBigIndex lengthAreaR_;
    01357
    01359   CoinFactorizationDouble *elementR_;
    01360
    01362   int *indexRowR_;
    01363
    01365   CoinBigIndexArrayWithLength startColumnR_;
    01366
    01368   double  * denseArea_;
    01369
    01371   int * densePermute_;
    01372
    01374   int numberDense_;
    01375
    01377   int denseThreshold_;
    01378
    01380   CoinFactorizationDoubleArrayWithLength workArea_;
    01381
    01383   CoinUnsignedIntArrayWithLength workArea2_;
    01384
    01386   CoinBigIndex numberCompressions_;
    01387
    01389   mutable double ftranCountInput_;
    01390   mutable double ftranCountAfterL_;
    01391   mutable double ftranCountAfterR_;
    01392   mutable double ftranCountAfterU_;
    01393   mutable double btranCountInput_;
    01394   mutable double btranCountAfterU_;
    01395   mutable double btranCountAfterR_;
    01396   mutable double btranCountAfterL_;
    01397
    01399   mutable int numberFtranCounts_;
    01400   mutable int numberBtranCounts_;
    01401
    01403   double ftranAverageAfterL_;
    01404   double ftranAverageAfterR_;
    01405   double ftranAverageAfterU_;
    01406   double btranAverageAfterU_;
    01407   double btranAverageAfterR_;
    01408   double btranAverageAfterL_;
    01409
    01411   mutable bool collectStatistics_;
    01412
    01414   int sparseThreshold_;
    01415
    01417   int sparseThreshold2_;
    01418
    01420   CoinBigIndexArrayWithLength startRowL_;
    01421
    01423   CoinIntArrayWithLength indexColumnL_;
    01424
    01426   CoinFactorizationDoubleArrayWithLength elementByRowL_;
    01427
    01429   mutable CoinIntArrayWithLength sparse_;
    01433   int biasLU_;
    01439   int persistenceFlag_;
    01441 };
    01442 // Dense coding
    01443 #ifdef COIN_HAS_LAPACK
    01444 #define DENSE_CODE 1
    01445 /* Type of Fortran integer translated into C */
    01446 #ifndef ipfint
    01447 //typedef ipfint FORTRAN_INTEGER_TYPE ;
    01448 typedef int ipfint;
    01449 typedef const int cipfint;
    01450 #endif
    01451 #endif
    01452 #endif
    01453 // Extra for ugly include
    01454 #ifdef UGLY_COIN_FACTOR_CODING
    01455 #define FAC_UNSET (FAC_SET+1)
    01456 {
    01457   goodPivot=false;
    01458   //store pivot columns (so can easily compress)
    01459   CoinBigIndex startColumnThis = startColumn[iPivotColumn];
    01460   CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
    01461   int put = 0;
    01462   CoinBigIndex startRowThis = startRow[iPivotRow];
    01463   CoinBigIndex endRow = startRowThis + numberDoRow + 1;
    01464   if ( pivotColumnPosition < 0 ) {
    01465     for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    01466       int iColumn = indexColumn[pivotColumnPosition];
    01467       if ( iColumn != iPivotColumn ) {
    01468         saveColumn[put++] = iColumn;
    01469       } else {
    01470         break;
    01471       }
    01472     }
    01473   } else {
    01474     for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
    01475       saveColumn[put++] = indexColumn[i];
    01476     }
    01477   }
    01478   assert (pivotColumnPosition<endRow);
    01479   assert (indexColumn[pivotColumnPosition]==iPivotColumn);
    01480   pivotColumnPosition++;
    01481   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    01482     saveColumn[put++] = indexColumn[pivotColumnPosition];
    01483   }
    01484   //take out this bit of indexColumn
    01485   int next = nextRow[iPivotRow];
    01486   int last = lastRow[iPivotRow];
    01487
    01488   nextRow[last] = next;
    01489   lastRow[next] = last;
    01490   nextRow[iPivotRow] = numberGoodU_;    //use for permute
    01491   lastRow[iPivotRow] = -2;
    01492   numberInRow[iPivotRow] = 0;
    01493   //store column in L, compress in U and take column out
    01494   CoinBigIndex l = lengthL_;
    01495   // **** HORRID coding coming up but a goto seems best!
    01496   {
    01497     if ( l + numberDoColumn > lengthAreaL_ ) {
    01498       //need more memory
    01499       if ((messageLevel_&4)!=0)
    01500         printf("more memory needed in middle of invert\n");
    01501       goto BAD_PIVOT;
    01502     }
    01503     //l+=currentAreaL_->elementByColumn-elementL;
    01504     CoinBigIndex lSave = l;
    01505
    01506     CoinBigIndex * startColumnL = startColumnL_.array();
    01507     startColumnL[numberGoodL_] = l;     //for luck and first time
    01508     numberGoodL_++;
    01509     startColumnL[numberGoodL_] = l + numberDoColumn;
    01510     lengthL_ += numberDoColumn;
    01511     if ( pivotRowPosition < 0 ) {
    01512       for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    01513         int iRow = indexRow[pivotRowPosition];
    01514         if ( iRow != iPivotRow ) {
    01515           indexRowL[l] = iRow;
    01516           elementL[l] = element[pivotRowPosition];
    01517           markRow[iRow] = l - lSave;
    01518           l++;
    01519           //take out of row list
    01520           CoinBigIndex start = startRow[iRow];
    01521           CoinBigIndex end = start + numberInRow[iRow];
    01522           CoinBigIndex where = start;
    01523
    01524           while ( indexColumn[where] != iPivotColumn ) {
    01525             where++;
    01526           }                     /* endwhile */
    01527 #if DEBUG_COIN
    01528           if ( where >= end ) {
    01529             abort (  );
    01530           }
    01531 #endif
    01532           indexColumn[where] = indexColumn[end - 1];
    01533           numberInRow[iRow]--;
    01534         } else {
    01535           break;
    01536         }
    01537       }
    01538     } else {
    01539       CoinBigIndex i;
    01540
    01541       for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
    01542         int iRow = indexRow[i];
    01543
    01544         markRow[iRow] = l - lSave;
    01545         indexRowL[l] = iRow;
    01546         elementL[l] = element[i];
    01547         l++;
    01548         //take out of row list
    01549         CoinBigIndex start = startRow[iRow];
    01550         CoinBigIndex end = start + numberInRow[iRow];
    01551         CoinBigIndex where = start;
    01552
    01553         while ( indexColumn[where] != iPivotColumn ) {
    01554           where++;
    01555         }                               /* endwhile */
    01556 #if DEBUG_COIN
    01557         if ( where >= end ) {
    01558           abort (  );
    01559         }
    01560 #endif
    01561         indexColumn[where] = indexColumn[end - 1];
    01562         numberInRow[iRow]--;
    01563         assert (numberInRow[iRow]>=0);
    01564       }
    01565     }
    01566     assert (pivotRowPosition<endColumn);
    01567     assert (indexRow[pivotRowPosition]==iPivotRow);
    01568     CoinFactorizationDouble pivotElement = element[pivotRowPosition];
    01569     CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
    01570
    01571     pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
    01572     pivotRowPosition++;
    01573     for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    01574       int iRow = indexRow[pivotRowPosition];
    01575
    01576       markRow[iRow] = l - lSave;
    01577       indexRowL[l] = iRow;
    01578       elementL[l] = element[pivotRowPosition];
    01579       l++;
    01580       //take out of row list
    01581       CoinBigIndex start = startRow[iRow];
    01582       CoinBigIndex end = start + numberInRow[iRow];
    01583       CoinBigIndex where = start;
    01584
    01585       while ( indexColumn[where] != iPivotColumn ) {
    01586         where++;
    01587       }                         /* endwhile */
    01588 #if DEBUG_COIN
    01589       if ( where >= end ) {
    01590         abort (  );
    01591       }
    01592 #endif
    01593       indexColumn[where] = indexColumn[end - 1];
    01594       numberInRow[iRow]--;
    01595       assert (numberInRow[iRow]>=0);
    01596     }
    01597     markRow[iPivotRow] = FAC_SET;
    01598     //compress pivot column (move pivot to front including saved)
    01599     numberInColumn[iPivotColumn] = 0;
    01600     //use end of L for temporary space
    01601     int *indexL = &indexRowL[lSave];
    01602     CoinFactorizationDouble *multipliersL = &elementL[lSave];
    01603
    01604     //adjust
    01605     int j;
    01606
    01607     for ( j = 0; j < numberDoColumn; j++ ) {
    01608       multipliersL[j] *= pivotMultiplier;
    01609     }
    01610     //zero out fill
    01611     CoinBigIndex iErase;
    01612     for ( iErase = 0; iErase < increment2 * numberDoRow;
    01613           iErase++ ) {
    01614       workArea2[iErase] = 0;
    01615     }
    01616     CoinBigIndex added = numberDoRow * numberDoColumn;
    01617     unsigned int *temp2 = workArea2;
    01618     int * nextColumn = nextColumn_.array();
    01619
    01620     //pack down and move to work
    01621     int jColumn;
    01622     for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01623       int iColumn = saveColumn[jColumn];
    01624       CoinBigIndex startColumnThis = startColumn[iColumn];
    01625       CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
    01626       int iRow = indexRow[startColumnThis];
    01627       CoinFactorizationDouble value = element[startColumnThis];
    01628       double largest;
    01629       CoinBigIndex put = startColumnThis;
    01630       CoinBigIndex positionLargest = -1;
    01631       CoinFactorizationDouble thisPivotValue = 0.0;
    01632
    01633       //compress column and find largest not updated
    01634       bool checkLargest;
    01635       int mark = markRow[iRow];
    01636
    01637       if ( mark == FAC_UNSET ) {
    01638         largest = fabs ( value );
    01639         positionLargest = put;
    01640         put++;
    01641         checkLargest = false;
    01642       } else {
    01643         //need to find largest
    01644         largest = 0.0;
    01645         checkLargest = true;
    01646         if ( mark != FAC_SET ) {
    01647           //will be updated
    01648           workArea[mark] = value;
    01649           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    01650           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    01651
    01652           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    01653           added--;
    01654         } else {
    01655           thisPivotValue = value;
    01656         }
    01657       }
    01658       CoinBigIndex i;
    01659       for ( i = startColumnThis + 1; i < endColumn; i++ ) {
    01660         iRow = indexRow[i];
    01661         value = element[i];
    01662         int mark = markRow[iRow];
    01663
    01664         if ( mark == FAC_UNSET ) {
    01665           //keep
    01666           indexRow[put] = iRow;
    01667           element[put] = value;
    01668           if ( checkLargest ) {
    01669             double absValue = fabs ( value );
    01670
    01671             if ( absValue > largest ) {
    01672               largest = absValue;
    01673               positionLargest = put;
    01674             }
    01675           }
    01676           put++;
    01677         } else if ( mark != FAC_SET ) {
    01678           //will be updated
    01679           workArea[mark] = value;
    01680           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    01681           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    01682
    01683           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    01684           added--;
    01685         } else {
    01686           thisPivotValue = value;
    01687         }
    01688       }
    01689       //slot in pivot
    01690       element[put] = element[startColumnThis];
    01691       indexRow[put] = indexRow[startColumnThis];
    01692       if ( positionLargest == startColumnThis ) {
    01693         positionLargest = put;  //follow if was largest
    01694       }
    01695       put++;
    01696       element[startColumnThis] = thisPivotValue;
    01697       indexRow[startColumnThis] = iPivotRow;
    01698       //clean up counts
    01699       startColumnThis++;
    01700       numberInColumn[iColumn] = put - startColumnThis;
    01701       int * numberInColumnPlus = numberInColumnPlus_.array();
    01702       numberInColumnPlus[iColumn]++;
    01703       startColumn[iColumn]++;
    01704       //how much space have we got
    01705       int next = nextColumn[iColumn];
    01706       CoinBigIndex space;
    01707
    01708       space = startColumn[next] - put - numberInColumnPlus[next];
    01709       //assume no zero elements
    01710       if ( numberDoColumn > space ) {
    01711         //getColumnSpace also moves fixed part
    01712         if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
    01713           goto BAD_PIVOT;
    01714         }
    01715         //redo starts
    01716         positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
    01717         startColumnThis = startColumn[iColumn];
    01718         put = startColumnThis + numberInColumn[iColumn];
    01719       }
    01720       double tolerance = zeroTolerance_;
    01721
    01722       int *nextCount = nextCount_.array();
    01723       for ( j = 0; j < numberDoColumn; j++ ) {
    01724         value = workArea[j] - thisPivotValue * multipliersL[j];
    01725         double absValue = fabs ( value );
    01726
    01727         if ( absValue > tolerance ) {
    01728           workArea[j] = 0.0;
    01729           element[put] = value;
    01730           indexRow[put] = indexL[j];
    01731           if ( absValue > largest ) {
    01732             largest = absValue;
    01733             positionLargest = put;
    01734           }
    01735           put++;
    01736         } else {
    01737           workArea[j] = 0.0;
    01738           added--;
    01739           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01740           int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01741
    01742           if ( temp2[word] & ( 1 << bit ) ) {
    01743             //take out of row list
    01744             iRow = indexL[j];
    01745             CoinBigIndex start = startRow[iRow];
    01746             CoinBigIndex end = start + numberInRow[iRow];
    01747             CoinBigIndex where = start;
    01748
    01749             while ( indexColumn[where] != iColumn ) {
    01750               where++;
    01751             }                   /* endwhile */
    01752 #if DEBUG_COIN
    01753             if ( where >= end ) {
    01754               abort (  );
    01755             }
    01756 #endif
    01757             indexColumn[where] = indexColumn[end - 1];
    01758             numberInRow[iRow]--;
    01759           } else {
    01760             //make sure won't be added
    01761             int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01762             int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01763
    01764             temp2[word] = temp2[word] | ( 1 << bit );   //say already in counts
    01765           }
    01766         }
    01767       }
    01768       numberInColumn[iColumn] = put - startColumnThis;
    01769       //move largest
    01770       if ( positionLargest >= 0 ) {
    01771         value = element[positionLargest];
    01772         iRow = indexRow[positionLargest];
    01773         element[positionLargest] = element[startColumnThis];
    01774         indexRow[positionLargest] = indexRow[startColumnThis];
    01775         element[startColumnThis] = value;
    01776         indexRow[startColumnThis] = iRow;
    01777       }
    01778       //linked list for column
    01779       if ( nextCount[iColumn + numberRows_] != -2 ) {
    01780         //modify linked list
    01781         deleteLink ( iColumn + numberRows_ );
    01782         addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
    01783       }
    01784       temp2 += increment2;
    01785     }
    01786     //get space for row list
    01787     unsigned int *putBase = workArea2;
    01788     int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
    01789     int i = 0;
    01790
    01791     // do linked lists and update counts
    01792     while ( bigLoops ) {
    01793       bigLoops--;
    01794       int bit;
    01795       for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
    01796         unsigned int *putThis = putBase;
    01797         int iRow = indexL[i];
    01798
    01799         //get space
    01800         int number = 0;
    01801         int jColumn;
    01802
    01803         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01804           unsigned int test = *putThis;
    01805
    01806           putThis += increment2;
    01807           test = 1 - ( ( test >> bit ) & 1 );
    01808           number += test;
    01809         }
    01810         int next = nextRow[iRow];
    01811         CoinBigIndex space;
    01812
    01813         space = startRow[next] - startRow[iRow];
    01814         number += numberInRow[iRow];
    01815         if ( space < number ) {
    01816           if ( !getRowSpace ( iRow, number ) ) {
    01817             goto BAD_PIVOT;
    01818           }
    01819         }
    01820         // now do
    01821         putThis = putBase;
    01822         next = nextRow[iRow];
    01823         number = numberInRow[iRow];
    01824         CoinBigIndex end = startRow[iRow] + number;
    01825         int saveIndex = indexColumn[startRow[next]];
    01826
    01827         //add in
    01828         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01829           unsigned int test = *putThis;
    01830
    01831           putThis += increment2;
    01832           test = 1 - ( ( test >> bit ) & 1 );
    01833           indexColumn[end] = saveColumn[jColumn];
    01834           end += test;
    01835         }
    01836         //put back next one in case zapped
    01837         indexColumn[startRow[next]] = saveIndex;
    01838         markRow[iRow] = FAC_UNSET;
    01839         number = end - startRow[iRow];
    01840         numberInRow[iRow] = number;
    01841         deleteLink ( iRow );
    01842         addLink ( iRow, number );
    01843       }
    01844       putBase++;
    01845     }                           /* endwhile */
    01846     int bit;
    01847
    01848     for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
    01849       unsigned int *putThis = putBase;
    01850       int iRow = indexL[i];
    01851
    01852       //get space
    01853       int number = 0;
    01854       int jColumn;
    01855
    01856       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01857         unsigned int test = *putThis;
    01858
    01859         putThis += increment2;
    01860         test = 1 - ( ( test >> bit ) & 1 );
    01861         number += test;
    01862       }
    01863       int next = nextRow[iRow];
    01864       CoinBigIndex space;
    01865
    01866       space = startRow[next] - startRow[iRow];
    01867       number += numberInRow[iRow];
    01868       if ( space < number ) {
    01869         if ( !getRowSpace ( iRow, number ) ) {
    01870           goto BAD_PIVOT;
    01871         }
    01872       }
    01873       // now do
    01874       putThis = putBase;
    01875       next = nextRow[iRow];
    01876       number = numberInRow[iRow];
    01877       CoinBigIndex end = startRow[iRow] + number;
    01878       int saveIndex;
    01879
    01880       saveIndex = indexColumn[startRow[next]];
    01881
    01882       //add in
    01883       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01884         unsigned int test = *putThis;
    01885
    01886         putThis += increment2;
    01887         test = 1 - ( ( test >> bit ) & 1 );
    01888
    01889         indexColumn[end] = saveColumn[jColumn];
    01890         end += test;
    01891       }
    01892       indexColumn[startRow[next]] = saveIndex;
    01893       markRow[iRow] = FAC_UNSET;
    01894       number = end - startRow[iRow];
    01895       numberInRow[iRow] = number;
    01896       deleteLink ( iRow );
    01897       addLink ( iRow, number );
    01898     }
    01899     markRow[iPivotRow] = FAC_UNSET;
    01900     //modify linked list for pivots
    01901     deleteLink ( iPivotRow );
    01902     deleteLink ( iPivotColumn + numberRows_ );
    01903     totalElements_ += added;
    01904     goodPivot= true;
    01905     // **** UGLY UGLY UGLY
    01906   }
    01907  BAD_PIVOT:
    01908
    01909   ;
    01910 }
    01911 #undef FAC_UNSET
    01912 #endif