- Main Page
- Namespaces
- Classes
- Files
- Directories

- Search for |

- File List
- File Members

# CoinFactorization.hpp

Go to the documentation of this file.

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