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

Generated on 12 Mar 2015 for Dip-All by  doxygen 1.6.1