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

- File List
- File Members

# /home/coin/svn-release/Bcps-0.93.11/CoinUtils/src/CoinFactorization.hpp

Go to the documentation of this file.

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