• Main Page
  • Related Pages
  • Modules
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

/build/buildd/coinutils-2.6.4/CoinUtils/src/CoinFactorization.hpp

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

Generated on Fri Oct 15 2010 18:21:02 by  doxygen 1.7.1