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

- File List
- File Members

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

Go to the documentation of this file.

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