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