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

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

Generated on Mon Sep 14 03:06:04 2009 by  doxygen 1.4.7