/home/coin/SVN-release/Osi-0.102.2/CoinUtils/src/CoinFactorization.hpp

Go to the documentation of this file.
00001 /* $Id: CoinFactorization.hpp 1215 2009-11-05 11:03:04Z forrest $ */
00002 // Copyright (C) 2002, International Business Machines
00003 // Corporation and others.  All Rights Reserved.
00004 
00005 /* 
00006    Authors
00007    
00008    John Forrest
00009 
00010  */
00011 #ifndef CoinFactorization_H
00012 #define CoinFactorization_H
00013 //#define COIN_ONE_ETA_COPY 100
00014 
00015 #include <iostream>
00016 #include <string>
00017 #include <cassert>
00018 #include <cstdio>
00019 #include "CoinFinite.hpp"
00020 #include "CoinIndexedVector.hpp"
00021 class CoinPackedMatrix;
00047 class CoinFactorization {
00048    friend void CoinFactorizationUnitTest( const std::string & mpsDir );
00049 
00050 public:
00051 
00054 
00055     CoinFactorization (  );
00057   CoinFactorization ( const CoinFactorization &other);
00058 
00060    ~CoinFactorization (  );
00062   void almostDestructor();
00064   void show_self (  ) const;
00066   int saveFactorization (const char * file  ) const;
00070   int restoreFactorization (const char * file  , bool factor=false) ;
00072   void sort (  ) const;
00074     CoinFactorization & operator = ( const CoinFactorization & other );
00076 
00086   int factorize ( const CoinPackedMatrix & matrix, 
00087                   int rowIsBasic[], int columnIsBasic[] , 
00088                   double areaFactor = 0.0 );
00099   int factorize ( int numberRows,
00100                   int numberColumns,
00101                   CoinBigIndex numberElements,
00102                   CoinBigIndex maximumL,
00103                   CoinBigIndex maximumU,
00104                   const int indicesRow[],
00105                   const int indicesColumn[], const double elements[] ,
00106                   int permutation[],
00107                   double areaFactor = 0.0);
00112   int factorizePart1 ( int numberRows,
00113                        int numberColumns,
00114                        CoinBigIndex estimateNumberElements,
00115                        int * indicesRow[],
00116                        int * indicesColumn[],
00117                        CoinFactorizationDouble * elements[],
00118                   double areaFactor = 0.0);
00125   int factorizePart2 (int permutation[],int exactNumberElements);
00127   double conditionNumber() const;
00128   
00130 
00133 
00134   inline int status (  ) const {
00135     return status_;
00136   }
00138   inline void setStatus (  int value)
00139   {  status_=value;  }
00141   inline int pivots (  ) const {
00142     return numberPivots_;
00143   }
00145   inline void setPivots (  int value ) 
00146   { numberPivots_=value; }
00148   inline int *permute (  ) const {
00149     return permute_.array();
00150   }
00152   inline int *pivotColumn (  ) const {
00153     return pivotColumn_.array();
00154   }
00156   inline CoinFactorizationDouble *pivotRegion (  ) const {
00157     return pivotRegion_.array();
00158   }
00160   inline int *permuteBack (  ) const {
00161     return permuteBack_.array();
00162   }
00165   inline int *pivotColumnBack (  ) const {
00166     //return firstCount_.array();
00167     return pivotColumnBack_.array();
00168   }
00170   inline CoinBigIndex * startRowL() const
00171   { return startRowL_.array();}
00172 
00174   inline CoinBigIndex * startColumnL() const
00175   { return startColumnL_.array();}
00176 
00178   inline int * indexColumnL() const
00179   { return indexColumnL_.array();}
00180 
00182   inline int * indexRowL() const
00183   { return indexRowL_.array();}
00184 
00186   inline CoinFactorizationDouble * elementByRowL() const
00187   { return elementByRowL_.array();}
00188 
00190   inline int numberRowsExtra (  ) const {
00191     return numberRowsExtra_;
00192   }
00194   inline void setNumberRows(int value)
00195   { numberRows_ = value; }
00197   inline int numberRows (  ) const {
00198     return numberRows_;
00199   }
00201   inline CoinBigIndex numberL() const
00202   { return numberL_;}
00203 
00205   inline CoinBigIndex baseL() const
00206   { return baseL_;}
00208   inline int maximumRowsExtra (  ) const {
00209     return maximumRowsExtra_;
00210   }
00212   inline int numberColumns (  ) const {
00213     return numberColumns_;
00214   }
00216   inline int numberElements (  ) const {
00217     return totalElements_;
00218   }
00220   inline int numberForrestTomlin (  ) const {
00221     return numberInColumn_.array()[numberColumnsExtra_];
00222   }
00224   inline int numberGoodColumns (  ) const {
00225     return numberGoodU_;
00226   }
00228   inline double areaFactor (  ) const {
00229     return areaFactor_;
00230   }
00231   inline void areaFactor ( double value ) {
00232     areaFactor_=value;
00233   }
00235   double adjustedAreaFactor() const;
00237   inline void relaxAccuracyCheck(double value)
00238   { relaxCheck_ = value;}
00239   inline double getAccuracyCheck() const
00240   { return relaxCheck_;}
00242   inline int messageLevel (  ) const {
00243     return messageLevel_ ;
00244   }
00245   void messageLevel (  int value );
00247   inline int maximumPivots (  ) const {
00248     return maximumPivots_ ;
00249   }
00250   void maximumPivots (  int value );
00251 
00253   inline int denseThreshold() const 
00254     { return denseThreshold_;}
00256   inline void setDenseThreshold(int value)
00257     { denseThreshold_ = value;}
00259   inline double pivotTolerance (  ) const {
00260     return pivotTolerance_ ;
00261   }
00262   void pivotTolerance (  double value );
00264   inline double zeroTolerance (  ) const {
00265     return zeroTolerance_ ;
00266   }
00267   void zeroTolerance (  double value );
00268 #ifndef COIN_FAST_CODE
00270   inline double slackValue (  ) const {
00271     return slackValue_ ;
00272   }
00273   void slackValue (  double value );
00274 #endif
00276   double maximumCoefficient() const;
00278   inline bool forrestTomlin() const
00279   { return doForrestTomlin_;}
00280   inline void setForrestTomlin(bool value)
00281   { doForrestTomlin_=value;}
00283   inline bool spaceForForrestTomlin() const
00284   {
00285     CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
00286     CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
00287     return (space>=0)&&doForrestTomlin_;
00288   }
00290 
00293 
00295   inline int numberDense() const
00296   { return numberDense_;}
00297 
00299   inline CoinBigIndex numberElementsU (  ) const {
00300     return lengthU_;
00301   }
00303   inline void setNumberElementsU(CoinBigIndex value)
00304   { lengthU_ = value; }
00306   inline CoinBigIndex lengthAreaU (  ) const {
00307     return lengthAreaU_;
00308   }
00310   inline CoinBigIndex numberElementsL (  ) const {
00311     return lengthL_;
00312   }
00314   inline CoinBigIndex lengthAreaL (  ) const {
00315     return lengthAreaL_;
00316   }
00318   inline CoinBigIndex numberElementsR (  ) const {
00319     return lengthR_;
00320   }
00322   inline CoinBigIndex numberCompressions() const
00323   { return numberCompressions_;}
00325   inline int * numberInRow() const
00326   { return numberInRow_.array();}
00328   inline int * numberInColumn() const
00329   { return numberInColumn_.array();}
00331   inline CoinFactorizationDouble * elementU() const
00332   { return elementU_.array();}
00334   inline int * indexRowU() const
00335   { return indexRowU_.array();}
00337   inline CoinBigIndex * startColumnU() const
00338   { return startColumnU_.array();}
00340   inline int maximumColumnsExtra()
00341   { return maximumColumnsExtra_;}
00345   inline int biasLU() const
00346   { return biasLU_;}
00347   inline void setBiasLU(int value)
00348   { biasLU_=value;}
00354   inline int persistenceFlag() const
00355   { return persistenceFlag_;}
00356   void setPersistenceFlag(int value);
00358 
00361 
00369   int replaceColumn ( CoinIndexedVector * regionSparse,
00370                       int pivotRow,
00371                       double pivotCheck ,
00372                       bool checkBeforeModifying=false,
00373                       double acceptablePivot=1.0e-8);
00378   void replaceColumnU ( CoinIndexedVector * regionSparse,
00379                         CoinBigIndex * deleted,
00380                         int internalPivotRow);
00382 
00392   int updateColumnFT ( CoinIndexedVector * regionSparse,
00393                        CoinIndexedVector * regionSparse2);
00396   int updateColumn ( CoinIndexedVector * regionSparse,
00397                      CoinIndexedVector * regionSparse2,
00398                      bool noPermute=false) const;
00404   int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
00405                            CoinIndexedVector * regionSparse2,
00406                            CoinIndexedVector * regionSparse3,
00407                            bool noPermuteRegion3=false) ;
00412   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00413                               CoinIndexedVector * regionSparse2) const;
00415   void goSparse();
00417   inline int sparseThreshold ( ) const
00418   { return sparseThreshold_;}
00420   void sparseThreshold ( int value );
00422 
00423 
00427 
00428   inline void clearArrays()
00429   { gutsOfDestructor();}
00431 
00434 
00437   int add ( CoinBigIndex numberElements,
00438                int indicesRow[],
00439                int indicesColumn[], double elements[] );
00440 
00443   int addColumn ( CoinBigIndex numberElements,
00444                      int indicesRow[], double elements[] );
00445 
00448   int addRow ( CoinBigIndex numberElements,
00449                   int indicesColumn[], double elements[] );
00450 
00452   int deleteColumn ( int Row );
00454   int deleteRow ( int Row );
00455 
00459   int replaceRow ( int whichRow, int numberElements,
00460                       const int indicesColumn[], const double elements[] );
00462   void emptyRows(int numberToEmpty, const int which[]);
00464 
00465 
00466   void checkSparse();
00468   inline bool collectStatistics() const
00469   { return collectStatistics_;}
00471   inline void setCollectStatistics(bool onOff) const
00472   { collectStatistics_ = onOff;}
00474   void gutsOfDestructor(int type=1);
00476   void gutsOfInitialize(int type);
00477   void gutsOfCopy(const CoinFactorization &other);
00478 
00480   void resetStatistics();
00481 
00482 
00484 
00486 
00487   void getAreas ( int numberRows,
00488                   int numberColumns,
00489                   CoinBigIndex maximumL,
00490                   CoinBigIndex maximumU );
00491 
00494   void preProcess ( int state,
00495                     int possibleDuplicates = -1 );
00497   int factor (  );
00498 protected:
00501   int factorSparse (  );
00504   int factorSparseSmall (  );
00507   int factorSparseLarge (  );
00510   int factorDense (  );
00511 
00513   bool pivotOneOtherRow ( int pivotRow,
00514                           int pivotColumn );
00516   bool pivotRowSingleton ( int pivotRow,
00517                            int pivotColumn );
00519   bool pivotColumnSingleton ( int pivotRow,
00520                               int pivotColumn );
00521 
00526   bool getColumnSpace ( int iColumn,
00527                         int extraNeeded );
00528 
00531   bool reorderU();
00535   bool getColumnSpaceIterateR ( int iColumn, double value,
00536                                int iRow);
00542   CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00543                                int iRow);
00547   bool getRowSpace ( int iRow, int extraNeeded );
00548 
00552   bool getRowSpaceIterate ( int iRow,
00553                             int extraNeeded );
00555   void checkConsistency (  );
00557   inline void addLink ( int index, int count ) {
00558     int *nextCount = nextCount_.array();
00559     int *firstCount = firstCount_.array();
00560     int *lastCount = lastCount_.array();
00561     int next = firstCount[count];
00562       lastCount[index] = -2 - count;
00563     if ( next < 0 ) {
00564       //first with that count
00565       firstCount[count] = index;
00566       nextCount[index] = -1;
00567     } else {
00568       firstCount[count] = index;
00569       nextCount[index] = next;
00570       lastCount[next] = index;
00571   }}
00573   inline void deleteLink ( int index ) {
00574     int *nextCount = nextCount_.array();
00575     int *firstCount = firstCount_.array();
00576     int *lastCount = lastCount_.array();
00577     int next = nextCount[index];
00578     int last = lastCount[index];
00579     if ( last >= 0 ) {
00580       nextCount[last] = next;
00581     } else {
00582       int count = -last - 2;
00583 
00584       firstCount[count] = next;
00585     }
00586     if ( next >= 0 ) {
00587       lastCount[next] = last;
00588     }
00589     nextCount[index] = -2;
00590     lastCount[index] = -2;
00591     return;
00592   }
00594   void separateLinks(int count,bool rowsFirst);
00596   void cleanup (  );
00597 
00599   void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00601   void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00603   void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00605   void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00606 
00608   void updateColumnR ( CoinIndexedVector * region ) const;
00611   void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00612 
00614   void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00615 
00617   void updateColumnUSparse ( CoinIndexedVector * regionSparse, 
00618                              int * indexIn) const;
00620   void updateColumnUSparsish ( CoinIndexedVector * regionSparse, 
00621                                int * indexIn) const;
00623   int updateColumnUDensish ( double * COIN_RESTRICT region, 
00624                              int * COIN_RESTRICT regionIndex) const;
00626   void updateTwoColumnsUDensish (
00627                                  int & numberNonZero1,
00628                                  double * COIN_RESTRICT region1, 
00629                                  int * COIN_RESTRICT index1,
00630                                  int & numberNonZero2,
00631                                  double * COIN_RESTRICT region2, 
00632                                  int * COIN_RESTRICT index2) const;
00634   void updateColumnPFI ( CoinIndexedVector * regionSparse) const; 
00636   void permuteBack ( CoinIndexedVector * regionSparse, 
00637                      CoinIndexedVector * outVector) const;
00638 
00640   void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00643   void updateColumnTransposeU ( CoinIndexedVector * region,
00644                                 int smallestIndex) const;
00647   void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00648                                         int smallestIndex) const;
00651   void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00652                                        int smallestIndex) const;
00655   void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00658   void updateColumnTransposeUByColumn ( CoinIndexedVector * region,
00659                                         int smallestIndex) const;
00660 
00662   void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00664   void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00666   void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00667 
00669   void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00671   void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00673   void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00675   void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00677   void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00678 public:
00683   int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00684                          int pivotRow, double alpha);
00685 protected:
00688   int checkPivot(double saveFromU, double oldPivot) const;
00689   /********************************* START LARGE TEMPLATE ********/
00690 #ifdef INT_IS_8
00691 #define COINFACTORIZATION_BITS_PER_INT 64
00692 #define COINFACTORIZATION_SHIFT_PER_INT 6
00693 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00694 #else
00695 #define COINFACTORIZATION_BITS_PER_INT 32
00696 #define COINFACTORIZATION_SHIFT_PER_INT 5
00697 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00698 #endif
00699   template <class T>  inline bool
00700   pivot ( int pivotRow,
00701           int pivotColumn,
00702           CoinBigIndex pivotRowPosition,
00703           CoinBigIndex pivotColumnPosition,
00704           CoinFactorizationDouble work[],
00705           unsigned int workArea2[],
00706           int increment2,
00707           T markRow[] ,
00708           int largeInteger)
00709 {
00710   int *indexColumnU = indexColumnU_.array();
00711   CoinBigIndex *startColumnU = startColumnU_.array();
00712   int *numberInColumn = numberInColumn_.array();
00713   CoinFactorizationDouble *elementU = elementU_.array();
00714   int *indexRowU = indexRowU_.array();
00715   CoinBigIndex *startRowU = startRowU_.array();
00716   int *numberInRow = numberInRow_.array();
00717   CoinFactorizationDouble *elementL = elementL_.array();
00718   int *indexRowL = indexRowL_.array();
00719   int *saveColumn = saveColumn_.array();
00720   int *nextRow = nextRow_.array();
00721   int *lastRow = lastRow_.array() ;
00722 
00723   //store pivot columns (so can easily compress)
00724   int numberInPivotRow = numberInRow[pivotRow] - 1;
00725   CoinBigIndex startColumn = startColumnU[pivotColumn];
00726   int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00727   CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00728   int put = 0;
00729   CoinBigIndex startRow = startRowU[pivotRow];
00730   CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00731 
00732   if ( pivotColumnPosition < 0 ) {
00733     for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00734       int iColumn = indexColumnU[pivotColumnPosition];
00735       if ( iColumn != pivotColumn ) {
00736         saveColumn[put++] = iColumn;
00737       } else {
00738         break;
00739       }
00740     }
00741   } else {
00742     for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00743       saveColumn[put++] = indexColumnU[i];
00744     }
00745   }
00746   assert (pivotColumnPosition<endRow);
00747   assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00748   pivotColumnPosition++;
00749   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00750     saveColumn[put++] = indexColumnU[pivotColumnPosition];
00751   }
00752   //take out this bit of indexColumnU
00753   int next = nextRow[pivotRow];
00754   int last = lastRow[pivotRow];
00755 
00756   nextRow[last] = next;
00757   lastRow[next] = last;
00758   nextRow[pivotRow] = numberGoodU_;     //use for permute
00759   lastRow[pivotRow] = -2;
00760   numberInRow[pivotRow] = 0;
00761   //store column in L, compress in U and take column out
00762   CoinBigIndex l = lengthL_;
00763 
00764   if ( l + numberInPivotColumn > lengthAreaL_ ) {
00765     //need more memory
00766     if ((messageLevel_&4)!=0) 
00767       printf("more memory needed in middle of invert\n");
00768     return false;
00769   }
00770   //l+=currentAreaL_->elementByColumn-elementL;
00771   CoinBigIndex lSave = l;
00772 
00773   CoinBigIndex * startColumnL = startColumnL_.array();
00774   startColumnL[numberGoodL_] = l;       //for luck and first time
00775   numberGoodL_++;
00776   startColumnL[numberGoodL_] = l + numberInPivotColumn;
00777   lengthL_ += numberInPivotColumn;
00778   if ( pivotRowPosition < 0 ) {
00779     for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00780       int iRow = indexRowU[pivotRowPosition];
00781       if ( iRow != pivotRow ) {
00782         indexRowL[l] = iRow;
00783         elementL[l] = elementU[pivotRowPosition];
00784         markRow[iRow] = static_cast<T>(l - lSave);
00785         l++;
00786         //take out of row list
00787         CoinBigIndex start = startRowU[iRow];
00788         CoinBigIndex end = start + numberInRow[iRow];
00789         CoinBigIndex where = start;
00790 
00791         while ( indexColumnU[where] != pivotColumn ) {
00792           where++;
00793         }                       /* endwhile */
00794 #if DEBUG_COIN
00795         if ( where >= end ) {
00796           abort (  );
00797         }
00798 #endif
00799         indexColumnU[where] = indexColumnU[end - 1];
00800         numberInRow[iRow]--;
00801       } else {
00802         break;
00803       }
00804     }
00805   } else {
00806     CoinBigIndex i;
00807 
00808     for ( i = startColumn; i < pivotRowPosition; i++ ) {
00809       int iRow = indexRowU[i];
00810 
00811       markRow[iRow] = static_cast<T>(l - lSave);
00812       indexRowL[l] = iRow;
00813       elementL[l] = elementU[i];
00814       l++;
00815       //take out of row list
00816       CoinBigIndex start = startRowU[iRow];
00817       CoinBigIndex end = start + numberInRow[iRow];
00818       CoinBigIndex where = start;
00819 
00820       while ( indexColumnU[where] != pivotColumn ) {
00821         where++;
00822       }                         /* endwhile */
00823 #if DEBUG_COIN
00824       if ( where >= end ) {
00825         abort (  );
00826       }
00827 #endif
00828       indexColumnU[where] = indexColumnU[end - 1];
00829       numberInRow[iRow]--;
00830       assert (numberInRow[iRow]>=0);
00831     }
00832   }
00833   assert (pivotRowPosition<endColumn);
00834   assert (indexRowU[pivotRowPosition]==pivotRow);
00835   CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
00836   CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
00837 
00838   pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00839   pivotRowPosition++;
00840   for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00841     int iRow = indexRowU[pivotRowPosition];
00842     
00843     markRow[iRow] = static_cast<T>(l - lSave);
00844     indexRowL[l] = iRow;
00845     elementL[l] = elementU[pivotRowPosition];
00846     l++;
00847     //take out of row list
00848     CoinBigIndex start = startRowU[iRow];
00849     CoinBigIndex end = start + numberInRow[iRow];
00850     CoinBigIndex where = start;
00851     
00852     while ( indexColumnU[where] != pivotColumn ) {
00853       where++;
00854     }                           /* endwhile */
00855 #if DEBUG_COIN
00856     if ( where >= end ) {
00857       abort (  );
00858     }
00859 #endif
00860     indexColumnU[where] = indexColumnU[end - 1];
00861     numberInRow[iRow]--;
00862     assert (numberInRow[iRow]>=0);
00863   }
00864   markRow[pivotRow] = static_cast<T>(largeInteger);
00865   //compress pivot column (move pivot to front including saved)
00866   numberInColumn[pivotColumn] = 0;
00867   //use end of L for temporary space
00868   int *indexL = &indexRowL[lSave];
00869   CoinFactorizationDouble *multipliersL = &elementL[lSave];
00870 
00871   //adjust
00872   int j;
00873 
00874   for ( j = 0; j < numberInPivotColumn; j++ ) {
00875     multipliersL[j] *= pivotMultiplier;
00876   }
00877   //zero out fill
00878   CoinBigIndex iErase;
00879   for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00880         iErase++ ) {
00881     workArea2[iErase] = 0;
00882   }
00883   CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00884   unsigned int *temp2 = workArea2;
00885   int * nextColumn = nextColumn_.array();
00886 
00887   //pack down and move to work
00888   int jColumn;
00889   for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
00890     int iColumn = saveColumn[jColumn];
00891     CoinBigIndex startColumn = startColumnU[iColumn];
00892     CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
00893     int iRow = indexRowU[startColumn];
00894     CoinFactorizationDouble value = elementU[startColumn];
00895     double largest;
00896     CoinBigIndex put = startColumn;
00897     CoinBigIndex positionLargest = -1;
00898     CoinFactorizationDouble thisPivotValue = 0.0;
00899 
00900     //compress column and find largest not updated
00901     bool checkLargest;
00902     int mark = markRow[iRow];
00903 
00904     if ( mark == largeInteger+1 ) {
00905       largest = fabs ( value );
00906       positionLargest = put;
00907       put++;
00908       checkLargest = false;
00909     } else {
00910       //need to find largest
00911       largest = 0.0;
00912       checkLargest = true;
00913       if ( mark != largeInteger ) {
00914         //will be updated
00915         work[mark] = value;
00916         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00917         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00918 
00919         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
00920         added--;
00921       } else {
00922         thisPivotValue = value;
00923       }
00924     }
00925     CoinBigIndex i;
00926     for ( i = startColumn + 1; i < endColumn; i++ ) {
00927       iRow = indexRowU[i];
00928       value = elementU[i];
00929       int mark = markRow[iRow];
00930 
00931       if ( mark == largeInteger+1 ) {
00932         //keep
00933         indexRowU[put] = iRow;
00934         elementU[put] = value;
00935         if ( checkLargest ) {
00936           double absValue = fabs ( value );
00937 
00938           if ( absValue > largest ) {
00939             largest = absValue;
00940             positionLargest = put;
00941           }
00942         }
00943         put++;
00944       } else if ( mark != largeInteger ) {
00945         //will be updated
00946         work[mark] = value;
00947         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00948         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00949 
00950         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
00951         added--;
00952       } else {
00953         thisPivotValue = value;
00954       }
00955     }
00956     //slot in pivot
00957     elementU[put] = elementU[startColumn];
00958     indexRowU[put] = indexRowU[startColumn];
00959     if ( positionLargest == startColumn ) {
00960       positionLargest = put;    //follow if was largest
00961     }
00962     put++;
00963     elementU[startColumn] = thisPivotValue;
00964     indexRowU[startColumn] = pivotRow;
00965     //clean up counts
00966     startColumn++;
00967     numberInColumn[iColumn] = put - startColumn;
00968     int * numberInColumnPlus = numberInColumnPlus_.array();
00969     numberInColumnPlus[iColumn]++;
00970     startColumnU[iColumn]++;
00971     //how much space have we got
00972     int next = nextColumn[iColumn];
00973     CoinBigIndex space;
00974 
00975     space = startColumnU[next] - put - numberInColumnPlus[next];
00976     //assume no zero elements
00977     if ( numberInPivotColumn > space ) {
00978       //getColumnSpace also moves fixed part
00979       if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
00980         return false;
00981       }
00982       //redo starts
00983       positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
00984       startColumn = startColumnU[iColumn];
00985       put = startColumn + numberInColumn[iColumn];
00986     }
00987     double tolerance = zeroTolerance_;
00988 
00989     int *nextCount = nextCount_.array();
00990     for ( j = 0; j < numberInPivotColumn; j++ ) {
00991       value = work[j] - thisPivotValue * multipliersL[j];
00992       double absValue = fabs ( value );
00993 
00994       if ( absValue > tolerance ) {
00995         work[j] = 0.0;
00996         elementU[put] = value;
00997         indexRowU[put] = indexL[j];
00998         if ( absValue > largest ) {
00999           largest = absValue;
01000           positionLargest = put;
01001         }
01002         put++;
01003       } else {
01004         work[j] = 0.0;
01005         added--;
01006         int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01007         int bit = j & COINFACTORIZATION_MASK_PER_INT;
01008 
01009         if ( temp2[word] & ( 1 << bit ) ) {
01010           //take out of row list
01011           iRow = indexL[j];
01012           CoinBigIndex start = startRowU[iRow];
01013           CoinBigIndex end = start + numberInRow[iRow];
01014           CoinBigIndex where = start;
01015 
01016           while ( indexColumnU[where] != iColumn ) {
01017             where++;
01018           }                     /* endwhile */
01019 #if DEBUG_COIN
01020           if ( where >= end ) {
01021             abort (  );
01022           }
01023 #endif
01024           indexColumnU[where] = indexColumnU[end - 1];
01025           numberInRow[iRow]--;
01026         } else {
01027           //make sure won't be added
01028           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01029           int bit = j & COINFACTORIZATION_MASK_PER_INT;
01030 
01031           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
01032         }
01033       }
01034     }
01035     numberInColumn[iColumn] = put - startColumn;
01036     //move largest
01037     if ( positionLargest >= 0 ) {
01038       value = elementU[positionLargest];
01039       iRow = indexRowU[positionLargest];
01040       elementU[positionLargest] = elementU[startColumn];
01041       indexRowU[positionLargest] = indexRowU[startColumn];
01042       elementU[startColumn] = value;
01043       indexRowU[startColumn] = iRow;
01044     }
01045     //linked list for column
01046     if ( nextCount[iColumn + numberRows_] != -2 ) {
01047       //modify linked list
01048       deleteLink ( iColumn + numberRows_ );
01049       addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01050     }
01051     temp2 += increment2;
01052   }
01053   //get space for row list
01054   unsigned int *putBase = workArea2;
01055   int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01056   int i = 0;
01057 
01058   // do linked lists and update counts
01059   while ( bigLoops ) {
01060     bigLoops--;
01061     int bit;
01062     for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01063       unsigned int *putThis = putBase;
01064       int iRow = indexL[i];
01065 
01066       //get space
01067       int number = 0;
01068       int jColumn;
01069 
01070       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01071         unsigned int test = *putThis;
01072 
01073         putThis += increment2;
01074         test = 1 - ( ( test >> bit ) & 1 );
01075         number += test;
01076       }
01077       int next = nextRow[iRow];
01078       CoinBigIndex space;
01079 
01080       space = startRowU[next] - startRowU[iRow];
01081       number += numberInRow[iRow];
01082       if ( space < number ) {
01083         if ( !getRowSpace ( iRow, number ) ) {
01084           return false;
01085         }
01086       }
01087       // now do
01088       putThis = putBase;
01089       next = nextRow[iRow];
01090       number = numberInRow[iRow];
01091       CoinBigIndex end = startRowU[iRow] + number;
01092       int saveIndex = indexColumnU[startRowU[next]];
01093 
01094       //add in
01095       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01096         unsigned int test = *putThis;
01097 
01098         putThis += increment2;
01099         test = 1 - ( ( test >> bit ) & 1 );
01100         indexColumnU[end] = saveColumn[jColumn];
01101         end += test;
01102       }
01103       //put back next one in case zapped
01104       indexColumnU[startRowU[next]] = saveIndex;
01105       markRow[iRow] = static_cast<T>(largeInteger+1);
01106       number = end - startRowU[iRow];
01107       numberInRow[iRow] = number;
01108       deleteLink ( iRow );
01109       addLink ( iRow, number );
01110     }
01111     putBase++;
01112   }                             /* endwhile */
01113   int bit;
01114 
01115   for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01116     unsigned int *putThis = putBase;
01117     int iRow = indexL[i];
01118 
01119     //get space
01120     int number = 0;
01121     int jColumn;
01122 
01123     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01124       unsigned int test = *putThis;
01125 
01126       putThis += increment2;
01127       test = 1 - ( ( test >> bit ) & 1 );
01128       number += test;
01129     }
01130     int next = nextRow[iRow];
01131     CoinBigIndex space;
01132 
01133     space = startRowU[next] - startRowU[iRow];
01134     number += numberInRow[iRow];
01135     if ( space < number ) {
01136       if ( !getRowSpace ( iRow, number ) ) {
01137         return false;
01138       }
01139     }
01140     // now do
01141     putThis = putBase;
01142     next = nextRow[iRow];
01143     number = numberInRow[iRow];
01144     CoinBigIndex end = startRowU[iRow] + number;
01145     int saveIndex;
01146 
01147     saveIndex = indexColumnU[startRowU[next]];
01148 
01149     //add in
01150     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01151       unsigned int test = *putThis;
01152 
01153       putThis += increment2;
01154       test = 1 - ( ( test >> bit ) & 1 );
01155 
01156       indexColumnU[end] = saveColumn[jColumn];
01157       end += test;
01158     }
01159     indexColumnU[startRowU[next]] = saveIndex;
01160     markRow[iRow] = static_cast<T>(largeInteger+1);
01161     number = end - startRowU[iRow];
01162     numberInRow[iRow] = number;
01163     deleteLink ( iRow );
01164     addLink ( iRow, number );
01165   }
01166   markRow[pivotRow] = static_cast<T>(largeInteger+1);
01167   //modify linked list for pivots
01168   deleteLink ( pivotRow );
01169   deleteLink ( pivotColumn + numberRows_ );
01170   totalElements_ += added;
01171   return true;
01172 }
01173 
01174   /********************************* END LARGE TEMPLATE ********/
01176 
01177 protected:
01178 
01181 
01182   double pivotTolerance_;
01184   double zeroTolerance_;
01185 #ifndef COIN_FAST_CODE
01187   double slackValue_;
01188 #else
01189 #ifndef slackValue_
01190 #define slackValue_ -1.0
01191 #endif
01192 #endif
01194   double areaFactor_;
01196   double relaxCheck_;
01198   int numberRows_;
01200   int numberRowsExtra_;
01202   int maximumRowsExtra_;
01204   int numberColumns_;
01206   int numberColumnsExtra_;
01208   int maximumColumnsExtra_;
01210   int numberGoodU_;
01212   int numberGoodL_;
01214   int maximumPivots_;
01216   int numberPivots_;
01219   CoinBigIndex totalElements_;
01221   CoinBigIndex factorElements_;
01223   CoinIntArrayWithLength pivotColumn_;
01225   CoinIntArrayWithLength permute_;
01227   CoinIntArrayWithLength permuteBack_;
01229   CoinIntArrayWithLength pivotColumnBack_;
01231   int status_;
01232 
01237   //int increasingRows_;
01238 
01240   int numberTrials_;
01242   CoinBigIndexArrayWithLength startRowU_;
01243 
01245   CoinIntArrayWithLength numberInRow_;
01246 
01248   CoinIntArrayWithLength numberInColumn_;
01249 
01251   CoinIntArrayWithLength numberInColumnPlus_;
01252 
01255   CoinIntArrayWithLength firstCount_;
01256 
01258   CoinIntArrayWithLength nextCount_;
01259 
01261   CoinIntArrayWithLength lastCount_;
01262 
01264   CoinIntArrayWithLength nextColumn_;
01265 
01267   CoinIntArrayWithLength lastColumn_;
01268 
01270   CoinIntArrayWithLength nextRow_;
01271 
01273   CoinIntArrayWithLength lastRow_;
01274 
01276   CoinIntArrayWithLength saveColumn_;
01277 
01279   CoinIntArrayWithLength markRow_;
01280 
01282   int messageLevel_;
01283 
01285   int biggerDimension_;
01286 
01288   CoinIntArrayWithLength indexColumnU_;
01289 
01291   CoinIntArrayWithLength pivotRowL_;
01292 
01294   CoinFactorizationDoubleArrayWithLength pivotRegion_;
01295 
01297   int numberSlacks_;
01298 
01300   int numberU_;
01301 
01303   CoinBigIndex maximumU_;
01304 
01306   //int baseU_;
01307 
01309   CoinBigIndex lengthU_;
01310 
01312   CoinBigIndex lengthAreaU_;
01313 
01315   CoinFactorizationDoubleArrayWithLength elementU_;
01316 
01318   CoinIntArrayWithLength indexRowU_;
01319 
01321   CoinBigIndexArrayWithLength startColumnU_;
01322 
01324   CoinBigIndexArrayWithLength convertRowToColumnU_;
01325 
01327   CoinBigIndex numberL_;
01328 
01330   CoinBigIndex baseL_;
01331 
01333   CoinBigIndex lengthL_;
01334 
01336   CoinBigIndex lengthAreaL_;
01337 
01339   CoinFactorizationDoubleArrayWithLength elementL_;
01340 
01342   CoinIntArrayWithLength indexRowL_;
01343 
01345   CoinBigIndexArrayWithLength startColumnL_;
01346 
01348   bool doForrestTomlin_;
01349 
01351   int numberR_;
01352 
01354   CoinBigIndex lengthR_;
01355 
01357   CoinBigIndex lengthAreaR_;
01358 
01360   CoinFactorizationDouble *elementR_;
01361 
01363   int *indexRowR_;
01364 
01366   CoinBigIndexArrayWithLength startColumnR_;
01367 
01369   double  * denseArea_;
01370 
01372   int * densePermute_;
01373 
01375   int numberDense_;
01376 
01378   int denseThreshold_;
01379 
01381   CoinFactorizationDoubleArrayWithLength workArea_;
01382 
01384   CoinUnsignedIntArrayWithLength workArea2_;
01385 
01387   CoinBigIndex numberCompressions_;
01388 
01390   mutable double ftranCountInput_;
01391   mutable double ftranCountAfterL_;
01392   mutable double ftranCountAfterR_;
01393   mutable double ftranCountAfterU_;
01394   mutable double btranCountInput_;
01395   mutable double btranCountAfterU_;
01396   mutable double btranCountAfterR_;
01397   mutable double btranCountAfterL_;
01398 
01400   mutable int numberFtranCounts_;
01401   mutable int numberBtranCounts_;
01402 
01404   double ftranAverageAfterL_;
01405   double ftranAverageAfterR_;
01406   double ftranAverageAfterU_;
01407   double btranAverageAfterU_;
01408   double btranAverageAfterR_;
01409   double btranAverageAfterL_;
01410 
01412   mutable bool collectStatistics_;
01413 
01415   int sparseThreshold_;
01416 
01418   int sparseThreshold2_;
01419 
01421   CoinBigIndexArrayWithLength startRowL_;
01422 
01424   CoinIntArrayWithLength indexColumnL_;
01425 
01427   CoinFactorizationDoubleArrayWithLength elementByRowL_;
01428 
01430   mutable CoinIntArrayWithLength sparse_;
01434   int biasLU_;
01440   int persistenceFlag_;
01442 };
01443 // Dense coding
01444 #ifdef COIN_HAS_LAPACK
01445 #define DENSE_CODE 1
01446 /* Type of Fortran integer translated into C */
01447 #ifndef ipfint
01448 //typedef ipfint FORTRAN_INTEGER_TYPE ;
01449 typedef int ipfint;
01450 typedef const int cipfint;
01451 #endif
01452 #endif
01453 #endif
01454 // Extra for ugly include
01455 #ifdef UGLY_COIN_FACTOR_CODING
01456 #define FAC_UNSET (FAC_SET+1)
01457 {
01458   goodPivot=false;
01459   //store pivot columns (so can easily compress)
01460   CoinBigIndex startColumnThis = startColumn[iPivotColumn];
01461   CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
01462   int put = 0;
01463   CoinBigIndex startRowThis = startRow[iPivotRow];
01464   CoinBigIndex endRow = startRowThis + numberDoRow + 1;
01465   if ( pivotColumnPosition < 0 ) {
01466     for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01467       int iColumn = indexColumn[pivotColumnPosition];
01468       if ( iColumn != iPivotColumn ) {
01469         saveColumn[put++] = iColumn;
01470       } else {
01471         break;
01472       }
01473     }
01474   } else {
01475     for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
01476       saveColumn[put++] = indexColumn[i];
01477     }
01478   }
01479   assert (pivotColumnPosition<endRow);
01480   assert (indexColumn[pivotColumnPosition]==iPivotColumn);
01481   pivotColumnPosition++;
01482   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01483     saveColumn[put++] = indexColumn[pivotColumnPosition];
01484   }
01485   //take out this bit of indexColumn
01486   int next = nextRow[iPivotRow];
01487   int last = lastRow[iPivotRow];
01488   
01489   nextRow[last] = next;
01490   lastRow[next] = last;
01491   nextRow[iPivotRow] = numberGoodU_;    //use for permute
01492   lastRow[iPivotRow] = -2;
01493   numberInRow[iPivotRow] = 0;
01494   //store column in L, compress in U and take column out
01495   CoinBigIndex l = lengthL_;
01496   // **** HORRID coding coming up but a goto seems best!
01497   {
01498     if ( l + numberDoColumn > lengthAreaL_ ) {
01499       //need more memory
01500       if ((messageLevel_&4)!=0) 
01501         printf("more memory needed in middle of invert\n");
01502       goto BAD_PIVOT;
01503     }
01504     //l+=currentAreaL_->elementByColumn-elementL;
01505     CoinBigIndex lSave = l;
01506     
01507     CoinBigIndex * startColumnL = startColumnL_.array();
01508     startColumnL[numberGoodL_] = l;     //for luck and first time
01509     numberGoodL_++;
01510     startColumnL[numberGoodL_] = l + numberDoColumn;
01511     lengthL_ += numberDoColumn;
01512     if ( pivotRowPosition < 0 ) {
01513       for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01514         int iRow = indexRow[pivotRowPosition];
01515         if ( iRow != iPivotRow ) {
01516           indexRowL[l] = iRow;
01517           elementL[l] = element[pivotRowPosition];
01518           markRow[iRow] = l - lSave;
01519           l++;
01520           //take out of row list
01521           CoinBigIndex start = startRow[iRow];
01522           CoinBigIndex end = start + numberInRow[iRow];
01523           CoinBigIndex where = start;
01524           
01525           while ( indexColumn[where] != iPivotColumn ) {
01526             where++;
01527           }                     /* endwhile */
01528 #if DEBUG_COIN
01529           if ( where >= end ) {
01530             abort (  );
01531           }
01532 #endif
01533           indexColumn[where] = indexColumn[end - 1];
01534           numberInRow[iRow]--;
01535         } else {
01536           break;
01537         }
01538       }
01539     } else {
01540       CoinBigIndex i;
01541       
01542       for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
01543         int iRow = indexRow[i];
01544         
01545         markRow[iRow] = l - lSave;
01546         indexRowL[l] = iRow;
01547         elementL[l] = element[i];
01548         l++;
01549         //take out of row list
01550         CoinBigIndex start = startRow[iRow];
01551         CoinBigIndex end = start + numberInRow[iRow];
01552         CoinBigIndex where = start;
01553         
01554         while ( indexColumn[where] != iPivotColumn ) {
01555           where++;
01556         }                               /* endwhile */
01557 #if DEBUG_COIN
01558         if ( where >= end ) {
01559           abort (  );
01560         }
01561 #endif
01562         indexColumn[where] = indexColumn[end - 1];
01563         numberInRow[iRow]--;
01564         assert (numberInRow[iRow]>=0);
01565       }
01566     }
01567     assert (pivotRowPosition<endColumn);
01568     assert (indexRow[pivotRowPosition]==iPivotRow);
01569     CoinFactorizationDouble pivotElement = element[pivotRowPosition];
01570     CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
01571     
01572     pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
01573     pivotRowPosition++;
01574     for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01575       int iRow = indexRow[pivotRowPosition];
01576       
01577       markRow[iRow] = l - lSave;
01578       indexRowL[l] = iRow;
01579       elementL[l] = element[pivotRowPosition];
01580       l++;
01581       //take out of row list
01582       CoinBigIndex start = startRow[iRow];
01583       CoinBigIndex end = start + numberInRow[iRow];
01584       CoinBigIndex where = start;
01585       
01586       while ( indexColumn[where] != iPivotColumn ) {
01587         where++;
01588       }                         /* endwhile */
01589 #if DEBUG_COIN
01590       if ( where >= end ) {
01591         abort (  );
01592       }
01593 #endif
01594       indexColumn[where] = indexColumn[end - 1];
01595       numberInRow[iRow]--;
01596       assert (numberInRow[iRow]>=0);
01597     }
01598     markRow[iPivotRow] = FAC_SET;
01599     //compress pivot column (move pivot to front including saved)
01600     numberInColumn[iPivotColumn] = 0;
01601     //use end of L for temporary space
01602     int *indexL = &indexRowL[lSave];
01603     CoinFactorizationDouble *multipliersL = &elementL[lSave];
01604     
01605     //adjust
01606     int j;
01607     
01608     for ( j = 0; j < numberDoColumn; j++ ) {
01609       multipliersL[j] *= pivotMultiplier;
01610     }
01611     //zero out fill
01612     CoinBigIndex iErase;
01613     for ( iErase = 0; iErase < increment2 * numberDoRow;
01614           iErase++ ) {
01615       workArea2[iErase] = 0;
01616     }
01617     CoinBigIndex added = numberDoRow * numberDoColumn;
01618     unsigned int *temp2 = workArea2;
01619     int * nextColumn = nextColumn_.array();
01620     
01621     //pack down and move to work
01622     int jColumn;
01623     for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01624       int iColumn = saveColumn[jColumn];
01625       CoinBigIndex startColumnThis = startColumn[iColumn];
01626       CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
01627       int iRow = indexRow[startColumnThis];
01628       CoinFactorizationDouble value = element[startColumnThis];
01629       double largest;
01630       CoinBigIndex put = startColumnThis;
01631       CoinBigIndex positionLargest = -1;
01632       CoinFactorizationDouble thisPivotValue = 0.0;
01633       
01634       //compress column and find largest not updated
01635       bool checkLargest;
01636       int mark = markRow[iRow];
01637       
01638       if ( mark == FAC_UNSET ) {
01639         largest = fabs ( value );
01640         positionLargest = put;
01641         put++;
01642         checkLargest = false;
01643       } else {
01644         //need to find largest
01645         largest = 0.0;
01646         checkLargest = true;
01647         if ( mark != FAC_SET ) {
01648           //will be updated
01649           workArea[mark] = value;
01650           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01651           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01652           
01653           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
01654           added--;
01655         } else {
01656           thisPivotValue = value;
01657         }
01658       }
01659       CoinBigIndex i;
01660       for ( i = startColumnThis + 1; i < endColumn; i++ ) {
01661         iRow = indexRow[i];
01662         value = element[i];
01663         int mark = markRow[iRow];
01664         
01665         if ( mark == FAC_UNSET ) {
01666           //keep
01667           indexRow[put] = iRow;
01668           element[put] = value;
01669           if ( checkLargest ) {
01670             double absValue = fabs ( value );
01671             
01672             if ( absValue > largest ) {
01673               largest = absValue;
01674               positionLargest = put;
01675             }
01676           }
01677           put++;
01678         } else if ( mark != FAC_SET ) {
01679           //will be updated
01680           workArea[mark] = value;
01681           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01682           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01683           
01684           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
01685           added--;
01686         } else {
01687           thisPivotValue = value;
01688         }
01689       }
01690       //slot in pivot
01691       element[put] = element[startColumnThis];
01692       indexRow[put] = indexRow[startColumnThis];
01693       if ( positionLargest == startColumnThis ) {
01694         positionLargest = put;  //follow if was largest
01695       }
01696       put++;
01697       element[startColumnThis] = thisPivotValue;
01698       indexRow[startColumnThis] = iPivotRow;
01699       //clean up counts
01700       startColumnThis++;
01701       numberInColumn[iColumn] = put - startColumnThis;
01702       int * numberInColumnPlus = numberInColumnPlus_.array();
01703       numberInColumnPlus[iColumn]++;
01704       startColumn[iColumn]++;
01705       //how much space have we got
01706       int next = nextColumn[iColumn];
01707       CoinBigIndex space;
01708       
01709       space = startColumn[next] - put - numberInColumnPlus[next];
01710       //assume no zero elements
01711       if ( numberDoColumn > space ) {
01712         //getColumnSpace also moves fixed part
01713         if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
01714           goto BAD_PIVOT;
01715         }
01716         //redo starts
01717         positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
01718         startColumnThis = startColumn[iColumn];
01719         put = startColumnThis + numberInColumn[iColumn];
01720       }
01721       double tolerance = zeroTolerance_;
01722       
01723       int *nextCount = nextCount_.array();
01724       for ( j = 0; j < numberDoColumn; j++ ) {
01725         value = workArea[j] - thisPivotValue * multipliersL[j];
01726         double absValue = fabs ( value );
01727         
01728         if ( absValue > tolerance ) {
01729           workArea[j] = 0.0;
01730           element[put] = value;
01731           indexRow[put] = indexL[j];
01732           if ( absValue > largest ) {
01733             largest = absValue;
01734             positionLargest = put;
01735           }
01736           put++;
01737         } else {
01738           workArea[j] = 0.0;
01739           added--;
01740           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01741           int bit = j & COINFACTORIZATION_MASK_PER_INT;
01742           
01743           if ( temp2[word] & ( 1 << bit ) ) {
01744             //take out of row list
01745             iRow = indexL[j];
01746             CoinBigIndex start = startRow[iRow];
01747             CoinBigIndex end = start + numberInRow[iRow];
01748             CoinBigIndex where = start;
01749             
01750             while ( indexColumn[where] != iColumn ) {
01751               where++;
01752             }                   /* endwhile */
01753 #if DEBUG_COIN
01754             if ( where >= end ) {
01755               abort (  );
01756             }
01757 #endif
01758             indexColumn[where] = indexColumn[end - 1];
01759             numberInRow[iRow]--;
01760           } else {
01761             //make sure won't be added
01762             int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01763             int bit = j & COINFACTORIZATION_MASK_PER_INT;
01764             
01765             temp2[word] = temp2[word] | ( 1 << bit );   //say already in counts
01766           }
01767         }
01768       }
01769       numberInColumn[iColumn] = put - startColumnThis;
01770       //move largest
01771       if ( positionLargest >= 0 ) {
01772         value = element[positionLargest];
01773         iRow = indexRow[positionLargest];
01774         element[positionLargest] = element[startColumnThis];
01775         indexRow[positionLargest] = indexRow[startColumnThis];
01776         element[startColumnThis] = value;
01777         indexRow[startColumnThis] = iRow;
01778       }
01779       //linked list for column
01780       if ( nextCount[iColumn + numberRows_] != -2 ) {
01781         //modify linked list
01782         deleteLink ( iColumn + numberRows_ );
01783         addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01784       }
01785       temp2 += increment2;
01786     }
01787     //get space for row list
01788     unsigned int *putBase = workArea2;
01789     int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01790     int i = 0;
01791     
01792     // do linked lists and update counts
01793     while ( bigLoops ) {
01794       bigLoops--;
01795       int bit;
01796       for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01797         unsigned int *putThis = putBase;
01798         int iRow = indexL[i];
01799         
01800         //get space
01801         int number = 0;
01802         int jColumn;
01803         
01804         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01805           unsigned int test = *putThis;
01806           
01807           putThis += increment2;
01808           test = 1 - ( ( test >> bit ) & 1 );
01809           number += test;
01810         }
01811         int next = nextRow[iRow];
01812         CoinBigIndex space;
01813         
01814         space = startRow[next] - startRow[iRow];
01815         number += numberInRow[iRow];
01816         if ( space < number ) {
01817           if ( !getRowSpace ( iRow, number ) ) {
01818             goto BAD_PIVOT;
01819           }
01820         }
01821         // now do
01822         putThis = putBase;
01823         next = nextRow[iRow];
01824         number = numberInRow[iRow];
01825         CoinBigIndex end = startRow[iRow] + number;
01826         int saveIndex = indexColumn[startRow[next]];
01827         
01828         //add in
01829         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01830           unsigned int test = *putThis;
01831           
01832           putThis += increment2;
01833           test = 1 - ( ( test >> bit ) & 1 );
01834           indexColumn[end] = saveColumn[jColumn];
01835           end += test;
01836         }
01837         //put back next one in case zapped
01838         indexColumn[startRow[next]] = saveIndex;
01839         markRow[iRow] = FAC_UNSET;
01840         number = end - startRow[iRow];
01841         numberInRow[iRow] = number;
01842         deleteLink ( iRow );
01843         addLink ( iRow, number );
01844       }
01845       putBase++;
01846     }                           /* endwhile */
01847     int bit;
01848     
01849     for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
01850       unsigned int *putThis = putBase;
01851       int iRow = indexL[i];
01852       
01853       //get space
01854       int number = 0;
01855       int jColumn;
01856       
01857       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01858         unsigned int test = *putThis;
01859         
01860         putThis += increment2;
01861         test = 1 - ( ( test >> bit ) & 1 );
01862         number += test;
01863       }
01864       int next = nextRow[iRow];
01865       CoinBigIndex space;
01866       
01867       space = startRow[next] - startRow[iRow];
01868       number += numberInRow[iRow];
01869       if ( space < number ) {
01870         if ( !getRowSpace ( iRow, number ) ) {
01871           goto BAD_PIVOT;
01872         }
01873       }
01874       // now do
01875       putThis = putBase;
01876       next = nextRow[iRow];
01877       number = numberInRow[iRow];
01878       CoinBigIndex end = startRow[iRow] + number;
01879       int saveIndex;
01880       
01881       saveIndex = indexColumn[startRow[next]];
01882       
01883       //add in
01884       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01885         unsigned int test = *putThis;
01886         
01887         putThis += increment2;
01888         test = 1 - ( ( test >> bit ) & 1 );
01889         
01890         indexColumn[end] = saveColumn[jColumn];
01891         end += test;
01892       }
01893       indexColumn[startRow[next]] = saveIndex;
01894       markRow[iRow] = FAC_UNSET;
01895       number = end - startRow[iRow];
01896       numberInRow[iRow] = number;
01897       deleteLink ( iRow );
01898       addLink ( iRow, number );
01899     }
01900     markRow[iPivotRow] = FAC_UNSET;
01901     //modify linked list for pivots
01902     deleteLink ( iPivotRow );
01903     deleteLink ( iPivotColumn + numberRows_ );
01904     totalElements_ += added;
01905     goodPivot= true;
01906     // **** UGLY UGLY UGLY
01907   }
01908  BAD_PIVOT:
01909 
01910   ;
01911 }
01912 #undef FAC_UNSET
01913 #endif

Generated on Wed Apr 7 03:06:52 2010 by  doxygen 1.4.7