/home/coin/SVN-release/CoinUtils-2.3.1/CoinUtils/src/CoinFactorization.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2002, International Business Machines
00002 // Corporation and others.  All Rights Reserved.
00003 
00004 /* 
00005    Authors
00006    
00007    John Forrest
00008 
00009  */
00010 #ifndef CoinFactorization_H
00011 #define CoinFactorization_H
00012 
00013 #include <iostream>
00014 #include <string>
00015 #include <cassert>
00016 #include "CoinFinite.hpp"
00017 #include "CoinIndexedVector.hpp"
00018 class CoinPackedMatrix;
00045 class CoinFactorization {
00046    friend void CoinFactorizationUnitTest( const std::string & mpsDir );
00047 
00048 public:
00049 
00052 
00053     CoinFactorization (  );
00055   CoinFactorization ( const CoinFactorization &other);
00056 
00058    ~CoinFactorization (  );
00060   void almostDestructor();
00062   void show_self (  ) const;
00064   int saveFactorization (const char * file  ) const;
00068   int restoreFactorization (const char * file  , bool factor=false) ;
00070   void sort (  ) const;
00072     CoinFactorization & operator = ( const CoinFactorization & other );
00074 
00084   int factorize ( const CoinPackedMatrix & matrix, 
00085                   int rowIsBasic[], int columnIsBasic[] , 
00086                   double areaFactor = 0.0 );
00097   int factorize ( int numberRows,
00098                   int numberColumns,
00099                   CoinBigIndex numberElements,
00100                   CoinBigIndex maximumL,
00101                   CoinBigIndex maximumU,
00102                   const int indicesRow[],
00103                   const int indicesColumn[], const double elements[] ,
00104                   int permutation[],
00105                   double areaFactor = 0.0);
00110   int factorizePart1 ( int numberRows,
00111                        int numberColumns,
00112                        CoinBigIndex estimateNumberElements,
00113                        int * indicesRow[],
00114                        int * indicesColumn[],
00115                        double * elements[],
00116                   double areaFactor = 0.0);
00123   int factorizePart2 (int permutation[],int exactNumberElements);
00125   double conditionNumber() const;
00126   
00128 
00131 
00132   inline int status (  ) const {
00133     return status_;
00134   }
00136   inline void setStatus (  int value)
00137   {  status_=value;  }
00139   inline int pivots (  ) const {
00140     return numberPivots_;
00141   }
00143   inline void setPivots (  int value ) 
00144   { numberPivots_=value; }
00146   inline int *permute (  ) const {
00147     return permute_.array();
00148   }
00150   inline int *pivotColumn (  ) const {
00151     return pivotColumn_.array();
00152   }
00154   inline double *pivotRegion (  ) const {
00155     return pivotRegion_.array();
00156   }
00158   inline int *permuteBack (  ) const {
00159     return permuteBack_.array();
00160   }
00162   inline int *pivotColumnBack (  ) const {
00163     return pivotColumnBack_.array();
00164   }
00166   inline CoinBigIndex * startRowL() const
00167   { return startRowL_.array();}
00168 
00170   inline CoinBigIndex * startColumnL() const
00171   { return startColumnL_.array();}
00172 
00174   inline int * indexColumnL() const
00175   { return indexColumnL_.array();}
00176 
00178   inline int * indexRowL() const
00179   { return indexRowL_.array();}
00180 
00182   inline double * elementByRowL() const
00183   { return elementByRowL_.array();}
00184 
00186   inline int numberRowsExtra (  ) const {
00187     return numberRowsExtra_;
00188   }
00190   inline void setNumberRows(int value)
00191   { numberRows_ = value; }
00193   inline int numberRows (  ) const {
00194     return numberRows_;
00195   }
00197   inline CoinBigIndex numberL() const
00198   { return numberL_;}
00199 
00201   inline CoinBigIndex baseL() const
00202   { return baseL_;}
00204   inline int maximumRowsExtra (  ) const {
00205     return maximumRowsExtra_;
00206   }
00208   inline int numberColumns (  ) const {
00209     return numberColumns_;
00210   }
00212   inline int numberElements (  ) const {
00213     return totalElements_;
00214   }
00216   inline int numberForrestTomlin (  ) const {
00217     return numberInColumn_.array()[numberColumnsExtra_];
00218   }
00220   inline int numberGoodColumns (  ) const {
00221     return numberGoodU_;
00222   }
00224   inline double areaFactor (  ) const {
00225     return areaFactor_;
00226   }
00227   inline void areaFactor ( double value ) {
00228     areaFactor_=value;
00229   }
00231   double adjustedAreaFactor() const;
00233   inline void relaxAccuracyCheck(double value)
00234   { relaxCheck_ = value;}
00235   inline double getAccuracyCheck() const
00236   { return relaxCheck_;}
00238   inline int messageLevel (  ) const {
00239     return messageLevel_ ;
00240   }
00241   void messageLevel (  int value );
00243   inline int maximumPivots (  ) const {
00244     return maximumPivots_ ;
00245   }
00246   void maximumPivots (  int value );
00247 
00249   inline int denseThreshold() const 
00250     { return denseThreshold_;}
00252   inline void setDenseThreshold(int value)
00253     { denseThreshold_ = value;}
00255   inline double pivotTolerance (  ) const {
00256     return pivotTolerance_ ;
00257   }
00258   void pivotTolerance (  double value );
00260   inline double zeroTolerance (  ) const {
00261     return zeroTolerance_ ;
00262   }
00263   void zeroTolerance (  double value );
00264 #ifndef COIN_FAST_CODE
00266   inline double slackValue (  ) const {
00267     return slackValue_ ;
00268   }
00269   void slackValue (  double value );
00270 #endif
00272   double maximumCoefficient() const;
00274   inline bool forrestTomlin() const
00275   { return doForrestTomlin_;}
00276   inline void setForrestTomlin(bool value)
00277   { doForrestTomlin_=value;}
00279   inline bool spaceForForrestTomlin() const
00280   {
00281     CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
00282     CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
00283     return (space>=0)&&doForrestTomlin_;
00284   }
00286 
00289 
00291   inline int numberDense() const
00292   { return numberDense_;}
00293 
00295   inline CoinBigIndex numberElementsU (  ) const {
00296     return lengthU_;
00297   }
00299   inline void setNumberElementsU(CoinBigIndex value)
00300   { lengthU_ = value; }
00302   inline CoinBigIndex lengthAreaU (  ) const {
00303     return lengthAreaU_;
00304   }
00306   inline CoinBigIndex numberElementsL (  ) const {
00307     return lengthL_;
00308   }
00310   inline CoinBigIndex lengthAreaL (  ) const {
00311     return lengthAreaL_;
00312   }
00314   inline CoinBigIndex numberElementsR (  ) const {
00315     return lengthR_;
00316   }
00318   inline CoinBigIndex numberCompressions() const
00319   { return numberCompressions_;}
00321   inline int * numberInRow() const
00322   { return numberInRow_.array();}
00324   inline int * numberInColumn() const
00325   { return numberInColumn_.array();}
00327   inline double * elementU() const
00328   { return elementU_.array();}
00330   inline int * indexRowU() const
00331   { return indexRowU_.array();}
00333   inline CoinBigIndex * startColumnU() const
00334   { return startColumnU_.array();}
00336   inline int maximumColumnsExtra()
00337   { return maximumColumnsExtra_;}
00341   inline int biasLU() const
00342   { return biasLU_;}
00343   inline void setBiasLU(int value)
00344   { biasLU_=value;}
00350   inline int persistenceFlag() const
00351   { return persistenceFlag_;}
00352   void setPersistenceFlag(int value);
00354 
00357 
00365   int replaceColumn ( CoinIndexedVector * regionSparse,
00366                       int pivotRow,
00367                       double pivotCheck ,
00368                       bool checkBeforeModifying=false);
00370 
00380   int updateColumnFT ( CoinIndexedVector * regionSparse,
00381                        CoinIndexedVector * regionSparse2);
00384   int updateColumn ( CoinIndexedVector * regionSparse,
00385                      CoinIndexedVector * regionSparse2,
00386                      bool noPermute=false) const;
00392   int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
00393                            CoinIndexedVector * regionSparse2,
00394                            CoinIndexedVector * regionSparse3,
00395                            bool noPermuteRegion3=false) ;
00400   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00401                               CoinIndexedVector * regionSparse2) const;
00403   void goSparse();
00405   inline int sparseThreshold ( ) const
00406   { return sparseThreshold_;}
00408   void sparseThreshold ( int value );
00410 
00411 
00415 
00416   inline void clearArrays()
00417   { gutsOfDestructor();}
00419 
00422 
00425   int add ( CoinBigIndex numberElements,
00426                int indicesRow[],
00427                int indicesColumn[], double elements[] );
00428 
00431   int addColumn ( CoinBigIndex numberElements,
00432                      int indicesRow[], double elements[] );
00433 
00436   int addRow ( CoinBigIndex numberElements,
00437                   int indicesColumn[], double elements[] );
00438 
00440   int deleteColumn ( int Row );
00442   int deleteRow ( int Row );
00443 
00447   int replaceRow ( int whichRow, int numberElements,
00448                       const int indicesColumn[], const double elements[] );
00450   void emptyRows(int numberToEmpty, const int which[]);
00452 
00453 
00454   void checkSparse();
00456   inline bool collectStatistics() const
00457   { return collectStatistics_;}
00459   inline void setCollectStatistics(bool onOff) const
00460   { collectStatistics_ = onOff;}
00462   void gutsOfDestructor(int type=1);
00464   void gutsOfInitialize(int type);
00465   void gutsOfCopy(const CoinFactorization &other);
00466 
00468   void resetStatistics();
00469 
00470 
00472 
00474 
00475   void getAreas ( int numberRows,
00476                   int numberColumns,
00477                   CoinBigIndex maximumL,
00478                   CoinBigIndex maximumU );
00479 
00482   void preProcess ( int state,
00483                     int possibleDuplicates = -1 );
00485   int factor (  );
00486 protected:
00489   int factorSparse (  );
00492   int factorDense (  );
00493 
00495   bool pivotOneOtherRow ( int pivotRow,
00496                           int pivotColumn );
00498   bool pivotRowSingleton ( int pivotRow,
00499                            int pivotColumn );
00501   bool pivotColumnSingleton ( int pivotRow,
00502                               int pivotColumn );
00503 
00508   bool getColumnSpace ( int iColumn,
00509                         int extraNeeded );
00510 
00514   bool getColumnSpaceIterateR ( int iColumn, double value,
00515                                int iRow);
00521   CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00522                                int iRow);
00526   bool getRowSpace ( int iRow, int extraNeeded );
00527 
00531   bool getRowSpaceIterate ( int iRow,
00532                             int extraNeeded );
00534   void checkConsistency (  );
00536   inline void addLink ( int index, int count ) {
00537     int *nextCount = nextCount_.array();
00538     int *firstCount = firstCount_.array();
00539     int *lastCount = lastCount_.array();
00540     int next = firstCount[count];
00541       lastCount[index] = -2 - count;
00542     if ( next < 0 ) {
00543       //first with that count
00544       firstCount[count] = index;
00545       nextCount[index] = -1;
00546     } else {
00547       firstCount[count] = index;
00548       nextCount[index] = next;
00549       lastCount[next] = index;
00550   }}
00552   inline void deleteLink ( int index ) {
00553     int *nextCount = nextCount_.array();
00554     int *firstCount = firstCount_.array();
00555     int *lastCount = lastCount_.array();
00556     int next = nextCount[index];
00557     int last = lastCount[index];
00558     if ( last >= 0 ) {
00559       nextCount[last] = next;
00560     } else {
00561       int count = -last - 2;
00562 
00563       firstCount[count] = next;
00564     }
00565     if ( next >= 0 ) {
00566       lastCount[next] = last;
00567     }
00568     nextCount[index] = -2;
00569     lastCount[index] = -2;
00570     return;
00571   }
00573   void separateLinks(int count,bool rowsFirst);
00575   void cleanup (  );
00576 
00578   void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00580   void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00582   void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00584   void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00585 
00587   void updateColumnR ( CoinIndexedVector * region ) const;
00590   void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00591 
00593   void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00594 
00596   void updateColumnUSparse ( CoinIndexedVector * regionSparse, 
00597                              int * indexIn) const;
00599   void updateColumnUSparsish ( CoinIndexedVector * regionSparse, 
00600                                int * indexIn) const;
00602   int updateColumnUDensish ( double * COIN_RESTRICT region, 
00603                              int * COIN_RESTRICT regionIndex) const;
00605   void updateTwoColumnsUDensish (
00606                                  int & numberNonZero1,
00607                                  double * COIN_RESTRICT region1, 
00608                                  int * COIN_RESTRICT index1,
00609                                  int & numberNonZero2,
00610                                  double * COIN_RESTRICT region2, 
00611                                  int * COIN_RESTRICT index2) const;
00613   void updateColumnPFI ( CoinIndexedVector * regionSparse) const; 
00615   void permuteBack ( CoinIndexedVector * regionSparse, 
00616                      CoinIndexedVector * outVector) const;
00617 
00619   void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00622   void updateColumnTransposeU ( CoinIndexedVector * region,
00623                                 int smallestIndex) const;
00626   void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00627                                         int smallestIndex) const;
00630   void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00631                                        int smallestIndex) const;
00634   void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00635 
00637   void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00639   void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00641   void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00642 
00644   void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00646   void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00648   void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00650   void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00652   void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00653 public:
00658   int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00659                          int pivotRow, double alpha);
00660 protected:
00663   int checkPivot(double saveFromU, double oldPivot) const;
00664   /********************************* START LARGE TEMPLATE ********/
00665 #ifdef INT_IS_8
00666 #define COINFACTORIZATION_BITS_PER_INT 64
00667 #define COINFACTORIZATION_SHIFT_PER_INT 6
00668 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00669 #else
00670 #define COINFACTORIZATION_BITS_PER_INT 32
00671 #define COINFACTORIZATION_SHIFT_PER_INT 5
00672 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00673 #endif
00674   template <class T>  inline bool
00675   pivot ( int pivotRow,
00676           int pivotColumn,
00677           CoinBigIndex pivotRowPosition,
00678           CoinBigIndex pivotColumnPosition,
00679           double work[],
00680           unsigned int workArea2[],
00681           int increment,
00682           int increment2,
00683           T markRow[] ,
00684           int largeInteger)
00685 {
00686   int *indexColumnU = indexColumnU_.array();
00687   CoinBigIndex *startColumnU = startColumnU_.array();
00688   int *numberInColumn = numberInColumn_.array();
00689   double *elementU = elementU_.array();
00690   int *indexRowU = indexRowU_.array();
00691   CoinBigIndex *startRowU = startRowU_.array();
00692   int *numberInRow = numberInRow_.array();
00693   double *elementL = elementL_.array();
00694   int *indexRowL = indexRowL_.array();
00695   int *saveColumn = saveColumn_.array();
00696   int *nextRow = nextRow_.array();
00697   int *lastRow = lastRow_.array() ;
00698 
00699   //store pivot columns (so can easily compress)
00700   int numberInPivotRow = numberInRow[pivotRow] - 1;
00701   CoinBigIndex startColumn = startColumnU[pivotColumn];
00702   int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00703   CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00704   int put = 0;
00705   CoinBigIndex startRow = startRowU[pivotRow];
00706   CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00707 
00708   if ( pivotColumnPosition < 0 ) {
00709     for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00710       int iColumn = indexColumnU[pivotColumnPosition];
00711       if ( iColumn != pivotColumn ) {
00712         saveColumn[put++] = iColumn;
00713       } else {
00714         break;
00715       }
00716     }
00717   } else {
00718     for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00719       saveColumn[put++] = indexColumnU[i];
00720     }
00721   }
00722   assert (pivotColumnPosition<endRow);
00723   assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00724   pivotColumnPosition++;
00725   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00726     saveColumn[put++] = indexColumnU[pivotColumnPosition];
00727   }
00728   //take out this bit of indexColumnU
00729   int next = nextRow[pivotRow];
00730   int last = lastRow[pivotRow];
00731 
00732   nextRow[last] = next;
00733   lastRow[next] = last;
00734   nextRow[pivotRow] = numberGoodU_;     //use for permute
00735   lastRow[pivotRow] = -2;
00736   numberInRow[pivotRow] = 0;
00737   //store column in L, compress in U and take column out
00738   CoinBigIndex l = lengthL_;
00739 
00740   if ( l + numberInPivotColumn > lengthAreaL_ ) {
00741     //need more memory
00742     printf("more memory needed in middle of invert\n");
00743     return false;
00744   }
00745   //l+=currentAreaL_->elementByColumn-elementL;
00746   CoinBigIndex lSave = l;
00747 
00748   pivotRowL_.array()[numberGoodL_] = pivotRow;
00749   CoinBigIndex * startColumnL = startColumnL_.array();
00750   startColumnL[numberGoodL_] = l;       //for luck and first time
00751   numberGoodL_++;
00752   startColumnL[numberGoodL_] = l + numberInPivotColumn;
00753   lengthL_ += numberInPivotColumn;
00754   if ( pivotRowPosition < 0 ) {
00755     for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00756       int iRow = indexRowU[pivotRowPosition];
00757       if ( iRow != pivotRow ) {
00758         indexRowL[l] = iRow;
00759         elementL[l] = elementU[pivotRowPosition];
00760         markRow[iRow] = l - lSave;
00761         l++;
00762         //take out of row list
00763         CoinBigIndex start = startRowU[iRow];
00764         CoinBigIndex end = start + numberInRow[iRow];
00765         CoinBigIndex where = start;
00766 
00767         while ( indexColumnU[where] != pivotColumn ) {
00768           where++;
00769         }                       /* endwhile */
00770 #if DEBUG_COIN
00771         if ( where >= end ) {
00772           abort (  );
00773         }
00774 #endif
00775         indexColumnU[where] = indexColumnU[end - 1];
00776         numberInRow[iRow]--;
00777       } else {
00778         break;
00779       }
00780     }
00781   } else {
00782     CoinBigIndex i;
00783 
00784     for ( i = startColumn; i < pivotRowPosition; i++ ) {
00785       int iRow = indexRowU[i];
00786 
00787       markRow[iRow] = l - lSave;
00788       indexRowL[l] = iRow;
00789       elementL[l] = elementU[i];
00790       l++;
00791       //take out of row list
00792       CoinBigIndex start = startRowU[iRow];
00793       CoinBigIndex end = start + numberInRow[iRow];
00794       CoinBigIndex where = start;
00795 
00796       while ( indexColumnU[where] != pivotColumn ) {
00797         where++;
00798       }                         /* endwhile */
00799 #if DEBUG_COIN
00800       if ( where >= end ) {
00801         abort (  );
00802       }
00803 #endif
00804       indexColumnU[where] = indexColumnU[end - 1];
00805       numberInRow[iRow]--;
00806       assert (numberInRow[iRow]>=0);
00807     }
00808   }
00809   assert (pivotRowPosition<endColumn);
00810   assert (indexRowU[pivotRowPosition]==pivotRow);
00811   double pivotElement = elementU[pivotRowPosition];
00812   double pivotMultiplier = 1.0 / pivotElement;
00813 
00814   pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00815   pivotRowPosition++;
00816   for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00817     int iRow = indexRowU[pivotRowPosition];
00818     
00819     markRow[iRow] = l - lSave;
00820     indexRowL[l] = iRow;
00821     elementL[l] = elementU[pivotRowPosition];
00822     l++;
00823     //take out of row list
00824     CoinBigIndex start = startRowU[iRow];
00825     CoinBigIndex end = start + numberInRow[iRow];
00826     CoinBigIndex where = start;
00827     
00828     while ( indexColumnU[where] != pivotColumn ) {
00829       where++;
00830     }                           /* endwhile */
00831 #if DEBUG_COIN
00832     if ( where >= end ) {
00833       abort (  );
00834     }
00835 #endif
00836     indexColumnU[where] = indexColumnU[end - 1];
00837     numberInRow[iRow]--;
00838     assert (numberInRow[iRow]>=0);
00839   }
00840   markRow[pivotRow] = largeInteger;
00841   //compress pivot column (move pivot to front including saved)
00842   numberInColumn[pivotColumn] = 0;
00843   //use end of L for temporary space
00844   int *indexL = &indexRowL[lSave];
00845   double *multipliersL = &elementL[lSave];
00846 
00847   //adjust
00848   int j;
00849 
00850   for ( j = 0; j < numberInPivotColumn; j++ ) {
00851     multipliersL[j] *= pivotMultiplier;
00852   }
00853   //zero out fill
00854   CoinBigIndex iErase;
00855   for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00856         iErase++ ) {
00857     workArea2[iErase] = 0;
00858   }
00859   CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00860   unsigned int *temp2 = workArea2;
00861   int * nextColumn = nextColumn_.array();
00862 
00863   //pack down and move to work
00864   int jColumn;
00865   for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
00866     int iColumn = saveColumn[jColumn];
00867     CoinBigIndex startColumn = startColumnU[iColumn];
00868     CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
00869     int iRow = indexRowU[startColumn];
00870     double value = elementU[startColumn];
00871     double largest;
00872     CoinBigIndex put = startColumn;
00873     CoinBigIndex positionLargest = -1;
00874     double thisPivotValue = 0.0;
00875 
00876     //compress column and find largest not updated
00877     bool checkLargest;
00878     int mark = markRow[iRow];
00879 
00880     if ( mark < 0 ) {
00881       largest = fabs ( value );
00882       positionLargest = put;
00883       put++;
00884       checkLargest = false;
00885     } else {
00886       //need to find largest
00887       largest = 0.0;
00888       checkLargest = true;
00889       if ( mark != largeInteger ) {
00890         //will be updated
00891         work[mark] = value;
00892         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00893         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00894 
00895         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
00896         added--;
00897       } else {
00898         thisPivotValue = value;
00899       }
00900     }
00901     CoinBigIndex i;
00902     for ( i = startColumn + 1; i < endColumn; i++ ) {
00903       iRow = indexRowU[i];
00904       value = elementU[i];
00905       int mark = markRow[iRow];
00906 
00907       if ( mark < 0 ) {
00908         //keep
00909         indexRowU[put] = iRow;
00910         elementU[put] = value;
00911         if ( checkLargest ) {
00912           double absValue = fabs ( value );
00913 
00914           if ( absValue > largest ) {
00915             largest = absValue;
00916             positionLargest = put;
00917           }
00918         }
00919         put++;
00920       } else if ( mark != largeInteger ) {
00921         //will be updated
00922         work[mark] = value;
00923         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00924         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00925 
00926         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
00927         added--;
00928       } else {
00929         thisPivotValue = value;
00930       }
00931     }
00932     //slot in pivot
00933     elementU[put] = elementU[startColumn];
00934     indexRowU[put] = indexRowU[startColumn];
00935     if ( positionLargest == startColumn ) {
00936       positionLargest = put;    //follow if was largest
00937     }
00938     put++;
00939     elementU[startColumn] = thisPivotValue;
00940     indexRowU[startColumn] = pivotRow;
00941     //clean up counts
00942     startColumn++;
00943     numberInColumn[iColumn] = put - startColumn;
00944     int * numberInColumnPlus = numberInColumnPlus_.array();
00945     numberInColumnPlus[iColumn]++;
00946     startColumnU[iColumn]++;
00947     //how much space have we got
00948     int next = nextColumn[iColumn];
00949     CoinBigIndex space;
00950 
00951     space = startColumnU[next] - put - numberInColumnPlus[next];
00952     //assume no zero elements
00953     if ( numberInPivotColumn > space ) {
00954       //getColumnSpace also moves fixed part
00955       if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
00956         return false;
00957       }
00958       //redo starts
00959       positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
00960       startColumn = startColumnU[iColumn];
00961       put = startColumn + numberInColumn[iColumn];
00962     }
00963     double tolerance = zeroTolerance_;
00964 
00965     int *nextCount = nextCount_.array();
00966     for ( j = 0; j < numberInPivotColumn; j++ ) {
00967       value = work[j] - thisPivotValue * multipliersL[j];
00968       double absValue = fabs ( value );
00969 
00970       if ( absValue > tolerance ) {
00971         work[j] = 0.0;
00972         elementU[put] = value;
00973         indexRowU[put] = indexL[j];
00974         if ( absValue > largest ) {
00975           largest = absValue;
00976           positionLargest = put;
00977         }
00978         put++;
00979       } else {
00980         work[j] = 0.0;
00981         added--;
00982         int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
00983         int bit = j & COINFACTORIZATION_MASK_PER_INT;
00984 
00985         if ( temp2[word] & ( 1 << bit ) ) {
00986           //take out of row list
00987           iRow = indexL[j];
00988           CoinBigIndex start = startRowU[iRow];
00989           CoinBigIndex end = start + numberInRow[iRow];
00990           CoinBigIndex where = start;
00991 
00992           while ( indexColumnU[where] != iColumn ) {
00993             where++;
00994           }                     /* endwhile */
00995 #if DEBUG_COIN
00996           if ( where >= end ) {
00997             abort (  );
00998           }
00999 #endif
01000           indexColumnU[where] = indexColumnU[end - 1];
01001           numberInRow[iRow]--;
01002         } else {
01003           //make sure won't be added
01004           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01005           int bit = j & COINFACTORIZATION_MASK_PER_INT;
01006 
01007           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
01008         }
01009       }
01010     }
01011     numberInColumn[iColumn] = put - startColumn;
01012     //move largest
01013     if ( positionLargest >= 0 ) {
01014       value = elementU[positionLargest];
01015       iRow = indexRowU[positionLargest];
01016       elementU[positionLargest] = elementU[startColumn];
01017       indexRowU[positionLargest] = indexRowU[startColumn];
01018       elementU[startColumn] = value;
01019       indexRowU[startColumn] = iRow;
01020     }
01021     //linked list for column
01022     if ( nextCount[iColumn + numberRows_] != -2 ) {
01023       //modify linked list
01024       deleteLink ( iColumn + numberRows_ );
01025       addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01026     }
01027     temp2 += increment2;
01028   }
01029   //get space for row list
01030   unsigned int *putBase = workArea2;
01031   int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01032   int i = 0;
01033 
01034   // do linked lists and update counts
01035   while ( bigLoops ) {
01036     bigLoops--;
01037     int bit;
01038     for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01039       unsigned int *putThis = putBase;
01040       int iRow = indexL[i];
01041 
01042       //get space
01043       int number = 0;
01044       int jColumn;
01045 
01046       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01047         unsigned int test = *putThis;
01048 
01049         putThis += increment2;
01050         test = 1 - ( ( test >> bit ) & 1 );
01051         number += test;
01052       }
01053       int next = nextRow[iRow];
01054       CoinBigIndex space;
01055 
01056       space = startRowU[next] - startRowU[iRow];
01057       number += numberInRow[iRow];
01058       if ( space < number ) {
01059         if ( !getRowSpace ( iRow, number ) ) {
01060           return false;
01061         }
01062       }
01063       // now do
01064       putThis = putBase;
01065       next = nextRow[iRow];
01066       number = numberInRow[iRow];
01067       CoinBigIndex end = startRowU[iRow] + number;
01068       int saveIndex = indexColumnU[startRowU[next]];
01069 
01070       //add in
01071       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01072         unsigned int test = *putThis;
01073 
01074         putThis += increment2;
01075         test = 1 - ( ( test >> bit ) & 1 );
01076         indexColumnU[end] = saveColumn[jColumn];
01077         end += test;
01078       }
01079       //put back next one in case zapped
01080       indexColumnU[startRowU[next]] = saveIndex;
01081       markRow[iRow] = -1;
01082       number = end - startRowU[iRow];
01083       numberInRow[iRow] = number;
01084       deleteLink ( iRow );
01085       addLink ( iRow, number );
01086     }
01087     putBase++;
01088   }                             /* endwhile */
01089   int bit;
01090 
01091   for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01092     unsigned int *putThis = putBase;
01093     int iRow = indexL[i];
01094 
01095     //get space
01096     int number = 0;
01097     int jColumn;
01098 
01099     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01100       unsigned int test = *putThis;
01101 
01102       putThis += increment2;
01103       test = 1 - ( ( test >> bit ) & 1 );
01104       number += test;
01105     }
01106     int next = nextRow[iRow];
01107     CoinBigIndex space;
01108 
01109     space = startRowU[next] - startRowU[iRow];
01110     number += numberInRow[iRow];
01111     if ( space < number ) {
01112       if ( !getRowSpace ( iRow, number ) ) {
01113         return false;
01114       }
01115     }
01116     // now do
01117     putThis = putBase;
01118     next = nextRow[iRow];
01119     number = numberInRow[iRow];
01120     CoinBigIndex end = startRowU[iRow] + number;
01121     int saveIndex;
01122 
01123     saveIndex = indexColumnU[startRowU[next]];
01124 
01125     //add in
01126     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01127       unsigned int test = *putThis;
01128 
01129       putThis += increment2;
01130       test = 1 - ( ( test >> bit ) & 1 );
01131 
01132       indexColumnU[end] = saveColumn[jColumn];
01133       end += test;
01134     }
01135     indexColumnU[startRowU[next]] = saveIndex;
01136     markRow[iRow] = -1;
01137     number = end - startRowU[iRow];
01138     numberInRow[iRow] = number;
01139     deleteLink ( iRow );
01140     addLink ( iRow, number );
01141   }
01142   markRow[pivotRow] = -1;
01143   //modify linked list for pivots
01144   deleteLink ( pivotRow );
01145   deleteLink ( pivotColumn + numberRows_ );
01146   totalElements_ += added;
01147   return true;
01148 }
01149 
01150   /********************************* END LARGE TEMPLATE ********/
01152 
01153 protected:
01154 
01157 
01158   double pivotTolerance_;
01160   double zeroTolerance_;
01161 #ifndef COIN_FAST_CODE
01163   double slackValue_;
01164 #else
01165 #ifndef slackValue_
01166 #define slackValue_ -1.0
01167 #endif
01168 #endif
01170   double areaFactor_;
01172   double relaxCheck_;
01174   int numberRows_;
01176   int numberRowsExtra_;
01178   int maximumRowsExtra_;
01180   int numberColumns_;
01182   int numberColumnsExtra_;
01184   int maximumColumnsExtra_;
01186   int numberGoodU_;
01188   int numberGoodL_;
01190   int maximumPivots_;
01192   int numberPivots_;
01195   CoinBigIndex totalElements_;
01197   CoinBigIndex factorElements_;
01199   CoinIntArrayWithLength pivotColumn_;
01201   CoinIntArrayWithLength permute_;
01203   CoinIntArrayWithLength permuteBack_;
01205   CoinIntArrayWithLength pivotColumnBack_;
01207   int status_;
01208 
01213   //int increasingRows_;
01214 
01216   int numberTrials_;
01218   CoinBigIndexArrayWithLength startRowU_;
01219 
01221   CoinIntArrayWithLength numberInRow_;
01222 
01224   CoinIntArrayWithLength numberInColumn_;
01225 
01227   CoinIntArrayWithLength numberInColumnPlus_;
01228 
01231   CoinIntArrayWithLength firstCount_;
01232 
01234   CoinIntArrayWithLength nextCount_;
01235 
01237   CoinIntArrayWithLength lastCount_;
01238 
01240   CoinIntArrayWithLength nextColumn_;
01241 
01243   CoinIntArrayWithLength lastColumn_;
01244 
01246   CoinIntArrayWithLength nextRow_;
01247 
01249   CoinIntArrayWithLength lastRow_;
01250 
01252   CoinIntArrayWithLength saveColumn_;
01253 
01255   CoinIntArrayWithLength markRow_;
01256 
01258   int messageLevel_;
01259 
01261   int biggerDimension_;
01262 
01264   CoinIntArrayWithLength indexColumnU_;
01265 
01267   CoinIntArrayWithLength pivotRowL_;
01268 
01270   CoinDoubleArrayWithLength pivotRegion_;
01271 
01273   int numberSlacks_;
01274 
01276   int numberU_;
01277 
01279   CoinBigIndex maximumU_;
01280 
01282   //int baseU_;
01283 
01285   CoinBigIndex lengthU_;
01286 
01288   CoinBigIndex lengthAreaU_;
01289 
01291   CoinDoubleArrayWithLength elementU_;
01292 
01294   CoinIntArrayWithLength indexRowU_;
01295 
01297   CoinBigIndexArrayWithLength startColumnU_;
01298 
01300   CoinBigIndexArrayWithLength convertRowToColumnU_;
01301 
01303   CoinBigIndex numberL_;
01304 
01306   CoinBigIndex baseL_;
01307 
01309   CoinBigIndex lengthL_;
01310 
01312   CoinBigIndex lengthAreaL_;
01313 
01315   CoinDoubleArrayWithLength elementL_;
01316 
01318   CoinIntArrayWithLength indexRowL_;
01319 
01321   CoinBigIndexArrayWithLength startColumnL_;
01322 
01324   bool doForrestTomlin_;
01325 
01327   int numberR_;
01328 
01330   CoinBigIndex lengthR_;
01331 
01333   CoinBigIndex lengthAreaR_;
01334 
01336   double *elementR_;
01337 
01339   int *indexRowR_;
01340 
01342   CoinBigIndexArrayWithLength startColumnR_;
01343 
01345   double  * denseArea_;
01346 
01348   int * densePermute_;
01349 
01351   int numberDense_;
01352 
01354   int denseThreshold_;
01355 
01357   CoinDoubleArrayWithLength workArea_;
01358 
01360   CoinUnsignedIntArrayWithLength workArea2_;
01361 
01363   CoinBigIndex numberCompressions_;
01364 
01366   mutable double ftranCountInput_;
01367   mutable double ftranCountAfterL_;
01368   mutable double ftranCountAfterR_;
01369   mutable double ftranCountAfterU_;
01370   mutable double btranCountInput_;
01371   mutable double btranCountAfterU_;
01372   mutable double btranCountAfterR_;
01373   mutable double btranCountAfterL_;
01374 
01376   mutable int numberFtranCounts_;
01377   mutable int numberBtranCounts_;
01378 
01380   double ftranAverageAfterL_;
01381   double ftranAverageAfterR_;
01382   double ftranAverageAfterU_;
01383   double btranAverageAfterU_;
01384   double btranAverageAfterR_;
01385   double btranAverageAfterL_;
01386 
01388   mutable bool collectStatistics_;
01389 
01391   int sparseThreshold_;
01392 
01394   int sparseThreshold2_;
01395 
01397   CoinBigIndexArrayWithLength startRowL_;
01398 
01400   CoinIntArrayWithLength indexColumnL_;
01401 
01403   CoinDoubleArrayWithLength elementByRowL_;
01404 
01406   mutable CoinIntArrayWithLength sparse_;
01410   int biasLU_;
01416   int persistenceFlag_;
01418 };
01419 // Dense coding
01420 #ifdef COIN_HAS_LAPACK
01421 #define DENSE_CODE 1
01422 /* Type of Fortran integer translated into C */
01423 #ifndef ipfint
01424 //typedef ipfint FORTRAN_INTEGER_TYPE ;
01425 typedef int ipfint;
01426 typedef const int cipfint;
01427 #endif
01428 #endif
01429 #endif

Generated on Thu Jul 17 03:01:06 2008 by  doxygen 1.4.7