CoinFactorization.hpp

Go to the documentation of this file.
00001 /* $Id: CoinFactorization.hpp 1767 2015-01-05 12:36:13Z forrest $ */
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   }
00167   inline int *lastRow (  ) const {
00168     return lastRow_.array();
00169   }
00172   inline int *pivotColumnBack (  ) const {
00173     //return firstCount_.array();
00174     return pivotColumnBack_.array();
00175   }
00177   inline CoinBigIndex * startRowL() const
00178   { return startRowL_.array();}
00179 
00181   inline CoinBigIndex * startColumnL() const
00182   { return startColumnL_.array();}
00183 
00185   inline int * indexColumnL() const
00186   { return indexColumnL_.array();}
00187 
00189   inline int * indexRowL() const
00190   { return indexRowL_.array();}
00191 
00193   inline CoinFactorizationDouble * elementByRowL() const
00194   { return elementByRowL_.array();}
00195 
00197   inline int numberRowsExtra (  ) const {
00198     return numberRowsExtra_;
00199   }
00201   inline void setNumberRows(int value)
00202   { numberRows_ = value; }
00204   inline int numberRows (  ) const {
00205     return numberRows_;
00206   }
00208   inline CoinBigIndex numberL() const
00209   { return numberL_;}
00210 
00212   inline CoinBigIndex baseL() const
00213   { return baseL_;}
00215   inline int maximumRowsExtra (  ) const {
00216     return maximumRowsExtra_;
00217   }
00219   inline int numberColumns (  ) const {
00220     return numberColumns_;
00221   }
00223   inline int numberElements (  ) const {
00224     return totalElements_;
00225   }
00227   inline int numberForrestTomlin (  ) const {
00228     return numberInColumn_.array()[numberColumnsExtra_];
00229   }
00231   inline int numberGoodColumns (  ) const {
00232     return numberGoodU_;
00233   }
00235   inline double areaFactor (  ) const {
00236     return areaFactor_;
00237   }
00238   inline void areaFactor ( double value ) {
00239     areaFactor_=value;
00240   }
00242   double adjustedAreaFactor() const;
00244   inline void relaxAccuracyCheck(double value)
00245   { relaxCheck_ = value;}
00246   inline double getAccuracyCheck() const
00247   { return relaxCheck_;}
00249   inline int messageLevel (  ) const {
00250     return messageLevel_ ;
00251   }
00252   void messageLevel (  int value );
00254   inline int maximumPivots (  ) const {
00255     return maximumPivots_ ;
00256   }
00257   void maximumPivots (  int value );
00258 
00260   inline int denseThreshold() const 
00261     { return denseThreshold_;}
00263   inline void setDenseThreshold(int value)
00264     { denseThreshold_ = value;}
00266   inline double pivotTolerance (  ) const {
00267     return pivotTolerance_ ;
00268   }
00269   void pivotTolerance (  double value );
00271   inline double zeroTolerance (  ) const {
00272     return zeroTolerance_ ;
00273   }
00274   void zeroTolerance (  double value );
00275 #ifndef COIN_FAST_CODE
00277   inline double slackValue (  ) const {
00278     return slackValue_ ;
00279   }
00280   void slackValue (  double value );
00281 #endif
00283   double maximumCoefficient() const;
00285   inline bool forrestTomlin() const
00286   { return doForrestTomlin_;}
00287   inline void setForrestTomlin(bool value)
00288   { doForrestTomlin_=value;}
00290   inline bool spaceForForrestTomlin() const
00291   {
00292     CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
00293     CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
00294     return (space>=0)&&doForrestTomlin_;
00295   }
00297 
00300 
00302   inline int numberDense() const
00303   { return numberDense_;}
00304 
00306   inline CoinBigIndex numberElementsU (  ) const {
00307     return lengthU_;
00308   }
00310   inline void setNumberElementsU(CoinBigIndex value)
00311   { lengthU_ = value; }
00313   inline CoinBigIndex lengthAreaU (  ) const {
00314     return lengthAreaU_;
00315   }
00317   inline CoinBigIndex numberElementsL (  ) const {
00318     return lengthL_;
00319   }
00321   inline CoinBigIndex lengthAreaL (  ) const {
00322     return lengthAreaL_;
00323   }
00325   inline CoinBigIndex numberElementsR (  ) const {
00326     return lengthR_;
00327   }
00329   inline CoinBigIndex numberCompressions() const
00330   { return numberCompressions_;}
00332   inline int * numberInRow() const
00333   { return numberInRow_.array();}
00335   inline int * numberInColumn() const
00336   { return numberInColumn_.array();}
00338   inline CoinFactorizationDouble * elementU() const
00339   { return elementU_.array();}
00341   inline int * indexRowU() const
00342   { return indexRowU_.array();}
00344   inline CoinBigIndex * startColumnU() const
00345   { return startColumnU_.array();}
00347   inline int maximumColumnsExtra()
00348   { return maximumColumnsExtra_;}
00352   inline int biasLU() const
00353   { return biasLU_;}
00354   inline void setBiasLU(int value)
00355   { biasLU_=value;}
00361   inline int persistenceFlag() const
00362   { return persistenceFlag_;}
00363   void setPersistenceFlag(int value);
00365 
00368 
00376   int replaceColumn ( CoinIndexedVector * regionSparse,
00377                       int pivotRow,
00378                       double pivotCheck ,
00379                       bool checkBeforeModifying=false,
00380                       double acceptablePivot=1.0e-8);
00385   void replaceColumnU ( CoinIndexedVector * regionSparse,
00386                         CoinBigIndex * deleted,
00387                         int internalPivotRow);
00388 #ifdef ABC_USE_COIN_FACTORIZATION
00389 
00391   CoinIndexedVector * fakeVector(CoinIndexedVector * vector,
00392                                  int already=0) const;
00393   void deleteFakeVector(CoinIndexedVector * vector,
00394                         CoinIndexedVector * fakeVector) const;
00399   double checkReplacePart1 (  CoinIndexedVector * regionSparse,
00400                                      int pivotRow);
00405   double checkReplacePart1 (  CoinIndexedVector * regionSparse,
00406                                       CoinIndexedVector * partialUpdate,
00407                                      int pivotRow);
00410   int checkReplacePart2 ( int pivotRow,
00411                                   double btranAlpha, 
00412                                   double ftranAlpha, 
00413                                   double ftAlpha,
00414                                   double acceptablePivot = 1.0e-8);
00417   void replaceColumnPart3 ( CoinIndexedVector * regionSparse,
00418                             int pivotRow,
00419                             double alpha );
00422   void replaceColumnPart3 ( CoinIndexedVector * regionSparse,
00423                             CoinIndexedVector * partialUpdate,
00424                             int pivotRow,
00425                             double alpha );
00433   int updateColumnFT ( CoinIndexedVector & regionSparse);
00434   int updateColumnFTPart1 ( CoinIndexedVector & regionSparse) ;
00435   void updateColumnFTPart2 ( CoinIndexedVector & regionSparse) ;
00439   void updateColumnFT ( CoinIndexedVector & regionSparseFT,
00440                         CoinIndexedVector & partialUpdate,
00441                         int which);
00443   int updateColumn ( CoinIndexedVector & regionSparse) const;
00448   int updateTwoColumnsFT ( CoinIndexedVector & regionSparseFT,
00449                            CoinIndexedVector & regionSparseOther);
00451   int updateColumnTranspose ( CoinIndexedVector & regionSparse) const;
00453   void updateColumnCpu ( CoinIndexedVector & regionSparse,int whichCpu) const;
00455   void updateColumnTransposeCpu ( CoinIndexedVector & regionSparse,int whichCpu) const;
00457   void updateFullColumn ( CoinIndexedVector & regionSparse) const;
00459   void updateFullColumnTranspose ( CoinIndexedVector & regionSparse) const;
00461   void updateWeights ( CoinIndexedVector & regionSparse) const;
00463   inline bool wantsTableauColumn() const
00464   {return false;}
00466   inline double minimumPivotTolerance (  ) const {
00467     return pivotTolerance_ ;
00468   }
00469   inline void minimumPivotTolerance (  double value )
00470   { pivotTolerance(value);}
00472   inline void setParallelMode(int value)
00473   { parallelMode_=value;}
00475   inline void setSolveMode(int value)
00476   { parallelMode_ &= 3;parallelMode_ |= (value<<2);}
00478   inline int solveMode() const
00479   { return parallelMode_ >> 2;}
00481   void updatePartialUpdate(CoinIndexedVector & partialUpdate);
00483   void makeNonSingular(int *  COIN_RESTRICT sequence);
00484 #endif
00485 
00486 
00496   int updateColumnFT ( CoinIndexedVector * regionSparse,
00497                        CoinIndexedVector * regionSparse2);
00500   int updateColumn ( CoinIndexedVector * regionSparse,
00501                      CoinIndexedVector * regionSparse2,
00502                      bool noPermute=false) const;
00508   int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
00509                            CoinIndexedVector * regionSparse2,
00510                            CoinIndexedVector * regionSparse3,
00511                            bool noPermuteRegion3=false) ;
00516   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00517                               CoinIndexedVector * regionSparse2) const;
00519   void goSparse();
00521   inline int sparseThreshold ( ) const
00522   { return sparseThreshold_;}
00524   void sparseThreshold ( int value );
00526 
00527 
00531 
00532   inline void clearArrays()
00533   { gutsOfDestructor();}
00535 
00538 
00541   int add ( CoinBigIndex numberElements,
00542                int indicesRow[],
00543                int indicesColumn[], double elements[] );
00544 
00547   int addColumn ( CoinBigIndex numberElements,
00548                      int indicesRow[], double elements[] );
00549 
00552   int addRow ( CoinBigIndex numberElements,
00553                   int indicesColumn[], double elements[] );
00554 
00556   int deleteColumn ( int Row );
00558   int deleteRow ( int Row );
00559 
00563   int replaceRow ( int whichRow, int numberElements,
00564                       const int indicesColumn[], const double elements[] );
00566   void emptyRows(int numberToEmpty, const int which[]);
00568 
00569 
00570   void checkSparse();
00572 #if 0 //def CLP_FACTORIZATION_INSTRUMENT
00573   inline bool collectStatistics() const
00574   { return collectStatistics_;}
00576   inline void setCollectStatistics(bool onOff) const
00577   { collectStatistics_ = onOff;}
00578 #else
00579   inline bool collectStatistics() const
00580   { return true;}
00582   inline void setCollectStatistics(bool onOff) const
00583   { }
00584 #endif
00586   void gutsOfDestructor(int type=1);
00588   void gutsOfInitialize(int type);
00589   void gutsOfCopy(const CoinFactorization &other);
00590 
00592   void resetStatistics();
00593 
00594 
00596 
00598 
00599   void getAreas ( int numberRows,
00600                   int numberColumns,
00601                   CoinBigIndex maximumL,
00602                   CoinBigIndex maximumU );
00603 
00606   void preProcess ( int state,
00607                     int possibleDuplicates = -1 );
00609   int factor (  );
00610 protected:
00613   int factorSparse (  );
00616   int factorSparseSmall (  );
00619   int factorSparseLarge (  );
00622   int factorDense (  );
00623 
00625   bool pivotOneOtherRow ( int pivotRow,
00626                           int pivotColumn );
00628   bool pivotRowSingleton ( int pivotRow,
00629                            int pivotColumn );
00631   bool pivotColumnSingleton ( int pivotRow,
00632                               int pivotColumn );
00633 
00638   bool getColumnSpace ( int iColumn,
00639                         int extraNeeded );
00640 
00643   bool reorderU();
00647   bool getColumnSpaceIterateR ( int iColumn, double value,
00648                                int iRow);
00654   CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00655                                int iRow);
00659   bool getRowSpace ( int iRow, int extraNeeded );
00660 
00664   bool getRowSpaceIterate ( int iRow,
00665                             int extraNeeded );
00667   void checkConsistency (  );
00669   inline void addLink ( int index, int count ) {
00670     int *nextCount = nextCount_.array();
00671     int *firstCount = firstCount_.array();
00672     int *lastCount = lastCount_.array();
00673     int next = firstCount[count];
00674       lastCount[index] = -2 - count;
00675     if ( next < 0 ) {
00676       //first with that count
00677       firstCount[count] = index;
00678       nextCount[index] = -1;
00679     } else {
00680       firstCount[count] = index;
00681       nextCount[index] = next;
00682       lastCount[next] = index;
00683   }}
00685   inline void deleteLink ( int index ) {
00686     int *nextCount = nextCount_.array();
00687     int *firstCount = firstCount_.array();
00688     int *lastCount = lastCount_.array();
00689     int next = nextCount[index];
00690     int last = lastCount[index];
00691     if ( last >= 0 ) {
00692       nextCount[last] = next;
00693     } else {
00694       int count = -last - 2;
00695 
00696       firstCount[count] = next;
00697     }
00698     if ( next >= 0 ) {
00699       lastCount[next] = last;
00700     }
00701     nextCount[index] = -2;
00702     lastCount[index] = -2;
00703     return;
00704   }
00706   void separateLinks(int count,bool rowsFirst);
00708   void cleanup (  );
00709 
00711   void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00713   void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00715   void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00717   void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00718 
00720   void updateColumnR ( CoinIndexedVector * region ) const;
00723   void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00724 
00726   void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00727 
00729   void updateColumnUSparse ( CoinIndexedVector * regionSparse, 
00730                              int * indexIn) const;
00732   void updateColumnUSparsish ( CoinIndexedVector * regionSparse, 
00733                                int * indexIn) const;
00735   int updateColumnUDensish ( double * COIN_RESTRICT region, 
00736                              int * COIN_RESTRICT regionIndex) const;
00738   void updateTwoColumnsUDensish (
00739                                  int & numberNonZero1,
00740                                  double * COIN_RESTRICT region1, 
00741                                  int * COIN_RESTRICT index1,
00742                                  int & numberNonZero2,
00743                                  double * COIN_RESTRICT region2, 
00744                                  int * COIN_RESTRICT index2) const;
00746   void updateColumnPFI ( CoinIndexedVector * regionSparse) const; 
00748   void permuteBack ( CoinIndexedVector * regionSparse, 
00749                      CoinIndexedVector * outVector) const;
00750 
00752   void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00755   void updateColumnTransposeU ( CoinIndexedVector * region,
00756                                 int smallestIndex) const;
00759   void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00760                                         int smallestIndex) const;
00763   void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00764                                        int smallestIndex) const;
00767   void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00770   void updateColumnTransposeUByColumn ( CoinIndexedVector * region,
00771                                         int smallestIndex) const;
00772 
00774   void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00776   void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00778   void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00779 
00781   void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00783   void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00785   void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00787   void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00789   void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00790 public:
00795   int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00796                          int pivotRow, double alpha);
00797 protected:
00800   int checkPivot(double saveFromU, double oldPivot) const;
00801   /********************************* START LARGE TEMPLATE ********/
00802 #ifdef INT_IS_8
00803 #define COINFACTORIZATION_BITS_PER_INT 64
00804 #define COINFACTORIZATION_SHIFT_PER_INT 6
00805 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00806 #else
00807 #define COINFACTORIZATION_BITS_PER_INT 32
00808 #define COINFACTORIZATION_SHIFT_PER_INT 5
00809 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00810 #endif
00811   template <class T>  inline bool
00812   pivot ( int pivotRow,
00813           int pivotColumn,
00814           CoinBigIndex pivotRowPosition,
00815           CoinBigIndex pivotColumnPosition,
00816           CoinFactorizationDouble work[],
00817           unsigned int workArea2[],
00818           int increment2,
00819           T markRow[] ,
00820           int largeInteger)
00821 {
00822   int *indexColumnU = indexColumnU_.array();
00823   CoinBigIndex *startColumnU = startColumnU_.array();
00824   int *numberInColumn = numberInColumn_.array();
00825   CoinFactorizationDouble *elementU = elementU_.array();
00826   int *indexRowU = indexRowU_.array();
00827   CoinBigIndex *startRowU = startRowU_.array();
00828   int *numberInRow = numberInRow_.array();
00829   CoinFactorizationDouble *elementL = elementL_.array();
00830   int *indexRowL = indexRowL_.array();
00831   int *saveColumn = saveColumn_.array();
00832   int *nextRow = nextRow_.array();
00833   int *lastRow = lastRow_.array() ;
00834 
00835   //store pivot columns (so can easily compress)
00836   int numberInPivotRow = numberInRow[pivotRow] - 1;
00837   CoinBigIndex startColumn = startColumnU[pivotColumn];
00838   int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00839   CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00840   int put = 0;
00841   CoinBigIndex startRow = startRowU[pivotRow];
00842   CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00843 
00844   if ( pivotColumnPosition < 0 ) {
00845     for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00846       int iColumn = indexColumnU[pivotColumnPosition];
00847       if ( iColumn != pivotColumn ) {
00848         saveColumn[put++] = iColumn;
00849       } else {
00850         break;
00851       }
00852     }
00853   } else {
00854     for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00855       saveColumn[put++] = indexColumnU[i];
00856     }
00857   }
00858   assert (pivotColumnPosition<endRow);
00859   assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00860   pivotColumnPosition++;
00861   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00862     saveColumn[put++] = indexColumnU[pivotColumnPosition];
00863   }
00864   //take out this bit of indexColumnU
00865   int next = nextRow[pivotRow];
00866   int last = lastRow[pivotRow];
00867 
00868   nextRow[last] = next;
00869   lastRow[next] = last;
00870   nextRow[pivotRow] = numberGoodU_;     //use for permute
00871   lastRow[pivotRow] = -2;
00872   numberInRow[pivotRow] = 0;
00873   //store column in L, compress in U and take column out
00874   CoinBigIndex l = lengthL_;
00875 
00876   if ( l + numberInPivotColumn > lengthAreaL_ ) {
00877     //need more memory
00878     if ((messageLevel_&4)!=0) 
00879       printf("more memory needed in middle of invert\n");
00880     return false;
00881   }
00882   //l+=currentAreaL_->elementByColumn-elementL;
00883   CoinBigIndex lSave = l;
00884 
00885   CoinBigIndex * startColumnL = startColumnL_.array();
00886   startColumnL[numberGoodL_] = l;       //for luck and first time
00887   numberGoodL_++;
00888   startColumnL[numberGoodL_] = l + numberInPivotColumn;
00889   lengthL_ += numberInPivotColumn;
00890   if ( pivotRowPosition < 0 ) {
00891     for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00892       int iRow = indexRowU[pivotRowPosition];
00893       if ( iRow != pivotRow ) {
00894         indexRowL[l] = iRow;
00895         elementL[l] = elementU[pivotRowPosition];
00896         markRow[iRow] = static_cast<T>(l - lSave);
00897         l++;
00898         //take out of row list
00899         CoinBigIndex start = startRowU[iRow];
00900         CoinBigIndex end = start + numberInRow[iRow];
00901         CoinBigIndex where = start;
00902 
00903         while ( indexColumnU[where] != pivotColumn ) {
00904           where++;
00905         }                       /* endwhile */
00906 #if DEBUG_COIN
00907         if ( where >= end ) {
00908           abort (  );
00909         }
00910 #endif
00911         indexColumnU[where] = indexColumnU[end - 1];
00912         numberInRow[iRow]--;
00913       } else {
00914         break;
00915       }
00916     }
00917   } else {
00918     CoinBigIndex i;
00919 
00920     for ( i = startColumn; i < pivotRowPosition; i++ ) {
00921       int iRow = indexRowU[i];
00922 
00923       markRow[iRow] = static_cast<T>(l - lSave);
00924       indexRowL[l] = iRow;
00925       elementL[l] = elementU[i];
00926       l++;
00927       //take out of row list
00928       CoinBigIndex start = startRowU[iRow];
00929       CoinBigIndex end = start + numberInRow[iRow];
00930       CoinBigIndex where = start;
00931 
00932       while ( indexColumnU[where] != pivotColumn ) {
00933         where++;
00934       }                         /* endwhile */
00935 #if DEBUG_COIN
00936       if ( where >= end ) {
00937         abort (  );
00938       }
00939 #endif
00940       indexColumnU[where] = indexColumnU[end - 1];
00941       numberInRow[iRow]--;
00942       assert (numberInRow[iRow]>=0);
00943     }
00944   }
00945   assert (pivotRowPosition<endColumn);
00946   assert (indexRowU[pivotRowPosition]==pivotRow);
00947   CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
00948   CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
00949 
00950   pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00951   pivotRowPosition++;
00952   for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00953     int iRow = indexRowU[pivotRowPosition];
00954     
00955     markRow[iRow] = static_cast<T>(l - lSave);
00956     indexRowL[l] = iRow;
00957     elementL[l] = elementU[pivotRowPosition];
00958     l++;
00959     //take out of row list
00960     CoinBigIndex start = startRowU[iRow];
00961     CoinBigIndex end = start + numberInRow[iRow];
00962     CoinBigIndex where = start;
00963     
00964     while ( indexColumnU[where] != pivotColumn ) {
00965       where++;
00966     }                           /* endwhile */
00967 #if DEBUG_COIN
00968     if ( where >= end ) {
00969       abort (  );
00970     }
00971 #endif
00972     indexColumnU[where] = indexColumnU[end - 1];
00973     numberInRow[iRow]--;
00974     assert (numberInRow[iRow]>=0);
00975   }
00976   markRow[pivotRow] = static_cast<T>(largeInteger);
00977   //compress pivot column (move pivot to front including saved)
00978   numberInColumn[pivotColumn] = 0;
00979   //use end of L for temporary space
00980   int *indexL = &indexRowL[lSave];
00981   CoinFactorizationDouble *multipliersL = &elementL[lSave];
00982 
00983   //adjust
00984   int j;
00985 
00986   for ( j = 0; j < numberInPivotColumn; j++ ) {
00987     multipliersL[j] *= pivotMultiplier;
00988   }
00989   //zero out fill
00990   CoinBigIndex iErase;
00991   for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00992         iErase++ ) {
00993     workArea2[iErase] = 0;
00994   }
00995   CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00996   unsigned int *temp2 = workArea2;
00997   int * nextColumn = nextColumn_.array();
00998 
00999   //pack down and move to work
01000   int jColumn;
01001   for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01002     int iColumn = saveColumn[jColumn];
01003     CoinBigIndex startColumn = startColumnU[iColumn];
01004     CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
01005     int iRow = indexRowU[startColumn];
01006     CoinFactorizationDouble value = elementU[startColumn];
01007     double largest;
01008     CoinBigIndex put = startColumn;
01009     CoinBigIndex positionLargest = -1;
01010     CoinFactorizationDouble thisPivotValue = 0.0;
01011 
01012     //compress column and find largest not updated
01013     bool checkLargest;
01014     int mark = markRow[iRow];
01015 
01016     if ( mark == largeInteger+1 ) {
01017       largest = fabs ( value );
01018       positionLargest = put;
01019       put++;
01020       checkLargest = false;
01021     } else {
01022       //need to find largest
01023       largest = 0.0;
01024       checkLargest = true;
01025       if ( mark != largeInteger ) {
01026         //will be updated
01027         work[mark] = value;
01028         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01029         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01030 
01031         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
01032         added--;
01033       } else {
01034         thisPivotValue = value;
01035       }
01036     }
01037     CoinBigIndex i;
01038     for ( i = startColumn + 1; i < endColumn; i++ ) {
01039       iRow = indexRowU[i];
01040       value = elementU[i];
01041       int mark = markRow[iRow];
01042 
01043       if ( mark == largeInteger+1 ) {
01044         //keep
01045         indexRowU[put] = iRow;
01046         elementU[put] = value;
01047         if ( checkLargest ) {
01048           double absValue = fabs ( value );
01049 
01050           if ( absValue > largest ) {
01051             largest = absValue;
01052             positionLargest = put;
01053           }
01054         }
01055         put++;
01056       } else if ( mark != largeInteger ) {
01057         //will be updated
01058         work[mark] = value;
01059         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01060         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01061 
01062         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
01063         added--;
01064       } else {
01065         thisPivotValue = value;
01066       }
01067     }
01068     //slot in pivot
01069     elementU[put] = elementU[startColumn];
01070     indexRowU[put] = indexRowU[startColumn];
01071     if ( positionLargest == startColumn ) {
01072       positionLargest = put;    //follow if was largest
01073     }
01074     put++;
01075     elementU[startColumn] = thisPivotValue;
01076     indexRowU[startColumn] = pivotRow;
01077     //clean up counts
01078     startColumn++;
01079     numberInColumn[iColumn] = put - startColumn;
01080     int * numberInColumnPlus = numberInColumnPlus_.array();
01081     numberInColumnPlus[iColumn]++;
01082     startColumnU[iColumn]++;
01083     //how much space have we got
01084     int next = nextColumn[iColumn];
01085     CoinBigIndex space;
01086 
01087     space = startColumnU[next] - put - numberInColumnPlus[next];
01088     //assume no zero elements
01089     if ( numberInPivotColumn > space ) {
01090       //getColumnSpace also moves fixed part
01091       if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
01092         return false;
01093       }
01094       //redo starts
01095       if (positionLargest >= 0)
01096          positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
01097       startColumn = startColumnU[iColumn];
01098       put = startColumn + numberInColumn[iColumn];
01099     }
01100     double tolerance = zeroTolerance_;
01101 
01102     int *nextCount = nextCount_.array();
01103     for ( j = 0; j < numberInPivotColumn; j++ ) {
01104       value = work[j] - thisPivotValue * multipliersL[j];
01105       double absValue = fabs ( value );
01106 
01107       if ( absValue > tolerance ) {
01108         work[j] = 0.0;
01109         assert (put<lengthAreaU_); 
01110         elementU[put] = value;
01111         indexRowU[put] = indexL[j];
01112         if ( absValue > largest ) {
01113           largest = absValue;
01114           positionLargest = put;
01115         }
01116         put++;
01117       } else {
01118         work[j] = 0.0;
01119         added--;
01120         int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01121         int bit = j & COINFACTORIZATION_MASK_PER_INT;
01122 
01123         if ( temp2[word] & ( 1 << bit ) ) {
01124           //take out of row list
01125           iRow = indexL[j];
01126           CoinBigIndex start = startRowU[iRow];
01127           CoinBigIndex end = start + numberInRow[iRow];
01128           CoinBigIndex where = start;
01129 
01130           while ( indexColumnU[where] != iColumn ) {
01131             where++;
01132           }                     /* endwhile */
01133 #if DEBUG_COIN
01134           if ( where >= end ) {
01135             abort (  );
01136           }
01137 #endif
01138           indexColumnU[where] = indexColumnU[end - 1];
01139           numberInRow[iRow]--;
01140         } else {
01141           //make sure won't be added
01142           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01143           int bit = j & COINFACTORIZATION_MASK_PER_INT;
01144 
01145           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
01146         }
01147       }
01148     }
01149     numberInColumn[iColumn] = put - startColumn;
01150     //move largest
01151     if ( positionLargest >= 0 ) {
01152       value = elementU[positionLargest];
01153       iRow = indexRowU[positionLargest];
01154       elementU[positionLargest] = elementU[startColumn];
01155       indexRowU[positionLargest] = indexRowU[startColumn];
01156       elementU[startColumn] = value;
01157       indexRowU[startColumn] = iRow;
01158     }
01159     //linked list for column
01160     if ( nextCount[iColumn + numberRows_] != -2 ) {
01161       //modify linked list
01162       deleteLink ( iColumn + numberRows_ );
01163       addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01164     }
01165     temp2 += increment2;
01166   }
01167   //get space for row list
01168   unsigned int *putBase = workArea2;
01169   int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01170   int i = 0;
01171 
01172   // do linked lists and update counts
01173   while ( bigLoops ) {
01174     bigLoops--;
01175     int bit;
01176     for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01177       unsigned int *putThis = putBase;
01178       int iRow = indexL[i];
01179 
01180       //get space
01181       int number = 0;
01182       int jColumn;
01183 
01184       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01185         unsigned int test = *putThis;
01186 
01187         putThis += increment2;
01188         test = 1 - ( ( test >> bit ) & 1 );
01189         number += test;
01190       }
01191       int next = nextRow[iRow];
01192       CoinBigIndex space;
01193 
01194       space = startRowU[next] - startRowU[iRow];
01195       number += numberInRow[iRow];
01196       if ( space < number ) {
01197         if ( !getRowSpace ( iRow, number ) ) {
01198           return false;
01199         }
01200       }
01201       // now do
01202       putThis = putBase;
01203       next = nextRow[iRow];
01204       number = numberInRow[iRow];
01205       CoinBigIndex end = startRowU[iRow] + number;
01206       int saveIndex = indexColumnU[startRowU[next]];
01207 
01208       //add in
01209       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01210         unsigned int test = *putThis;
01211 
01212         putThis += increment2;
01213         test = 1 - ( ( test >> bit ) & 1 );
01214         indexColumnU[end] = saveColumn[jColumn];
01215         end += test;
01216       }
01217       //put back next one in case zapped
01218       indexColumnU[startRowU[next]] = saveIndex;
01219       markRow[iRow] = static_cast<T>(largeInteger+1);
01220       number = end - startRowU[iRow];
01221       numberInRow[iRow] = number;
01222       deleteLink ( iRow );
01223       addLink ( iRow, number );
01224     }
01225     putBase++;
01226   }                             /* endwhile */
01227   int bit;
01228 
01229   for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01230     unsigned int *putThis = putBase;
01231     int iRow = indexL[i];
01232 
01233     //get space
01234     int number = 0;
01235     int jColumn;
01236 
01237     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01238       unsigned int test = *putThis;
01239 
01240       putThis += increment2;
01241       test = 1 - ( ( test >> bit ) & 1 );
01242       number += test;
01243     }
01244     int next = nextRow[iRow];
01245     CoinBigIndex space;
01246 
01247     space = startRowU[next] - startRowU[iRow];
01248     number += numberInRow[iRow];
01249     if ( space < number ) {
01250       if ( !getRowSpace ( iRow, number ) ) {
01251         return false;
01252       }
01253     }
01254     // now do
01255     putThis = putBase;
01256     next = nextRow[iRow];
01257     number = numberInRow[iRow];
01258     CoinBigIndex end = startRowU[iRow] + number;
01259     int saveIndex;
01260 
01261     saveIndex = indexColumnU[startRowU[next]];
01262 
01263     //add in
01264     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01265       unsigned int test = *putThis;
01266 
01267       putThis += increment2;
01268       test = 1 - ( ( test >> bit ) & 1 );
01269 
01270       indexColumnU[end] = saveColumn[jColumn];
01271       end += test;
01272     }
01273     indexColumnU[startRowU[next]] = saveIndex;
01274     markRow[iRow] = static_cast<T>(largeInteger+1);
01275     number = end - startRowU[iRow];
01276     numberInRow[iRow] = number;
01277     deleteLink ( iRow );
01278     addLink ( iRow, number );
01279   }
01280   markRow[pivotRow] = static_cast<T>(largeInteger+1);
01281   //modify linked list for pivots
01282   deleteLink ( pivotRow );
01283   deleteLink ( pivotColumn + numberRows_ );
01284   totalElements_ += added;
01285   return true;
01286 }
01287 
01288   /********************************* END LARGE TEMPLATE ********/
01290 
01291 protected:
01292 
01295 
01296   double pivotTolerance_;
01298   double zeroTolerance_;
01299 #ifndef COIN_FAST_CODE
01301   double slackValue_;
01302 #else
01303 #ifndef slackValue_
01304 #define slackValue_ -1.0
01305 #endif
01306 #endif
01308   double areaFactor_;
01310   double relaxCheck_;
01312   int numberRows_;
01314   int numberRowsExtra_;
01316   int maximumRowsExtra_;
01318   int numberColumns_;
01320   int numberColumnsExtra_;
01322   int maximumColumnsExtra_;
01324   int numberGoodU_;
01326   int numberGoodL_;
01328   int maximumPivots_;
01330   int numberPivots_;
01333   CoinBigIndex totalElements_;
01335   CoinBigIndex factorElements_;
01337   CoinIntArrayWithLength pivotColumn_;
01339   CoinIntArrayWithLength permute_;
01341   CoinIntArrayWithLength permuteBack_;
01343   CoinIntArrayWithLength pivotColumnBack_;
01345   int status_;
01346 
01351   //int increasingRows_;
01352 
01354   int numberTrials_;
01356   CoinBigIndexArrayWithLength startRowU_;
01357 
01359   CoinIntArrayWithLength numberInRow_;
01360 
01362   CoinIntArrayWithLength numberInColumn_;
01363 
01365   CoinIntArrayWithLength numberInColumnPlus_;
01366 
01369   CoinIntArrayWithLength firstCount_;
01370 
01372   CoinIntArrayWithLength nextCount_;
01373 
01375   CoinIntArrayWithLength lastCount_;
01376 
01378   CoinIntArrayWithLength nextColumn_;
01379 
01381   CoinIntArrayWithLength lastColumn_;
01382 
01384   CoinIntArrayWithLength nextRow_;
01385 
01387   CoinIntArrayWithLength lastRow_;
01388 
01390   CoinIntArrayWithLength saveColumn_;
01391 
01393   CoinIntArrayWithLength markRow_;
01394 
01396   int messageLevel_;
01397 
01399   int biggerDimension_;
01400 
01402   CoinIntArrayWithLength indexColumnU_;
01403 
01405   CoinIntArrayWithLength pivotRowL_;
01406 
01408   CoinFactorizationDoubleArrayWithLength pivotRegion_;
01409 
01411   int numberSlacks_;
01412 
01414   int numberU_;
01415 
01417   CoinBigIndex maximumU_;
01418 
01420   //int baseU_;
01421 
01423   CoinBigIndex lengthU_;
01424 
01426   CoinBigIndex lengthAreaU_;
01427 
01429   CoinFactorizationDoubleArrayWithLength elementU_;
01430 
01432   CoinIntArrayWithLength indexRowU_;
01433 
01435   CoinBigIndexArrayWithLength startColumnU_;
01436 
01438   CoinBigIndexArrayWithLength convertRowToColumnU_;
01439 
01441   CoinBigIndex numberL_;
01442 
01444   CoinBigIndex baseL_;
01445 
01447   CoinBigIndex lengthL_;
01448 
01450   CoinBigIndex lengthAreaL_;
01451 
01453   CoinFactorizationDoubleArrayWithLength elementL_;
01454 
01456   CoinIntArrayWithLength indexRowL_;
01457 
01459   CoinBigIndexArrayWithLength startColumnL_;
01460 
01462   bool doForrestTomlin_;
01463 
01465   int numberR_;
01466 
01468   CoinBigIndex lengthR_;
01469 
01471   CoinBigIndex lengthAreaR_;
01472 
01474   CoinFactorizationDouble *elementR_;
01475 
01477   int *indexRowR_;
01478 
01480   CoinBigIndexArrayWithLength startColumnR_;
01481 
01483   double  * denseArea_;
01484 
01486   double  * denseAreaAddress_;
01487 
01489   int * densePermute_;
01490 
01492   int numberDense_;
01493 
01495   int denseThreshold_;
01496 
01498   CoinFactorizationDoubleArrayWithLength workArea_;
01499 
01501   CoinUnsignedIntArrayWithLength workArea2_;
01502 
01504   CoinBigIndex numberCompressions_;
01505 
01506 public:
01508   mutable double ftranCountInput_;
01509   mutable double ftranCountAfterL_;
01510   mutable double ftranCountAfterR_;
01511   mutable double ftranCountAfterU_;
01512   mutable double btranCountInput_;
01513   mutable double btranCountAfterU_;
01514   mutable double btranCountAfterR_;
01515   mutable double btranCountAfterL_;
01516 
01518   mutable int numberFtranCounts_;
01519   mutable int numberBtranCounts_;
01520 
01522   double ftranAverageAfterL_;
01523   double ftranAverageAfterR_;
01524   double ftranAverageAfterU_;
01525   double btranAverageAfterU_;
01526   double btranAverageAfterR_;
01527   double btranAverageAfterL_;
01528 protected:
01529 
01531 #if 0
01532   mutable bool collectStatistics_;
01533 #else
01534 #define collectStatistics_ 1
01535 #endif
01536 
01538   int sparseThreshold_;
01539 
01541   int sparseThreshold2_;
01542 
01544   CoinBigIndexArrayWithLength startRowL_;
01545 
01547   CoinIntArrayWithLength indexColumnL_;
01548 
01550   CoinFactorizationDoubleArrayWithLength elementByRowL_;
01551 
01553   mutable CoinIntArrayWithLength sparse_;
01557   int biasLU_;
01563   int persistenceFlag_;
01564 #ifdef ABC_USE_COIN_FACTORIZATION
01566   int parallelMode_;
01567 #endif
01568 
01569 };
01570 // Dense coding
01571 #ifdef COIN_HAS_LAPACK
01572 #ifndef COIN_FACTORIZATION_DENSE_CODE
01573 #define COIN_FACTORIZATION_DENSE_CODE 1
01574 #endif
01575 #endif
01576 #ifdef COIN_FACTORIZATION_DENSE_CODE
01577 /* Type of Fortran integer translated into C */
01578 #ifndef ipfint
01579 //typedef ipfint FORTRAN_INTEGER_TYPE ;
01580 typedef int ipfint;
01581 typedef const int cipfint;
01582 #endif
01583 #endif
01584 #endif
01585 // Extra for ugly include
01586 #ifdef UGLY_COIN_FACTOR_CODING
01587 #define FAC_UNSET (FAC_SET+1)
01588 {
01589   goodPivot=false;
01590   //store pivot columns (so can easily compress)
01591   CoinBigIndex startColumnThis = startColumn[iPivotColumn];
01592   CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
01593   int put = 0;
01594   CoinBigIndex startRowThis = startRow[iPivotRow];
01595   CoinBigIndex endRow = startRowThis + numberDoRow + 1;
01596   if ( pivotColumnPosition < 0 ) {
01597     for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01598       int iColumn = indexColumn[pivotColumnPosition];
01599       if ( iColumn != iPivotColumn ) {
01600         saveColumn[put++] = iColumn;
01601       } else {
01602         break;
01603       }
01604     }
01605   } else {
01606     for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
01607       saveColumn[put++] = indexColumn[i];
01608     }
01609   }
01610   assert (pivotColumnPosition<endRow);
01611   assert (indexColumn[pivotColumnPosition]==iPivotColumn);
01612   pivotColumnPosition++;
01613   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01614     saveColumn[put++] = indexColumn[pivotColumnPosition];
01615   }
01616   //take out this bit of indexColumn
01617   int next = nextRow[iPivotRow];
01618   int last = lastRow[iPivotRow];
01619   
01620   nextRow[last] = next;
01621   lastRow[next] = last;
01622   nextRow[iPivotRow] = numberGoodU_;    //use for permute
01623   lastRow[iPivotRow] = -2;
01624   numberInRow[iPivotRow] = 0;
01625   //store column in L, compress in U and take column out
01626   CoinBigIndex l = lengthL_;
01627   // **** HORRID coding coming up but a goto seems best!
01628   {
01629     if ( l + numberDoColumn > lengthAreaL_ ) {
01630       //need more memory
01631       if ((messageLevel_&4)!=0) 
01632         printf("more memory needed in middle of invert\n");
01633       goto BAD_PIVOT;
01634     }
01635     //l+=currentAreaL_->elementByColumn-elementL;
01636     CoinBigIndex lSave = l;
01637     
01638     CoinBigIndex * startColumnL = startColumnL_.array();
01639     startColumnL[numberGoodL_] = l;     //for luck and first time
01640     numberGoodL_++;
01641     startColumnL[numberGoodL_] = l + numberDoColumn;
01642     lengthL_ += numberDoColumn;
01643     if ( pivotRowPosition < 0 ) {
01644       for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01645         int iRow = indexRow[pivotRowPosition];
01646         if ( iRow != iPivotRow ) {
01647           indexRowL[l] = iRow;
01648           elementL[l] = element[pivotRowPosition];
01649           markRow[iRow] = l - lSave;
01650           l++;
01651           //take out of row list
01652           CoinBigIndex start = startRow[iRow];
01653           CoinBigIndex end = start + numberInRow[iRow];
01654           CoinBigIndex where = start;
01655           
01656           while ( indexColumn[where] != iPivotColumn ) {
01657             where++;
01658           }                     /* endwhile */
01659 #if DEBUG_COIN
01660           if ( where >= end ) {
01661             abort (  );
01662           }
01663 #endif
01664           indexColumn[where] = indexColumn[end - 1];
01665           numberInRow[iRow]--;
01666         } else {
01667           break;
01668         }
01669       }
01670     } else {
01671       CoinBigIndex i;
01672       
01673       for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
01674         int iRow = indexRow[i];
01675         
01676         markRow[iRow] = l - lSave;
01677         indexRowL[l] = iRow;
01678         elementL[l] = element[i];
01679         l++;
01680         //take out of row list
01681         CoinBigIndex start = startRow[iRow];
01682         CoinBigIndex end = start + numberInRow[iRow];
01683         CoinBigIndex where = start;
01684         
01685         while ( indexColumn[where] != iPivotColumn ) {
01686           where++;
01687         }                               /* endwhile */
01688 #if DEBUG_COIN
01689         if ( where >= end ) {
01690           abort (  );
01691         }
01692 #endif
01693         indexColumn[where] = indexColumn[end - 1];
01694         numberInRow[iRow]--;
01695         assert (numberInRow[iRow]>=0);
01696       }
01697     }
01698     assert (pivotRowPosition<endColumn);
01699     assert (indexRow[pivotRowPosition]==iPivotRow);
01700     CoinFactorizationDouble pivotElement = element[pivotRowPosition];
01701     CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
01702     
01703     pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
01704     pivotRowPosition++;
01705     for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01706       int iRow = indexRow[pivotRowPosition];
01707       
01708       markRow[iRow] = l - lSave;
01709       indexRowL[l] = iRow;
01710       elementL[l] = element[pivotRowPosition];
01711       l++;
01712       //take out of row list
01713       CoinBigIndex start = startRow[iRow];
01714       CoinBigIndex end = start + numberInRow[iRow];
01715       CoinBigIndex where = start;
01716       
01717       while ( indexColumn[where] != iPivotColumn ) {
01718         where++;
01719       }                         /* endwhile */
01720 #if DEBUG_COIN
01721       if ( where >= end ) {
01722         abort (  );
01723       }
01724 #endif
01725       indexColumn[where] = indexColumn[end - 1];
01726       numberInRow[iRow]--;
01727       assert (numberInRow[iRow]>=0);
01728     }
01729     markRow[iPivotRow] = FAC_SET;
01730     //compress pivot column (move pivot to front including saved)
01731     numberInColumn[iPivotColumn] = 0;
01732     //use end of L for temporary space
01733     int *indexL = &indexRowL[lSave];
01734     CoinFactorizationDouble *multipliersL = &elementL[lSave];
01735     
01736     //adjust
01737     int j;
01738     
01739     for ( j = 0; j < numberDoColumn; j++ ) {
01740       multipliersL[j] *= pivotMultiplier;
01741     }
01742     //zero out fill
01743     CoinBigIndex iErase;
01744     for ( iErase = 0; iErase < increment2 * numberDoRow;
01745           iErase++ ) {
01746       workArea2[iErase] = 0;
01747     }
01748     CoinBigIndex added = numberDoRow * numberDoColumn;
01749     unsigned int *temp2 = workArea2;
01750     int * nextColumn = nextColumn_.array();
01751     
01752     //pack down and move to work
01753     int jColumn;
01754     for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01755       int iColumn = saveColumn[jColumn];
01756       CoinBigIndex startColumnThis = startColumn[iColumn];
01757       CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
01758       int iRow = indexRow[startColumnThis];
01759       CoinFactorizationDouble value = element[startColumnThis];
01760       double largest;
01761       CoinBigIndex put = startColumnThis;
01762       CoinBigIndex positionLargest = -1;
01763       CoinFactorizationDouble thisPivotValue = 0.0;
01764       
01765       //compress column and find largest not updated
01766       bool checkLargest;
01767       int mark = markRow[iRow];
01768       
01769       if ( mark == FAC_UNSET ) {
01770         largest = fabs ( value );
01771         positionLargest = put;
01772         put++;
01773         checkLargest = false;
01774       } else {
01775         //need to find largest
01776         largest = 0.0;
01777         checkLargest = true;
01778         if ( mark != FAC_SET ) {
01779           //will be updated
01780           workArea[mark] = value;
01781           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01782           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01783           
01784           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
01785           added--;
01786         } else {
01787           thisPivotValue = value;
01788         }
01789       }
01790       CoinBigIndex i;
01791       for ( i = startColumnThis + 1; i < endColumn; i++ ) {
01792         iRow = indexRow[i];
01793         value = element[i];
01794         int mark = markRow[iRow];
01795         
01796         if ( mark == FAC_UNSET ) {
01797           //keep
01798           indexRow[put] = iRow;
01799           element[put] = value;
01800           if ( checkLargest ) {
01801             double absValue = fabs ( value );
01802             
01803             if ( absValue > largest ) {
01804               largest = absValue;
01805               positionLargest = put;
01806             }
01807           }
01808           put++;
01809         } else if ( mark != FAC_SET ) {
01810           //will be updated
01811           workArea[mark] = value;
01812           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01813           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01814           
01815           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
01816           added--;
01817         } else {
01818           thisPivotValue = value;
01819         }
01820       }
01821       //slot in pivot
01822       element[put] = element[startColumnThis];
01823       indexRow[put] = indexRow[startColumnThis];
01824       if ( positionLargest == startColumnThis ) {
01825         positionLargest = put;  //follow if was largest
01826       }
01827       put++;
01828       element[startColumnThis] = thisPivotValue;
01829       indexRow[startColumnThis] = iPivotRow;
01830       //clean up counts
01831       startColumnThis++;
01832       numberInColumn[iColumn] = put - startColumnThis;
01833       int * numberInColumnPlus = numberInColumnPlus_.array();
01834       numberInColumnPlus[iColumn]++;
01835       startColumn[iColumn]++;
01836       //how much space have we got
01837       int next = nextColumn[iColumn];
01838       CoinBigIndex space;
01839       
01840       space = startColumn[next] - put - numberInColumnPlus[next];
01841       //assume no zero elements
01842       if ( numberDoColumn > space ) {
01843         //getColumnSpace also moves fixed part
01844         if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
01845           goto BAD_PIVOT;
01846         }
01847         //redo starts
01848         positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
01849         startColumnThis = startColumn[iColumn];
01850         put = startColumnThis + numberInColumn[iColumn];
01851       }
01852       double tolerance = zeroTolerance_;
01853       
01854       int *nextCount = nextCount_.array();
01855       for ( j = 0; j < numberDoColumn; j++ ) {
01856         value = workArea[j] - thisPivotValue * multipliersL[j];
01857         double absValue = fabs ( value );
01858         
01859         if ( absValue > tolerance ) {
01860           workArea[j] = 0.0;
01861           element[put] = value;
01862           indexRow[put] = indexL[j];
01863           if ( absValue > largest ) {
01864             largest = absValue;
01865             positionLargest = put;
01866           }
01867           put++;
01868         } else {
01869           workArea[j] = 0.0;
01870           added--;
01871           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01872           int bit = j & COINFACTORIZATION_MASK_PER_INT;
01873           
01874           if ( temp2[word] & ( 1 << bit ) ) {
01875             //take out of row list
01876             iRow = indexL[j];
01877             CoinBigIndex start = startRow[iRow];
01878             CoinBigIndex end = start + numberInRow[iRow];
01879             CoinBigIndex where = start;
01880             
01881             while ( indexColumn[where] != iColumn ) {
01882               where++;
01883             }                   /* endwhile */
01884 #if DEBUG_COIN
01885             if ( where >= end ) {
01886               abort (  );
01887             }
01888 #endif
01889             indexColumn[where] = indexColumn[end - 1];
01890             numberInRow[iRow]--;
01891           } else {
01892             //make sure won't be added
01893             int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01894             int bit = j & COINFACTORIZATION_MASK_PER_INT;
01895             
01896             temp2[word] = temp2[word] | ( 1 << bit );   //say already in counts
01897           }
01898         }
01899       }
01900       numberInColumn[iColumn] = put - startColumnThis;
01901       //move largest
01902       if ( positionLargest >= 0 ) {
01903         value = element[positionLargest];
01904         iRow = indexRow[positionLargest];
01905         element[positionLargest] = element[startColumnThis];
01906         indexRow[positionLargest] = indexRow[startColumnThis];
01907         element[startColumnThis] = value;
01908         indexRow[startColumnThis] = iRow;
01909       }
01910       //linked list for column
01911       if ( nextCount[iColumn + numberRows_] != -2 ) {
01912         //modify linked list
01913         deleteLink ( iColumn + numberRows_ );
01914         addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01915       }
01916       temp2 += increment2;
01917     }
01918     //get space for row list
01919     unsigned int *putBase = workArea2;
01920     int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01921     int i = 0;
01922     
01923     // do linked lists and update counts
01924     while ( bigLoops ) {
01925       bigLoops--;
01926       int bit;
01927       for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01928         unsigned int *putThis = putBase;
01929         int iRow = indexL[i];
01930         
01931         //get space
01932         int number = 0;
01933         int jColumn;
01934         
01935         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01936           unsigned int test = *putThis;
01937           
01938           putThis += increment2;
01939           test = 1 - ( ( test >> bit ) & 1 );
01940           number += test;
01941         }
01942         int next = nextRow[iRow];
01943         CoinBigIndex space;
01944         
01945         space = startRow[next] - startRow[iRow];
01946         number += numberInRow[iRow];
01947         if ( space < number ) {
01948           if ( !getRowSpace ( iRow, number ) ) {
01949             goto BAD_PIVOT;
01950           }
01951         }
01952         // now do
01953         putThis = putBase;
01954         next = nextRow[iRow];
01955         number = numberInRow[iRow];
01956         CoinBigIndex end = startRow[iRow] + number;
01957         int saveIndex = indexColumn[startRow[next]];
01958         
01959         //add in
01960         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01961           unsigned int test = *putThis;
01962           
01963           putThis += increment2;
01964           test = 1 - ( ( test >> bit ) & 1 );
01965           indexColumn[end] = saveColumn[jColumn];
01966           end += test;
01967         }
01968         //put back next one in case zapped
01969         indexColumn[startRow[next]] = saveIndex;
01970         markRow[iRow] = FAC_UNSET;
01971         number = end - startRow[iRow];
01972         numberInRow[iRow] = number;
01973         deleteLink ( iRow );
01974         addLink ( iRow, number );
01975       }
01976       putBase++;
01977     }                           /* endwhile */
01978     int bit;
01979     
01980     for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
01981       unsigned int *putThis = putBase;
01982       int iRow = indexL[i];
01983       
01984       //get space
01985       int number = 0;
01986       int jColumn;
01987       
01988       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01989         unsigned int test = *putThis;
01990         
01991         putThis += increment2;
01992         test = 1 - ( ( test >> bit ) & 1 );
01993         number += test;
01994       }
01995       int next = nextRow[iRow];
01996       CoinBigIndex space;
01997       
01998       space = startRow[next] - startRow[iRow];
01999       number += numberInRow[iRow];
02000       if ( space < number ) {
02001         if ( !getRowSpace ( iRow, number ) ) {
02002           goto BAD_PIVOT;
02003         }
02004       }
02005       // now do
02006       putThis = putBase;
02007       next = nextRow[iRow];
02008       number = numberInRow[iRow];
02009       CoinBigIndex end = startRow[iRow] + number;
02010       int saveIndex;
02011       
02012       saveIndex = indexColumn[startRow[next]];
02013       
02014       //add in
02015       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
02016         unsigned int test = *putThis;
02017         
02018         putThis += increment2;
02019         test = 1 - ( ( test >> bit ) & 1 );
02020         
02021         indexColumn[end] = saveColumn[jColumn];
02022         end += test;
02023       }
02024       indexColumn[startRow[next]] = saveIndex;
02025       markRow[iRow] = FAC_UNSET;
02026       number = end - startRow[iRow];
02027       numberInRow[iRow] = number;
02028       deleteLink ( iRow );
02029       addLink ( iRow, number );
02030     }
02031     markRow[iPivotRow] = FAC_UNSET;
02032     //modify linked list for pivots
02033     deleteLink ( iPivotRow );
02034     deleteLink ( iPivotColumn + numberRows_ );
02035     totalElements_ += added;
02036     goodPivot= true;
02037     // **** UGLY UGLY UGLY
02038   }
02039  BAD_PIVOT:
02040 
02041   ;
02042 }
02043 #undef FAC_UNSET
02044 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 28 Aug 2016 for CoinUtils by  doxygen 1.6.1