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