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

- File List
- File Members

# /home/coin/svn-release/Alps-1.5.0/CoinUtils/src/CoinFactorization.hpp

Go to the documentation of this file.

    00001 /* $Id: CoinFactorization.hpp 1767 2015-01-05 12:36:13Z forrest $ */
    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   }
    00167   inline int *lastRow (  ) const {
    00168     return lastRow_.array();
    00169   }
    00172   inline int *pivotColumnBack (  ) const {
    00173     //return firstCount_.array();
    00174     return pivotColumnBack_.array();
    00175   }
    00177   inline CoinBigIndex * startRowL() const
    00178   { return startRowL_.array();}
    00179
    00181   inline CoinBigIndex * startColumnL() const
    00182   { return startColumnL_.array();}
    00183
    00185   inline int * indexColumnL() const
    00186   { return indexColumnL_.array();}
    00187
    00189   inline int * indexRowL() const
    00190   { return indexRowL_.array();}
    00191
    00193   inline CoinFactorizationDouble * elementByRowL() const
    00194   { return elementByRowL_.array();}
    00195
    00197   inline int numberRowsExtra (  ) const {
    00198     return numberRowsExtra_;
    00199   }
    00201   inline void setNumberRows(int value)
    00202   { numberRows_ = value; }
    00204   inline int numberRows (  ) const {
    00205     return numberRows_;
    00206   }
    00208   inline CoinBigIndex numberL() const
    00209   { return numberL_;}
    00210
    00212   inline CoinBigIndex baseL() const
    00213   { return baseL_;}
    00215   inline int maximumRowsExtra (  ) const {
    00216     return maximumRowsExtra_;
    00217   }
    00219   inline int numberColumns (  ) const {
    00220     return numberColumns_;
    00221   }
    00223   inline int numberElements (  ) const {
    00224     return totalElements_;
    00225   }
    00227   inline int numberForrestTomlin (  ) const {
    00228     return numberInColumn_.array()[numberColumnsExtra_];
    00229   }
    00231   inline int numberGoodColumns (  ) const {
    00232     return numberGoodU_;
    00233   }
    00235   inline double areaFactor (  ) const {
    00236     return areaFactor_;
    00237   }
    00238   inline void areaFactor ( double value ) {
    00239     areaFactor_=value;
    00240   }
    00242   double adjustedAreaFactor() const;
    00244   inline void relaxAccuracyCheck(double value)
    00245   { relaxCheck_ = value;}
    00246   inline double getAccuracyCheck() const
    00247   { return relaxCheck_;}
    00249   inline int messageLevel (  ) const {
    00250     return messageLevel_ ;
    00251   }
    00252   void messageLevel (  int value );
    00254   inline int maximumPivots (  ) const {
    00255     return maximumPivots_ ;
    00256   }
    00257   void maximumPivots (  int value );
    00258
    00260   inline int denseThreshold() const
    00261     { return denseThreshold_;}
    00263   inline void setDenseThreshold(int value)
    00264     { denseThreshold_ = value;}
    00266   inline double pivotTolerance (  ) const {
    00267     return pivotTolerance_ ;
    00268   }
    00269   void pivotTolerance (  double value );
    00271   inline double zeroTolerance (  ) const {
    00272     return zeroTolerance_ ;
    00273   }
    00274   void zeroTolerance (  double value );
    00275 #ifndef COIN_FAST_CODE
    00277   inline double slackValue (  ) const {
    00278     return slackValue_ ;
    00279   }
    00280   void slackValue (  double value );
    00281 #endif
    00283   double maximumCoefficient() const;
    00285   inline bool forrestTomlin() const
    00286   { return doForrestTomlin_;}
    00287   inline void setForrestTomlin(bool value)
    00288   { doForrestTomlin_=value;}
    00290   inline bool spaceForForrestTomlin() const
    00291   {
    00292     CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
    00293     CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
    00294     return (space>=0)&&doForrestTomlin_;
    00295   }
    00297
    00300
    00302   inline int numberDense() const
    00303   { return numberDense_;}
    00304
    00306   inline CoinBigIndex numberElementsU (  ) const {
    00307     return lengthU_;
    00308   }
    00310   inline void setNumberElementsU(CoinBigIndex value)
    00311   { lengthU_ = value; }
    00313   inline CoinBigIndex lengthAreaU (  ) const {
    00314     return lengthAreaU_;
    00315   }
    00317   inline CoinBigIndex numberElementsL (  ) const {
    00318     return lengthL_;
    00319   }
    00321   inline CoinBigIndex lengthAreaL (  ) const {
    00322     return lengthAreaL_;
    00323   }
    00325   inline CoinBigIndex numberElementsR (  ) const {
    00326     return lengthR_;
    00327   }
    00329   inline CoinBigIndex numberCompressions() const
    00330   { return numberCompressions_;}
    00332   inline int * numberInRow() const
    00333   { return numberInRow_.array();}
    00335   inline int * numberInColumn() const
    00336   { return numberInColumn_.array();}
    00338   inline CoinFactorizationDouble * elementU() const
    00339   { return elementU_.array();}
    00341   inline int * indexRowU() const
    00342   { return indexRowU_.array();}
    00344   inline CoinBigIndex * startColumnU() const
    00345   { return startColumnU_.array();}
    00347   inline int maximumColumnsExtra()
    00348   { return maximumColumnsExtra_;}
    00352   inline int biasLU() const
    00353   { return biasLU_;}
    00354   inline void setBiasLU(int value)
    00355   { biasLU_=value;}
    00361   inline int persistenceFlag() const
    00362   { return persistenceFlag_;}
    00363   void setPersistenceFlag(int value);
    00365
    00368
    00376   int replaceColumn ( CoinIndexedVector * regionSparse,
    00377                       int pivotRow,
    00378                       double pivotCheck ,
    00379                       bool checkBeforeModifying=false,
    00380                       double acceptablePivot=1.0e-8);
    00385   void replaceColumnU ( CoinIndexedVector * regionSparse,
    00386                         CoinBigIndex * deleted,
    00387                         int internalPivotRow);
    00388 #ifdef ABC_USE_COIN_FACTORIZATION
    00389
    00391   CoinIndexedVector * fakeVector(CoinIndexedVector * vector,
    00392                                  int already=0) const;
    00393   void deleteFakeVector(CoinIndexedVector * vector,
    00394                         CoinIndexedVector * fakeVector) const;
    00399   double checkReplacePart1 (  CoinIndexedVector * regionSparse,
    00400                                      int pivotRow);
    00405   double checkReplacePart1 (  CoinIndexedVector * regionSparse,
    00406                                       CoinIndexedVector * partialUpdate,
    00407                                      int pivotRow);
    00410   int checkReplacePart2 ( int pivotRow,
    00411                                   double btranAlpha,
    00412                                   double ftranAlpha,
    00413                                   double ftAlpha,
    00414                                   double acceptablePivot = 1.0e-8);
    00417   void replaceColumnPart3 ( CoinIndexedVector * regionSparse,
    00418                             int pivotRow,
    00419                             double alpha );
    00422   void replaceColumnPart3 ( CoinIndexedVector * regionSparse,
    00423                             CoinIndexedVector * partialUpdate,
    00424                             int pivotRow,
    00425                             double alpha );
    00433   int updateColumnFT ( CoinIndexedVector & regionSparse);
    00434   int updateColumnFTPart1 ( CoinIndexedVector & regionSparse) ;
    00435   void updateColumnFTPart2 ( CoinIndexedVector & regionSparse) ;
    00439   void updateColumnFT ( CoinIndexedVector & regionSparseFT,
    00440                         CoinIndexedVector & partialUpdate,
    00441                         int which);
    00443   int updateColumn ( CoinIndexedVector & regionSparse) const;
    00448   int updateTwoColumnsFT ( CoinIndexedVector & regionSparseFT,
    00449                            CoinIndexedVector & regionSparseOther);
    00451   int updateColumnTranspose ( CoinIndexedVector & regionSparse) const;
    00453   void updateColumnCpu ( CoinIndexedVector & regionSparse,int whichCpu) const;
    00455   void updateColumnTransposeCpu ( CoinIndexedVector & regionSparse,int whichCpu) const;
    00457   void updateFullColumn ( CoinIndexedVector & regionSparse) const;
    00459   void updateFullColumnTranspose ( CoinIndexedVector & regionSparse) const;
    00461   void updateWeights ( CoinIndexedVector & regionSparse) const;
    00463   inline bool wantsTableauColumn() const
    00464   {return false;}
    00466   inline double minimumPivotTolerance (  ) const {
    00467     return pivotTolerance_ ;
    00468   }
    00469   inline void minimumPivotTolerance (  double value )
    00470   { pivotTolerance(value);}
    00472   inline void setParallelMode(int value)
    00473   { parallelMode_=value;}
    00475   inline void setSolveMode(int value)
    00476   { parallelMode_ &= 3;parallelMode_ |= (value<<2);}
    00478   inline int solveMode() const
    00479   { return parallelMode_ >> 2;}
    00481   void updatePartialUpdate(CoinIndexedVector & partialUpdate);
    00483   void makeNonSingular(int *  COIN_RESTRICT sequence);
    00484 #endif
    00485
    00486
    00496   int updateColumnFT ( CoinIndexedVector * regionSparse,
    00497                        CoinIndexedVector * regionSparse2);
    00500   int updateColumn ( CoinIndexedVector * regionSparse,
    00501                      CoinIndexedVector * regionSparse2,
    00502                      bool noPermute=false) const;
    00508   int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
    00509                            CoinIndexedVector * regionSparse2,
    00510                            CoinIndexedVector * regionSparse3,
    00511                            bool noPermuteRegion3=false) ;
    00516   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
    00517                               CoinIndexedVector * regionSparse2) const;
    00519   void goSparse();
    00521   inline int sparseThreshold ( ) const
    00522   { return sparseThreshold_;}
    00524   void sparseThreshold ( int value );
    00526
    00527
    00531
    00532   inline void clearArrays()
    00533   { gutsOfDestructor();}
    00535
    00538
    00541   int add ( CoinBigIndex numberElements,
    00542                int indicesRow[],
    00543                int indicesColumn[], double elements[] );
    00544
    00547   int addColumn ( CoinBigIndex numberElements,
    00548                      int indicesRow[], double elements[] );
    00549
    00552   int addRow ( CoinBigIndex numberElements,
    00553                   int indicesColumn[], double elements[] );
    00554
    00556   int deleteColumn ( int Row );
    00558   int deleteRow ( int Row );
    00559
    00563   int replaceRow ( int whichRow, int numberElements,
    00564                       const int indicesColumn[], const double elements[] );
    00566   void emptyRows(int numberToEmpty, const int which[]);
    00568
    00569
    00570   void checkSparse();
    00572 #if 0 //def CLP_FACTORIZATION_INSTRUMENT
    00573   inline bool collectStatistics() const
    00574   { return collectStatistics_;}
    00576   inline void setCollectStatistics(bool onOff) const
    00577   { collectStatistics_ = onOff;}
    00578 #else
    00579   inline bool collectStatistics() const
    00580   { return true;}
    00582   inline void setCollectStatistics(bool onOff) const
    00583   { }
    00584 #endif
    00586   void gutsOfDestructor(int type=1);
    00588   void gutsOfInitialize(int type);
    00589   void gutsOfCopy(const CoinFactorization &other);
    00590
    00592   void resetStatistics();
    00593
    00594
    00596
    00598
    00599   void getAreas ( int numberRows,
    00600                   int numberColumns,
    00601                   CoinBigIndex maximumL,
    00602                   CoinBigIndex maximumU );
    00603
    00606   void preProcess ( int state,
    00607                     int possibleDuplicates = -1 );
    00609   int factor (  );
    00610 protected:
    00613   int factorSparse (  );
    00616   int factorSparseSmall (  );
    00619   int factorSparseLarge (  );
    00622   int factorDense (  );
    00623
    00625   bool pivotOneOtherRow ( int pivotRow,
    00626                           int pivotColumn );
    00628   bool pivotRowSingleton ( int pivotRow,
    00629                            int pivotColumn );
    00631   bool pivotColumnSingleton ( int pivotRow,
    00632                               int pivotColumn );
    00633
    00638   bool getColumnSpace ( int iColumn,
    00639                         int extraNeeded );
    00640
    00643   bool reorderU();
    00647   bool getColumnSpaceIterateR ( int iColumn, double value,
    00648                                int iRow);
    00654   CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
    00655                                int iRow);
    00659   bool getRowSpace ( int iRow, int extraNeeded );
    00660
    00664   bool getRowSpaceIterate ( int iRow,
    00665                             int extraNeeded );
    00667   void checkConsistency (  );
    00669   inline void addLink ( int index, int count ) {
    00670     int *nextCount = nextCount_.array();
    00671     int *firstCount = firstCount_.array();
    00672     int *lastCount = lastCount_.array();
    00673     int next = firstCount[count];
    00674       lastCount[index] = -2 - count;
    00675     if ( next < 0 ) {
    00676       //first with that count
    00677       firstCount[count] = index;
    00678       nextCount[index] = -1;
    00679     } else {
    00680       firstCount[count] = index;
    00681       nextCount[index] = next;
    00682       lastCount[next] = index;
    00683   }}
    00685   inline void deleteLink ( int index ) {
    00686     int *nextCount = nextCount_.array();
    00687     int *firstCount = firstCount_.array();
    00688     int *lastCount = lastCount_.array();
    00689     int next = nextCount[index];
    00690     int last = lastCount[index];
    00691     if ( last >= 0 ) {
    00692       nextCount[last] = next;
    00693     } else {
    00694       int count = -last - 2;
    00695
    00696       firstCount[count] = next;
    00697     }
    00698     if ( next >= 0 ) {
    00699       lastCount[next] = last;
    00700     }
    00701     nextCount[index] = -2;
    00702     lastCount[index] = -2;
    00703     return;
    00704   }
    00706   void separateLinks(int count,bool rowsFirst);
    00708   void cleanup (  );
    00709
    00711   void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
    00713   void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
    00715   void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
    00717   void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
    00718
    00720   void updateColumnR ( CoinIndexedVector * region ) const;
    00723   void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
    00724
    00726   void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
    00727
    00729   void updateColumnUSparse ( CoinIndexedVector * regionSparse,
    00730                              int * indexIn) const;
    00732   void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
    00733                                int * indexIn) const;
    00735   int updateColumnUDensish ( double * COIN_RESTRICT region,
    00736                              int * COIN_RESTRICT regionIndex) const;
    00738   void updateTwoColumnsUDensish (
    00739                                  int & numberNonZero1,
    00740                                  double * COIN_RESTRICT region1,
    00741                                  int * COIN_RESTRICT index1,
    00742                                  int & numberNonZero2,
    00743                                  double * COIN_RESTRICT region2,
    00744                                  int * COIN_RESTRICT index2) const;
    00746   void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
    00748   void permuteBack ( CoinIndexedVector * regionSparse,
    00749                      CoinIndexedVector * outVector) const;
    00750
    00752   void updateColumnTransposePFI ( CoinIndexedVector * region) const;
    00755   void updateColumnTransposeU ( CoinIndexedVector * region,
    00756                                 int smallestIndex) const;
    00759   void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
    00760                                         int smallestIndex) const;
    00763   void updateColumnTransposeUDensish ( CoinIndexedVector * region,
    00764                                        int smallestIndex) const;
    00767   void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
    00770   void updateColumnTransposeUByColumn ( CoinIndexedVector * region,
    00771                                         int smallestIndex) const;
    00772
    00774   void updateColumnTransposeR ( CoinIndexedVector * region ) const;
    00776   void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
    00778   void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
    00779
    00781   void updateColumnTransposeL ( CoinIndexedVector * region ) const;
    00783   void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
    00785   void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
    00787   void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
    00789   void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
    00790 public:
    00795   int replaceColumnPFI ( CoinIndexedVector * regionSparse,
    00796                          int pivotRow, double alpha);
    00797 protected:
    00800   int checkPivot(double saveFromU, double oldPivot) const;
    00801   /********************************* START LARGE TEMPLATE ********/
    00802 #ifdef INT_IS_8
    00803 #define COINFACTORIZATION_BITS_PER_INT 64
    00804 #define COINFACTORIZATION_SHIFT_PER_INT 6
    00805 #define COINFACTORIZATION_MASK_PER_INT 0x3f
    00806 #else
    00807 #define COINFACTORIZATION_BITS_PER_INT 32
    00808 #define COINFACTORIZATION_SHIFT_PER_INT 5
    00809 #define COINFACTORIZATION_MASK_PER_INT 0x1f
    00810 #endif
    00811   template <class T>  inline bool
    00812   pivot ( int pivotRow,
    00813           int pivotColumn,
    00814           CoinBigIndex pivotRowPosition,
    00815           CoinBigIndex pivotColumnPosition,
    00816           CoinFactorizationDouble work[],
    00817           unsigned int workArea2[],
    00818           int increment2,
    00819           T markRow[] ,
    00820           int largeInteger)
    00821 {
    00822   int *indexColumnU = indexColumnU_.array();
    00823   CoinBigIndex *startColumnU = startColumnU_.array();
    00824   int *numberInColumn = numberInColumn_.array();
    00825   CoinFactorizationDouble *elementU = elementU_.array();
    00826   int *indexRowU = indexRowU_.array();
    00827   CoinBigIndex *startRowU = startRowU_.array();
    00828   int *numberInRow = numberInRow_.array();
    00829   CoinFactorizationDouble *elementL = elementL_.array();
    00830   int *indexRowL = indexRowL_.array();
    00831   int *saveColumn = saveColumn_.array();
    00832   int *nextRow = nextRow_.array();
    00833   int *lastRow = lastRow_.array() ;
    00834
    00835   //store pivot columns (so can easily compress)
    00836   int numberInPivotRow = numberInRow[pivotRow] - 1;
    00837   CoinBigIndex startColumn = startColumnU[pivotColumn];
    00838   int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
    00839   CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
    00840   int put = 0;
    00841   CoinBigIndex startRow = startRowU[pivotRow];
    00842   CoinBigIndex endRow = startRow + numberInPivotRow + 1;
    00843
    00844   if ( pivotColumnPosition < 0 ) {
    00845     for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    00846       int iColumn = indexColumnU[pivotColumnPosition];
    00847       if ( iColumn != pivotColumn ) {
    00848         saveColumn[put++] = iColumn;
    00849       } else {
    00850         break;
    00851       }
    00852     }
    00853   } else {
    00854     for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
    00855       saveColumn[put++] = indexColumnU[i];
    00856     }
    00857   }
    00858   assert (pivotColumnPosition<endRow);
    00859   assert (indexColumnU[pivotColumnPosition]==pivotColumn);
    00860   pivotColumnPosition++;
    00861   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    00862     saveColumn[put++] = indexColumnU[pivotColumnPosition];
    00863   }
    00864   //take out this bit of indexColumnU
    00865   int next = nextRow[pivotRow];
    00866   int last = lastRow[pivotRow];
    00867
    00868   nextRow[last] = next;
    00869   lastRow[next] = last;
    00870   nextRow[pivotRow] = numberGoodU_;     //use for permute
    00871   lastRow[pivotRow] = -2;
    00872   numberInRow[pivotRow] = 0;
    00873   //store column in L, compress in U and take column out
    00874   CoinBigIndex l = lengthL_;
    00875
    00876   if ( l + numberInPivotColumn > lengthAreaL_ ) {
    00877     //need more memory
    00878     if ((messageLevel_&4)!=0)
    00879       printf("more memory needed in middle of invert\n");
    00880     return false;
    00881   }
    00882   //l+=currentAreaL_->elementByColumn-elementL;
    00883   CoinBigIndex lSave = l;
    00884
    00885   CoinBigIndex * startColumnL = startColumnL_.array();
    00886   startColumnL[numberGoodL_] = l;       //for luck and first time
    00887   numberGoodL_++;
    00888   startColumnL[numberGoodL_] = l + numberInPivotColumn;
    00889   lengthL_ += numberInPivotColumn;
    00890   if ( pivotRowPosition < 0 ) {
    00891     for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    00892       int iRow = indexRowU[pivotRowPosition];
    00893       if ( iRow != pivotRow ) {
    00894         indexRowL[l] = iRow;
    00895         elementL[l] = elementU[pivotRowPosition];
    00896         markRow[iRow] = static_cast<T>(l - lSave);
    00897         l++;
    00898         //take out of row list
    00899         CoinBigIndex start = startRowU[iRow];
    00900         CoinBigIndex end = start + numberInRow[iRow];
    00901         CoinBigIndex where = start;
    00902
    00903         while ( indexColumnU[where] != pivotColumn ) {
    00904           where++;
    00905         }                       /* endwhile */
    00906 #if DEBUG_COIN
    00907         if ( where >= end ) {
    00908           abort (  );
    00909         }
    00910 #endif
    00911         indexColumnU[where] = indexColumnU[end - 1];
    00912         numberInRow[iRow]--;
    00913       } else {
    00914         break;
    00915       }
    00916     }
    00917   } else {
    00918     CoinBigIndex i;
    00919
    00920     for ( i = startColumn; i < pivotRowPosition; i++ ) {
    00921       int iRow = indexRowU[i];
    00922
    00923       markRow[iRow] = static_cast<T>(l - lSave);
    00924       indexRowL[l] = iRow;
    00925       elementL[l] = elementU[i];
    00926       l++;
    00927       //take out of row list
    00928       CoinBigIndex start = startRowU[iRow];
    00929       CoinBigIndex end = start + numberInRow[iRow];
    00930       CoinBigIndex where = start;
    00931
    00932       while ( indexColumnU[where] != pivotColumn ) {
    00933         where++;
    00934       }                         /* endwhile */
    00935 #if DEBUG_COIN
    00936       if ( where >= end ) {
    00937         abort (  );
    00938       }
    00939 #endif
    00940       indexColumnU[where] = indexColumnU[end - 1];
    00941       numberInRow[iRow]--;
    00942       assert (numberInRow[iRow]>=0);
    00943     }
    00944   }
    00945   assert (pivotRowPosition<endColumn);
    00946   assert (indexRowU[pivotRowPosition]==pivotRow);
    00947   CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
    00948   CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
    00949
    00950   pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
    00951   pivotRowPosition++;
    00952   for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    00953     int iRow = indexRowU[pivotRowPosition];
    00954
    00955     markRow[iRow] = static_cast<T>(l - lSave);
    00956     indexRowL[l] = iRow;
    00957     elementL[l] = elementU[pivotRowPosition];
    00958     l++;
    00959     //take out of row list
    00960     CoinBigIndex start = startRowU[iRow];
    00961     CoinBigIndex end = start + numberInRow[iRow];
    00962     CoinBigIndex where = start;
    00963
    00964     while ( indexColumnU[where] != pivotColumn ) {
    00965       where++;
    00966     }                           /* endwhile */
    00967 #if DEBUG_COIN
    00968     if ( where >= end ) {
    00969       abort (  );
    00970     }
    00971 #endif
    00972     indexColumnU[where] = indexColumnU[end - 1];
    00973     numberInRow[iRow]--;
    00974     assert (numberInRow[iRow]>=0);
    00975   }
    00976   markRow[pivotRow] = static_cast<T>(largeInteger);
    00977   //compress pivot column (move pivot to front including saved)
    00978   numberInColumn[pivotColumn] = 0;
    00979   //use end of L for temporary space
    00980   int *indexL = &indexRowL[lSave];
    00981   CoinFactorizationDouble *multipliersL = &elementL[lSave];
    00982
    00983   //adjust
    00984   int j;
    00985
    00986   for ( j = 0; j < numberInPivotColumn; j++ ) {
    00987     multipliersL[j] *= pivotMultiplier;
    00988   }
    00989   //zero out fill
    00990   CoinBigIndex iErase;
    00991   for ( iErase = 0; iErase < increment2 * numberInPivotRow;
    00992         iErase++ ) {
    00993     workArea2[iErase] = 0;
    00994   }
    00995   CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
    00996   unsigned int *temp2 = workArea2;
    00997   int * nextColumn = nextColumn_.array();
    00998
    00999   //pack down and move to work
    01000   int jColumn;
    01001   for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01002     int iColumn = saveColumn[jColumn];
    01003     CoinBigIndex startColumn = startColumnU[iColumn];
    01004     CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
    01005     int iRow = indexRowU[startColumn];
    01006     CoinFactorizationDouble value = elementU[startColumn];
    01007     double largest;
    01008     CoinBigIndex put = startColumn;
    01009     CoinBigIndex positionLargest = -1;
    01010     CoinFactorizationDouble thisPivotValue = 0.0;
    01011
    01012     //compress column and find largest not updated
    01013     bool checkLargest;
    01014     int mark = markRow[iRow];
    01015
    01016     if ( mark == largeInteger+1 ) {
    01017       largest = fabs ( value );
    01018       positionLargest = put;
    01019       put++;
    01020       checkLargest = false;
    01021     } else {
    01022       //need to find largest
    01023       largest = 0.0;
    01024       checkLargest = true;
    01025       if ( mark != largeInteger ) {
    01026         //will be updated
    01027         work[mark] = value;
    01028         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    01029         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    01030
    01031         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
    01032         added--;
    01033       } else {
    01034         thisPivotValue = value;
    01035       }
    01036     }
    01037     CoinBigIndex i;
    01038     for ( i = startColumn + 1; i < endColumn; i++ ) {
    01039       iRow = indexRowU[i];
    01040       value = elementU[i];
    01041       int mark = markRow[iRow];
    01042
    01043       if ( mark == largeInteger+1 ) {
    01044         //keep
    01045         indexRowU[put] = iRow;
    01046         elementU[put] = value;
    01047         if ( checkLargest ) {
    01048           double absValue = fabs ( value );
    01049
    01050           if ( absValue > largest ) {
    01051             largest = absValue;
    01052             positionLargest = put;
    01053           }
    01054         }
    01055         put++;
    01056       } else if ( mark != largeInteger ) {
    01057         //will be updated
    01058         work[mark] = value;
    01059         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    01060         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    01061
    01062         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
    01063         added--;
    01064       } else {
    01065         thisPivotValue = value;
    01066       }
    01067     }
    01068     //slot in pivot
    01069     elementU[put] = elementU[startColumn];
    01070     indexRowU[put] = indexRowU[startColumn];
    01071     if ( positionLargest == startColumn ) {
    01072       positionLargest = put;    //follow if was largest
    01073     }
    01074     put++;
    01075     elementU[startColumn] = thisPivotValue;
    01076     indexRowU[startColumn] = pivotRow;
    01077     //clean up counts
    01078     startColumn++;
    01079     numberInColumn[iColumn] = put - startColumn;
    01080     int * numberInColumnPlus = numberInColumnPlus_.array();
    01081     numberInColumnPlus[iColumn]++;
    01082     startColumnU[iColumn]++;
    01083     //how much space have we got
    01084     int next = nextColumn[iColumn];
    01085     CoinBigIndex space;
    01086
    01087     space = startColumnU[next] - put - numberInColumnPlus[next];
    01088     //assume no zero elements
    01089     if ( numberInPivotColumn > space ) {
    01090       //getColumnSpace also moves fixed part
    01091       if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
    01092         return false;
    01093       }
    01094       //redo starts
    01095       if (positionLargest >= 0)
    01096          positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
    01097       startColumn = startColumnU[iColumn];
    01098       put = startColumn + numberInColumn[iColumn];
    01099     }
    01100     double tolerance = zeroTolerance_;
    01101
    01102     int *nextCount = nextCount_.array();
    01103     for ( j = 0; j < numberInPivotColumn; j++ ) {
    01104       value = work[j] - thisPivotValue * multipliersL[j];
    01105       double absValue = fabs ( value );
    01106
    01107       if ( absValue > tolerance ) {
    01108         work[j] = 0.0;
    01109         assert (put<lengthAreaU_);
    01110         elementU[put] = value;
    01111         indexRowU[put] = indexL[j];
    01112         if ( absValue > largest ) {
    01113           largest = absValue;
    01114           positionLargest = put;
    01115         }
    01116         put++;
    01117       } else {
    01118         work[j] = 0.0;
    01119         added--;
    01120         int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01121         int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01122
    01123         if ( temp2[word] & ( 1 << bit ) ) {
    01124           //take out of row list
    01125           iRow = indexL[j];
    01126           CoinBigIndex start = startRowU[iRow];
    01127           CoinBigIndex end = start + numberInRow[iRow];
    01128           CoinBigIndex where = start;
    01129
    01130           while ( indexColumnU[where] != iColumn ) {
    01131             where++;
    01132           }                     /* endwhile */
    01133 #if DEBUG_COIN
    01134           if ( where >= end ) {
    01135             abort (  );
    01136           }
    01137 #endif
    01138           indexColumnU[where] = indexColumnU[end - 1];
    01139           numberInRow[iRow]--;
    01140         } else {
    01141           //make sure won't be added
    01142           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01143           int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01144
    01145           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    01146         }
    01147       }
    01148     }
    01149     numberInColumn[iColumn] = put - startColumn;
    01150     //move largest
    01151     if ( positionLargest >= 0 ) {
    01152       value = elementU[positionLargest];
    01153       iRow = indexRowU[positionLargest];
    01154       elementU[positionLargest] = elementU[startColumn];
    01155       indexRowU[positionLargest] = indexRowU[startColumn];
    01156       elementU[startColumn] = value;
    01157       indexRowU[startColumn] = iRow;
    01158     }
    01159     //linked list for column
    01160     if ( nextCount[iColumn + numberRows_] != -2 ) {
    01161       //modify linked list
    01162       deleteLink ( iColumn + numberRows_ );
    01163       addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
    01164     }
    01165     temp2 += increment2;
    01166   }
    01167   //get space for row list
    01168   unsigned int *putBase = workArea2;
    01169   int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
    01170   int i = 0;
    01171
    01172   // do linked lists and update counts
    01173   while ( bigLoops ) {
    01174     bigLoops--;
    01175     int bit;
    01176     for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
    01177       unsigned int *putThis = putBase;
    01178       int iRow = indexL[i];
    01179
    01180       //get space
    01181       int number = 0;
    01182       int jColumn;
    01183
    01184       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01185         unsigned int test = *putThis;
    01186
    01187         putThis += increment2;
    01188         test = 1 - ( ( test >> bit ) & 1 );
    01189         number += test;
    01190       }
    01191       int next = nextRow[iRow];
    01192       CoinBigIndex space;
    01193
    01194       space = startRowU[next] - startRowU[iRow];
    01195       number += numberInRow[iRow];
    01196       if ( space < number ) {
    01197         if ( !getRowSpace ( iRow, number ) ) {
    01198           return false;
    01199         }
    01200       }
    01201       // now do
    01202       putThis = putBase;
    01203       next = nextRow[iRow];
    01204       number = numberInRow[iRow];
    01205       CoinBigIndex end = startRowU[iRow] + number;
    01206       int saveIndex = indexColumnU[startRowU[next]];
    01207
    01208       //add in
    01209       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01210         unsigned int test = *putThis;
    01211
    01212         putThis += increment2;
    01213         test = 1 - ( ( test >> bit ) & 1 );
    01214         indexColumnU[end] = saveColumn[jColumn];
    01215         end += test;
    01216       }
    01217       //put back next one in case zapped
    01218       indexColumnU[startRowU[next]] = saveIndex;
    01219       markRow[iRow] = static_cast<T>(largeInteger+1);
    01220       number = end - startRowU[iRow];
    01221       numberInRow[iRow] = number;
    01222       deleteLink ( iRow );
    01223       addLink ( iRow, number );
    01224     }
    01225     putBase++;
    01226   }                             /* endwhile */
    01227   int bit;
    01228
    01229   for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
    01230     unsigned int *putThis = putBase;
    01231     int iRow = indexL[i];
    01232
    01233     //get space
    01234     int number = 0;
    01235     int jColumn;
    01236
    01237     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01238       unsigned int test = *putThis;
    01239
    01240       putThis += increment2;
    01241       test = 1 - ( ( test >> bit ) & 1 );
    01242       number += test;
    01243     }
    01244     int next = nextRow[iRow];
    01245     CoinBigIndex space;
    01246
    01247     space = startRowU[next] - startRowU[iRow];
    01248     number += numberInRow[iRow];
    01249     if ( space < number ) {
    01250       if ( !getRowSpace ( iRow, number ) ) {
    01251         return false;
    01252       }
    01253     }
    01254     // now do
    01255     putThis = putBase;
    01256     next = nextRow[iRow];
    01257     number = numberInRow[iRow];
    01258     CoinBigIndex end = startRowU[iRow] + number;
    01259     int saveIndex;
    01260
    01261     saveIndex = indexColumnU[startRowU[next]];
    01262
    01263     //add in
    01264     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    01265       unsigned int test = *putThis;
    01266
    01267       putThis += increment2;
    01268       test = 1 - ( ( test >> bit ) & 1 );
    01269
    01270       indexColumnU[end] = saveColumn[jColumn];
    01271       end += test;
    01272     }
    01273     indexColumnU[startRowU[next]] = saveIndex;
    01274     markRow[iRow] = static_cast<T>(largeInteger+1);
    01275     number = end - startRowU[iRow];
    01276     numberInRow[iRow] = number;
    01277     deleteLink ( iRow );
    01278     addLink ( iRow, number );
    01279   }
    01280   markRow[pivotRow] = static_cast<T>(largeInteger+1);
    01281   //modify linked list for pivots
    01282   deleteLink ( pivotRow );
    01283   deleteLink ( pivotColumn + numberRows_ );
    01284   totalElements_ += added;
    01285   return true;
    01286 }
    01287
    01288   /********************************* END LARGE TEMPLATE ********/
    01290
    01291 protected:
    01292
    01295
    01296   double pivotTolerance_;
    01298   double zeroTolerance_;
    01299 #ifndef COIN_FAST_CODE
    01301   double slackValue_;
    01302 #else
    01303 #ifndef slackValue_
    01304 #define slackValue_ -1.0
    01305 #endif
    01306 #endif
    01308   double areaFactor_;
    01310   double relaxCheck_;
    01312   int numberRows_;
    01314   int numberRowsExtra_;
    01316   int maximumRowsExtra_;
    01318   int numberColumns_;
    01320   int numberColumnsExtra_;
    01322   int maximumColumnsExtra_;
    01324   int numberGoodU_;
    01326   int numberGoodL_;
    01328   int maximumPivots_;
    01330   int numberPivots_;
    01333   CoinBigIndex totalElements_;
    01335   CoinBigIndex factorElements_;
    01337   CoinIntArrayWithLength pivotColumn_;
    01339   CoinIntArrayWithLength permute_;
    01341   CoinIntArrayWithLength permuteBack_;
    01343   CoinIntArrayWithLength pivotColumnBack_;
    01345   int status_;
    01346
    01351   //int increasingRows_;
    01352
    01354   int numberTrials_;
    01356   CoinBigIndexArrayWithLength startRowU_;
    01357
    01359   CoinIntArrayWithLength numberInRow_;
    01360
    01362   CoinIntArrayWithLength numberInColumn_;
    01363
    01365   CoinIntArrayWithLength numberInColumnPlus_;
    01366
    01369   CoinIntArrayWithLength firstCount_;
    01370
    01372   CoinIntArrayWithLength nextCount_;
    01373
    01375   CoinIntArrayWithLength lastCount_;
    01376
    01378   CoinIntArrayWithLength nextColumn_;
    01379
    01381   CoinIntArrayWithLength lastColumn_;
    01382
    01384   CoinIntArrayWithLength nextRow_;
    01385
    01387   CoinIntArrayWithLength lastRow_;
    01388
    01390   CoinIntArrayWithLength saveColumn_;
    01391
    01393   CoinIntArrayWithLength markRow_;
    01394
    01396   int messageLevel_;
    01397
    01399   int biggerDimension_;
    01400
    01402   CoinIntArrayWithLength indexColumnU_;
    01403
    01405   CoinIntArrayWithLength pivotRowL_;
    01406
    01408   CoinFactorizationDoubleArrayWithLength pivotRegion_;
    01409
    01411   int numberSlacks_;
    01412
    01414   int numberU_;
    01415
    01417   CoinBigIndex maximumU_;
    01418
    01420   //int baseU_;
    01421
    01423   CoinBigIndex lengthU_;
    01424
    01426   CoinBigIndex lengthAreaU_;
    01427
    01429   CoinFactorizationDoubleArrayWithLength elementU_;
    01430
    01432   CoinIntArrayWithLength indexRowU_;
    01433
    01435   CoinBigIndexArrayWithLength startColumnU_;
    01436
    01438   CoinBigIndexArrayWithLength convertRowToColumnU_;
    01439
    01441   CoinBigIndex numberL_;
    01442
    01444   CoinBigIndex baseL_;
    01445
    01447   CoinBigIndex lengthL_;
    01448
    01450   CoinBigIndex lengthAreaL_;
    01451
    01453   CoinFactorizationDoubleArrayWithLength elementL_;
    01454
    01456   CoinIntArrayWithLength indexRowL_;
    01457
    01459   CoinBigIndexArrayWithLength startColumnL_;
    01460
    01462   bool doForrestTomlin_;
    01463
    01465   int numberR_;
    01466
    01468   CoinBigIndex lengthR_;
    01469
    01471   CoinBigIndex lengthAreaR_;
    01472
    01474   CoinFactorizationDouble *elementR_;
    01475
    01477   int *indexRowR_;
    01478
    01480   CoinBigIndexArrayWithLength startColumnR_;
    01481
    01483   double  * denseArea_;
    01484
    01486   double  * denseAreaAddress_;
    01487
    01489   int * densePermute_;
    01490
    01492   int numberDense_;
    01493
    01495   int denseThreshold_;
    01496
    01498   CoinFactorizationDoubleArrayWithLength workArea_;
    01499
    01501   CoinUnsignedIntArrayWithLength workArea2_;
    01502
    01504   CoinBigIndex numberCompressions_;
    01505
    01506 public:
    01508   mutable double ftranCountInput_;
    01509   mutable double ftranCountAfterL_;
    01510   mutable double ftranCountAfterR_;
    01511   mutable double ftranCountAfterU_;
    01512   mutable double btranCountInput_;
    01513   mutable double btranCountAfterU_;
    01514   mutable double btranCountAfterR_;
    01515   mutable double btranCountAfterL_;
    01516
    01518   mutable int numberFtranCounts_;
    01519   mutable int numberBtranCounts_;
    01520
    01522   double ftranAverageAfterL_;
    01523   double ftranAverageAfterR_;
    01524   double ftranAverageAfterU_;
    01525   double btranAverageAfterU_;
    01526   double btranAverageAfterR_;
    01527   double btranAverageAfterL_;
    01528 protected:
    01529
    01531 #if 0
    01532   mutable bool collectStatistics_;
    01533 #else
    01534 #define collectStatistics_ 1
    01535 #endif
    01536
    01538   int sparseThreshold_;
    01539
    01541   int sparseThreshold2_;
    01542
    01544   CoinBigIndexArrayWithLength startRowL_;
    01545
    01547   CoinIntArrayWithLength indexColumnL_;
    01548
    01550   CoinFactorizationDoubleArrayWithLength elementByRowL_;
    01551
    01553   mutable CoinIntArrayWithLength sparse_;
    01557   int biasLU_;
    01563   int persistenceFlag_;
    01564 #ifdef ABC_USE_COIN_FACTORIZATION
    01566   int parallelMode_;
    01567 #endif
    01568
    01569 };
    01570 // Dense coding
    01571 #ifdef COIN_HAS_LAPACK
    01572 #ifndef COIN_FACTORIZATION_DENSE_CODE
    01573 #define COIN_FACTORIZATION_DENSE_CODE 1
    01574 #endif
    01575 #endif
    01576 #ifdef COIN_FACTORIZATION_DENSE_CODE
    01577 /* Type of Fortran integer translated into C */
    01578 #ifndef ipfint
    01579 //typedef ipfint FORTRAN_INTEGER_TYPE ;
    01580 typedef int ipfint;
    01581 typedef const int cipfint;
    01582 #endif
    01583 #endif
    01584 #endif
    01585 // Extra for ugly include
    01586 #ifdef UGLY_COIN_FACTOR_CODING
    01587 #define FAC_UNSET (FAC_SET+1)
    01588 {
    01589   goodPivot=false;
    01590   //store pivot columns (so can easily compress)
    01591   CoinBigIndex startColumnThis = startColumn[iPivotColumn];
    01592   CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
    01593   int put = 0;
    01594   CoinBigIndex startRowThis = startRow[iPivotRow];
    01595   CoinBigIndex endRow = startRowThis + numberDoRow + 1;
    01596   if ( pivotColumnPosition < 0 ) {
    01597     for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    01598       int iColumn = indexColumn[pivotColumnPosition];
    01599       if ( iColumn != iPivotColumn ) {
    01600         saveColumn[put++] = iColumn;
    01601       } else {
    01602         break;
    01603       }
    01604     }
    01605   } else {
    01606     for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
    01607       saveColumn[put++] = indexColumn[i];
    01608     }
    01609   }
    01610   assert (pivotColumnPosition<endRow);
    01611   assert (indexColumn[pivotColumnPosition]==iPivotColumn);
    01612   pivotColumnPosition++;
    01613   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    01614     saveColumn[put++] = indexColumn[pivotColumnPosition];
    01615   }
    01616   //take out this bit of indexColumn
    01617   int next = nextRow[iPivotRow];
    01618   int last = lastRow[iPivotRow];
    01619
    01620   nextRow[last] = next;
    01621   lastRow[next] = last;
    01622   nextRow[iPivotRow] = numberGoodU_;    //use for permute
    01623   lastRow[iPivotRow] = -2;
    01624   numberInRow[iPivotRow] = 0;
    01625   //store column in L, compress in U and take column out
    01626   CoinBigIndex l = lengthL_;
    01627   // **** HORRID coding coming up but a goto seems best!
    01628   {
    01629     if ( l + numberDoColumn > lengthAreaL_ ) {
    01630       //need more memory
    01631       if ((messageLevel_&4)!=0)
    01632         printf("more memory needed in middle of invert\n");
    01633       goto BAD_PIVOT;
    01634     }
    01635     //l+=currentAreaL_->elementByColumn-elementL;
    01636     CoinBigIndex lSave = l;
    01637
    01638     CoinBigIndex * startColumnL = startColumnL_.array();
    01639     startColumnL[numberGoodL_] = l;     //for luck and first time
    01640     numberGoodL_++;
    01641     startColumnL[numberGoodL_] = l + numberDoColumn;
    01642     lengthL_ += numberDoColumn;
    01643     if ( pivotRowPosition < 0 ) {
    01644       for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    01645         int iRow = indexRow[pivotRowPosition];
    01646         if ( iRow != iPivotRow ) {
    01647           indexRowL[l] = iRow;
    01648           elementL[l] = element[pivotRowPosition];
    01649           markRow[iRow] = l - lSave;
    01650           l++;
    01651           //take out of row list
    01652           CoinBigIndex start = startRow[iRow];
    01653           CoinBigIndex end = start + numberInRow[iRow];
    01654           CoinBigIndex where = start;
    01655
    01656           while ( indexColumn[where] != iPivotColumn ) {
    01657             where++;
    01658           }                     /* endwhile */
    01659 #if DEBUG_COIN
    01660           if ( where >= end ) {
    01661             abort (  );
    01662           }
    01663 #endif
    01664           indexColumn[where] = indexColumn[end - 1];
    01665           numberInRow[iRow]--;
    01666         } else {
    01667           break;
    01668         }
    01669       }
    01670     } else {
    01671       CoinBigIndex i;
    01672
    01673       for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
    01674         int iRow = indexRow[i];
    01675
    01676         markRow[iRow] = l - lSave;
    01677         indexRowL[l] = iRow;
    01678         elementL[l] = element[i];
    01679         l++;
    01680         //take out of row list
    01681         CoinBigIndex start = startRow[iRow];
    01682         CoinBigIndex end = start + numberInRow[iRow];
    01683         CoinBigIndex where = start;
    01684
    01685         while ( indexColumn[where] != iPivotColumn ) {
    01686           where++;
    01687         }                               /* endwhile */
    01688 #if DEBUG_COIN
    01689         if ( where >= end ) {
    01690           abort (  );
    01691         }
    01692 #endif
    01693         indexColumn[where] = indexColumn[end - 1];
    01694         numberInRow[iRow]--;
    01695         assert (numberInRow[iRow]>=0);
    01696       }
    01697     }
    01698     assert (pivotRowPosition<endColumn);
    01699     assert (indexRow[pivotRowPosition]==iPivotRow);
    01700     CoinFactorizationDouble pivotElement = element[pivotRowPosition];
    01701     CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
    01702
    01703     pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
    01704     pivotRowPosition++;
    01705     for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    01706       int iRow = indexRow[pivotRowPosition];
    01707
    01708       markRow[iRow] = l - lSave;
    01709       indexRowL[l] = iRow;
    01710       elementL[l] = element[pivotRowPosition];
    01711       l++;
    01712       //take out of row list
    01713       CoinBigIndex start = startRow[iRow];
    01714       CoinBigIndex end = start + numberInRow[iRow];
    01715       CoinBigIndex where = start;
    01716
    01717       while ( indexColumn[where] != iPivotColumn ) {
    01718         where++;
    01719       }                         /* endwhile */
    01720 #if DEBUG_COIN
    01721       if ( where >= end ) {
    01722         abort (  );
    01723       }
    01724 #endif
    01725       indexColumn[where] = indexColumn[end - 1];
    01726       numberInRow[iRow]--;
    01727       assert (numberInRow[iRow]>=0);
    01728     }
    01729     markRow[iPivotRow] = FAC_SET;
    01730     //compress pivot column (move pivot to front including saved)
    01731     numberInColumn[iPivotColumn] = 0;
    01732     //use end of L for temporary space
    01733     int *indexL = &indexRowL[lSave];
    01734     CoinFactorizationDouble *multipliersL = &elementL[lSave];
    01735
    01736     //adjust
    01737     int j;
    01738
    01739     for ( j = 0; j < numberDoColumn; j++ ) {
    01740       multipliersL[j] *= pivotMultiplier;
    01741     }
    01742     //zero out fill
    01743     CoinBigIndex iErase;
    01744     for ( iErase = 0; iErase < increment2 * numberDoRow;
    01745           iErase++ ) {
    01746       workArea2[iErase] = 0;
    01747     }
    01748     CoinBigIndex added = numberDoRow * numberDoColumn;
    01749     unsigned int *temp2 = workArea2;
    01750     int * nextColumn = nextColumn_.array();
    01751
    01752     //pack down and move to work
    01753     int jColumn;
    01754     for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01755       int iColumn = saveColumn[jColumn];
    01756       CoinBigIndex startColumnThis = startColumn[iColumn];
    01757       CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
    01758       int iRow = indexRow[startColumnThis];
    01759       CoinFactorizationDouble value = element[startColumnThis];
    01760       double largest;
    01761       CoinBigIndex put = startColumnThis;
    01762       CoinBigIndex positionLargest = -1;
    01763       CoinFactorizationDouble thisPivotValue = 0.0;
    01764
    01765       //compress column and find largest not updated
    01766       bool checkLargest;
    01767       int mark = markRow[iRow];
    01768
    01769       if ( mark == FAC_UNSET ) {
    01770         largest = fabs ( value );
    01771         positionLargest = put;
    01772         put++;
    01773         checkLargest = false;
    01774       } else {
    01775         //need to find largest
    01776         largest = 0.0;
    01777         checkLargest = true;
    01778         if ( mark != FAC_SET ) {
    01779           //will be updated
    01780           workArea[mark] = value;
    01781           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    01782           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    01783
    01784           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    01785           added--;
    01786         } else {
    01787           thisPivotValue = value;
    01788         }
    01789       }
    01790       CoinBigIndex i;
    01791       for ( i = startColumnThis + 1; i < endColumn; i++ ) {
    01792         iRow = indexRow[i];
    01793         value = element[i];
    01794         int mark = markRow[iRow];
    01795
    01796         if ( mark == FAC_UNSET ) {
    01797           //keep
    01798           indexRow[put] = iRow;
    01799           element[put] = value;
    01800           if ( checkLargest ) {
    01801             double absValue = fabs ( value );
    01802
    01803             if ( absValue > largest ) {
    01804               largest = absValue;
    01805               positionLargest = put;
    01806             }
    01807           }
    01808           put++;
    01809         } else if ( mark != FAC_SET ) {
    01810           //will be updated
    01811           workArea[mark] = value;
    01812           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
    01813           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
    01814
    01815           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
    01816           added--;
    01817         } else {
    01818           thisPivotValue = value;
    01819         }
    01820       }
    01821       //slot in pivot
    01822       element[put] = element[startColumnThis];
    01823       indexRow[put] = indexRow[startColumnThis];
    01824       if ( positionLargest == startColumnThis ) {
    01825         positionLargest = put;  //follow if was largest
    01826       }
    01827       put++;
    01828       element[startColumnThis] = thisPivotValue;
    01829       indexRow[startColumnThis] = iPivotRow;
    01830       //clean up counts
    01831       startColumnThis++;
    01832       numberInColumn[iColumn] = put - startColumnThis;
    01833       int * numberInColumnPlus = numberInColumnPlus_.array();
    01834       numberInColumnPlus[iColumn]++;
    01835       startColumn[iColumn]++;
    01836       //how much space have we got
    01837       int next = nextColumn[iColumn];
    01838       CoinBigIndex space;
    01839
    01840       space = startColumn[next] - put - numberInColumnPlus[next];
    01841       //assume no zero elements
    01842       if ( numberDoColumn > space ) {
    01843         //getColumnSpace also moves fixed part
    01844         if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
    01845           goto BAD_PIVOT;
    01846         }
    01847         //redo starts
    01848         positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
    01849         startColumnThis = startColumn[iColumn];
    01850         put = startColumnThis + numberInColumn[iColumn];
    01851       }
    01852       double tolerance = zeroTolerance_;
    01853
    01854       int *nextCount = nextCount_.array();
    01855       for ( j = 0; j < numberDoColumn; j++ ) {
    01856         value = workArea[j] - thisPivotValue * multipliersL[j];
    01857         double absValue = fabs ( value );
    01858
    01859         if ( absValue > tolerance ) {
    01860           workArea[j] = 0.0;
    01861           element[put] = value;
    01862           indexRow[put] = indexL[j];
    01863           if ( absValue > largest ) {
    01864             largest = absValue;
    01865             positionLargest = put;
    01866           }
    01867           put++;
    01868         } else {
    01869           workArea[j] = 0.0;
    01870           added--;
    01871           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01872           int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01873
    01874           if ( temp2[word] & ( 1 << bit ) ) {
    01875             //take out of row list
    01876             iRow = indexL[j];
    01877             CoinBigIndex start = startRow[iRow];
    01878             CoinBigIndex end = start + numberInRow[iRow];
    01879             CoinBigIndex where = start;
    01880
    01881             while ( indexColumn[where] != iColumn ) {
    01882               where++;
    01883             }                   /* endwhile */
    01884 #if DEBUG_COIN
    01885             if ( where >= end ) {
    01886               abort (  );
    01887             }
    01888 #endif
    01889             indexColumn[where] = indexColumn[end - 1];
    01890             numberInRow[iRow]--;
    01891           } else {
    01892             //make sure won't be added
    01893             int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
    01894             int bit = j & COINFACTORIZATION_MASK_PER_INT;
    01895
    01896             temp2[word] = temp2[word] | ( 1 << bit );   //say already in counts
    01897           }
    01898         }
    01899       }
    01900       numberInColumn[iColumn] = put - startColumnThis;
    01901       //move largest
    01902       if ( positionLargest >= 0 ) {
    01903         value = element[positionLargest];
    01904         iRow = indexRow[positionLargest];
    01905         element[positionLargest] = element[startColumnThis];
    01906         indexRow[positionLargest] = indexRow[startColumnThis];
    01907         element[startColumnThis] = value;
    01908         indexRow[startColumnThis] = iRow;
    01909       }
    01910       //linked list for column
    01911       if ( nextCount[iColumn + numberRows_] != -2 ) {
    01912         //modify linked list
    01913         deleteLink ( iColumn + numberRows_ );
    01914         addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
    01915       }
    01916       temp2 += increment2;
    01917     }
    01918     //get space for row list
    01919     unsigned int *putBase = workArea2;
    01920     int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
    01921     int i = 0;
    01922
    01923     // do linked lists and update counts
    01924     while ( bigLoops ) {
    01925       bigLoops--;
    01926       int bit;
    01927       for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
    01928         unsigned int *putThis = putBase;
    01929         int iRow = indexL[i];
    01930
    01931         //get space
    01932         int number = 0;
    01933         int jColumn;
    01934
    01935         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01936           unsigned int test = *putThis;
    01937
    01938           putThis += increment2;
    01939           test = 1 - ( ( test >> bit ) & 1 );
    01940           number += test;
    01941         }
    01942         int next = nextRow[iRow];
    01943         CoinBigIndex space;
    01944
    01945         space = startRow[next] - startRow[iRow];
    01946         number += numberInRow[iRow];
    01947         if ( space < number ) {
    01948           if ( !getRowSpace ( iRow, number ) ) {
    01949             goto BAD_PIVOT;
    01950           }
    01951         }
    01952         // now do
    01953         putThis = putBase;
    01954         next = nextRow[iRow];
    01955         number = numberInRow[iRow];
    01956         CoinBigIndex end = startRow[iRow] + number;
    01957         int saveIndex = indexColumn[startRow[next]];
    01958
    01959         //add in
    01960         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01961           unsigned int test = *putThis;
    01962
    01963           putThis += increment2;
    01964           test = 1 - ( ( test >> bit ) & 1 );
    01965           indexColumn[end] = saveColumn[jColumn];
    01966           end += test;
    01967         }
    01968         //put back next one in case zapped
    01969         indexColumn[startRow[next]] = saveIndex;
    01970         markRow[iRow] = FAC_UNSET;
    01971         number = end - startRow[iRow];
    01972         numberInRow[iRow] = number;
    01973         deleteLink ( iRow );
    01974         addLink ( iRow, number );
    01975       }
    01976       putBase++;
    01977     }                           /* endwhile */
    01978     int bit;
    01979
    01980     for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
    01981       unsigned int *putThis = putBase;
    01982       int iRow = indexL[i];
    01983
    01984       //get space
    01985       int number = 0;
    01986       int jColumn;
    01987
    01988       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    01989         unsigned int test = *putThis;
    01990
    01991         putThis += increment2;
    01992         test = 1 - ( ( test >> bit ) & 1 );
    01993         number += test;
    01994       }
    01995       int next = nextRow[iRow];
    01996       CoinBigIndex space;
    01997
    01998       space = startRow[next] - startRow[iRow];
    01999       number += numberInRow[iRow];
    02000       if ( space < number ) {
    02001         if ( !getRowSpace ( iRow, number ) ) {
    02002           goto BAD_PIVOT;
    02003         }
    02004       }
    02005       // now do
    02006       putThis = putBase;
    02007       next = nextRow[iRow];
    02008       number = numberInRow[iRow];
    02009       CoinBigIndex end = startRow[iRow] + number;
    02010       int saveIndex;
    02011
    02012       saveIndex = indexColumn[startRow[next]];
    02013
    02014       //add in
    02015       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
    02016         unsigned int test = *putThis;
    02017
    02018         putThis += increment2;
    02019         test = 1 - ( ( test >> bit ) & 1 );
    02020
    02021         indexColumn[end] = saveColumn[jColumn];
    02022         end += test;
    02023       }
    02024       indexColumn[startRow[next]] = saveIndex;
    02025       markRow[iRow] = FAC_UNSET;
    02026       number = end - startRow[iRow];
    02027       numberInRow[iRow] = number;
    02028       deleteLink ( iRow );
    02029       addLink ( iRow, number );
    02030     }
    02031     markRow[iPivotRow] = FAC_UNSET;
    02032     //modify linked list for pivots
    02033     deleteLink ( iPivotRow );
    02034     deleteLink ( iPivotColumn + numberRows_ );
    02035     totalElements_ += added;
    02036     goodPivot= true;
    02037     // **** UGLY UGLY UGLY
    02038   }
    02039  BAD_PIVOT:
    02040
    02041   ;
    02042 }
    02043 #undef FAC_UNSET
    02044 #endif