/home/coin/SVN-release/Bcp-1.2.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 <cstdio>
00017 #include "CoinFinite.hpp"
00018 #include "CoinIndexedVector.hpp"
00019 class CoinPackedMatrix;
00046 class CoinFactorization {
00047    friend void CoinFactorizationUnitTest( const std::string & mpsDir );
00048 
00049 public:
00050 
00053 
00054     CoinFactorization (  );
00056   CoinFactorization ( const CoinFactorization &other);
00057 
00059    ~CoinFactorization (  );
00061   void almostDestructor();
00063   void show_self (  ) const;
00065   int saveFactorization (const char * file  ) const;
00069   int restoreFactorization (const char * file  , bool factor=false) ;
00071   void sort (  ) const;
00073     CoinFactorization & operator = ( const CoinFactorization & other );
00075 
00085   int factorize ( const CoinPackedMatrix & matrix, 
00086                   int rowIsBasic[], int columnIsBasic[] , 
00087                   double areaFactor = 0.0 );
00098   int factorize ( int numberRows,
00099                   int numberColumns,
00100                   CoinBigIndex numberElements,
00101                   CoinBigIndex maximumL,
00102                   CoinBigIndex maximumU,
00103                   const int indicesRow[],
00104                   const int indicesColumn[], const double elements[] ,
00105                   int permutation[],
00106                   double areaFactor = 0.0);
00111   int factorizePart1 ( int numberRows,
00112                        int numberColumns,
00113                        CoinBigIndex estimateNumberElements,
00114                        int * indicesRow[],
00115                        int * indicesColumn[],
00116                        double * elements[],
00117                   double areaFactor = 0.0);
00124   int factorizePart2 (int permutation[],int exactNumberElements);
00126   double conditionNumber() const;
00127   
00129 
00132 
00133   inline int status (  ) const {
00134     return status_;
00135   }
00137   inline void setStatus (  int value)
00138   {  status_=value;  }
00140   inline int pivots (  ) const {
00141     return numberPivots_;
00142   }
00144   inline void setPivots (  int value ) 
00145   { numberPivots_=value; }
00147   inline int *permute (  ) const {
00148     return permute_.array();
00149   }
00151   inline int *pivotColumn (  ) const {
00152     return pivotColumn_.array();
00153   }
00155   inline double *pivotRegion (  ) const {
00156     return pivotRegion_.array();
00157   }
00159   inline int *permuteBack (  ) const {
00160     return permuteBack_.array();
00161   }
00163   inline int *pivotColumnBack (  ) const {
00164     return pivotColumnBack_.array();
00165   }
00167   inline CoinBigIndex * startRowL() const
00168   { return startRowL_.array();}
00169 
00171   inline CoinBigIndex * startColumnL() const
00172   { return startColumnL_.array();}
00173 
00175   inline int * indexColumnL() const
00176   { return indexColumnL_.array();}
00177 
00179   inline int * indexRowL() const
00180   { return indexRowL_.array();}
00181 
00183   inline double * elementByRowL() const
00184   { return elementByRowL_.array();}
00185 
00187   inline int numberRowsExtra (  ) const {
00188     return numberRowsExtra_;
00189   }
00191   inline void setNumberRows(int value)
00192   { numberRows_ = value; }
00194   inline int numberRows (  ) const {
00195     return numberRows_;
00196   }
00198   inline CoinBigIndex numberL() const
00199   { return numberL_;}
00200 
00202   inline CoinBigIndex baseL() const
00203   { return baseL_;}
00205   inline int maximumRowsExtra (  ) const {
00206     return maximumRowsExtra_;
00207   }
00209   inline int numberColumns (  ) const {
00210     return numberColumns_;
00211   }
00213   inline int numberElements (  ) const {
00214     return totalElements_;
00215   }
00217   inline int numberForrestTomlin (  ) const {
00218     return numberInColumn_.array()[numberColumnsExtra_];
00219   }
00221   inline int numberGoodColumns (  ) const {
00222     return numberGoodU_;
00223   }
00225   inline double areaFactor (  ) const {
00226     return areaFactor_;
00227   }
00228   inline void areaFactor ( double value ) {
00229     areaFactor_=value;
00230   }
00232   double adjustedAreaFactor() const;
00234   inline void relaxAccuracyCheck(double value)
00235   { relaxCheck_ = value;}
00236   inline double getAccuracyCheck() const
00237   { return relaxCheck_;}
00239   inline int messageLevel (  ) const {
00240     return messageLevel_ ;
00241   }
00242   void messageLevel (  int value );
00244   inline int maximumPivots (  ) const {
00245     return maximumPivots_ ;
00246   }
00247   void maximumPivots (  int value );
00248 
00250   inline int denseThreshold() const 
00251     { return denseThreshold_;}
00253   inline void setDenseThreshold(int value)
00254     { denseThreshold_ = value;}
00256   inline double pivotTolerance (  ) const {
00257     return pivotTolerance_ ;
00258   }
00259   void pivotTolerance (  double value );
00261   inline double zeroTolerance (  ) const {
00262     return zeroTolerance_ ;
00263   }
00264   void zeroTolerance (  double value );
00265 #ifndef COIN_FAST_CODE
00267   inline double slackValue (  ) const {
00268     return slackValue_ ;
00269   }
00270   void slackValue (  double value );
00271 #endif
00273   double maximumCoefficient() const;
00275   inline bool forrestTomlin() const
00276   { return doForrestTomlin_;}
00277   inline void setForrestTomlin(bool value)
00278   { doForrestTomlin_=value;}
00280   inline bool spaceForForrestTomlin() const
00281   {
00282     CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
00283     CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
00284     return (space>=0)&&doForrestTomlin_;
00285   }
00287 
00290 
00292   inline int numberDense() const
00293   { return numberDense_;}
00294 
00296   inline CoinBigIndex numberElementsU (  ) const {
00297     return lengthU_;
00298   }
00300   inline void setNumberElementsU(CoinBigIndex value)
00301   { lengthU_ = value; }
00303   inline CoinBigIndex lengthAreaU (  ) const {
00304     return lengthAreaU_;
00305   }
00307   inline CoinBigIndex numberElementsL (  ) const {
00308     return lengthL_;
00309   }
00311   inline CoinBigIndex lengthAreaL (  ) const {
00312     return lengthAreaL_;
00313   }
00315   inline CoinBigIndex numberElementsR (  ) const {
00316     return lengthR_;
00317   }
00319   inline CoinBigIndex numberCompressions() const
00320   { return numberCompressions_;}
00322   inline int * numberInRow() const
00323   { return numberInRow_.array();}
00325   inline int * numberInColumn() const
00326   { return numberInColumn_.array();}
00328   inline double * elementU() const
00329   { return elementU_.array();}
00331   inline int * indexRowU() const
00332   { return indexRowU_.array();}
00334   inline CoinBigIndex * startColumnU() const
00335   { return startColumnU_.array();}
00337   inline int maximumColumnsExtra()
00338   { return maximumColumnsExtra_;}
00342   inline int biasLU() const
00343   { return biasLU_;}
00344   inline void setBiasLU(int value)
00345   { biasLU_=value;}
00351   inline int persistenceFlag() const
00352   { return persistenceFlag_;}
00353   void setPersistenceFlag(int value);
00355 
00358 
00366   int replaceColumn ( CoinIndexedVector * regionSparse,
00367                       int pivotRow,
00368                       double pivotCheck ,
00369                       bool checkBeforeModifying=false);
00371 
00381   int updateColumnFT ( CoinIndexedVector * regionSparse,
00382                        CoinIndexedVector * regionSparse2);
00385   int updateColumn ( CoinIndexedVector * regionSparse,
00386                      CoinIndexedVector * regionSparse2,
00387                      bool noPermute=false) const;
00393   int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
00394                            CoinIndexedVector * regionSparse2,
00395                            CoinIndexedVector * regionSparse3,
00396                            bool noPermuteRegion3=false) ;
00401   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
00402                               CoinIndexedVector * regionSparse2) const;
00404   void goSparse();
00406   inline int sparseThreshold ( ) const
00407   { return sparseThreshold_;}
00409   void sparseThreshold ( int value );
00411 
00412 
00416 
00417   inline void clearArrays()
00418   { gutsOfDestructor();}
00420 
00423 
00426   int add ( CoinBigIndex numberElements,
00427                int indicesRow[],
00428                int indicesColumn[], double elements[] );
00429 
00432   int addColumn ( CoinBigIndex numberElements,
00433                      int indicesRow[], double elements[] );
00434 
00437   int addRow ( CoinBigIndex numberElements,
00438                   int indicesColumn[], double elements[] );
00439 
00441   int deleteColumn ( int Row );
00443   int deleteRow ( int Row );
00444 
00448   int replaceRow ( int whichRow, int numberElements,
00449                       const int indicesColumn[], const double elements[] );
00451   void emptyRows(int numberToEmpty, const int which[]);
00453 
00454 
00455   void checkSparse();
00457   inline bool collectStatistics() const
00458   { return collectStatistics_;}
00460   inline void setCollectStatistics(bool onOff) const
00461   { collectStatistics_ = onOff;}
00463   void gutsOfDestructor(int type=1);
00465   void gutsOfInitialize(int type);
00466   void gutsOfCopy(const CoinFactorization &other);
00467 
00469   void resetStatistics();
00470 
00471 
00473 
00475 
00476   void getAreas ( int numberRows,
00477                   int numberColumns,
00478                   CoinBigIndex maximumL,
00479                   CoinBigIndex maximumU );
00480 
00483   void preProcess ( int state,
00484                     int possibleDuplicates = -1 );
00486   int factor (  );
00487 protected:
00490   int factorSparse (  );
00493   int factorSparseSmall (  );
00496   int factorSparseLarge (  );
00499   int factorDense (  );
00500 
00502   bool pivotOneOtherRow ( int pivotRow,
00503                           int pivotColumn );
00505   bool pivotRowSingleton ( int pivotRow,
00506                            int pivotColumn );
00508   bool pivotColumnSingleton ( int pivotRow,
00509                               int pivotColumn );
00510 
00515   bool getColumnSpace ( int iColumn,
00516                         int extraNeeded );
00517 
00521   bool getColumnSpaceIterateR ( int iColumn, double value,
00522                                int iRow);
00528   CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
00529                                int iRow);
00533   bool getRowSpace ( int iRow, int extraNeeded );
00534 
00538   bool getRowSpaceIterate ( int iRow,
00539                             int extraNeeded );
00541   void checkConsistency (  );
00543   inline void addLink ( int index, int count ) {
00544     int *nextCount = nextCount_.array();
00545     int *firstCount = firstCount_.array();
00546     int *lastCount = lastCount_.array();
00547     int next = firstCount[count];
00548       lastCount[index] = -2 - count;
00549     if ( next < 0 ) {
00550       //first with that count
00551       firstCount[count] = index;
00552       nextCount[index] = -1;
00553     } else {
00554       firstCount[count] = index;
00555       nextCount[index] = next;
00556       lastCount[next] = index;
00557   }}
00559   inline void deleteLink ( int index ) {
00560     int *nextCount = nextCount_.array();
00561     int *firstCount = firstCount_.array();
00562     int *lastCount = lastCount_.array();
00563     int next = nextCount[index];
00564     int last = lastCount[index];
00565     if ( last >= 0 ) {
00566       nextCount[last] = next;
00567     } else {
00568       int count = -last - 2;
00569 
00570       firstCount[count] = next;
00571     }
00572     if ( next >= 0 ) {
00573       lastCount[next] = last;
00574     }
00575     nextCount[index] = -2;
00576     lastCount[index] = -2;
00577     return;
00578   }
00580   void separateLinks(int count,bool rowsFirst);
00582   void cleanup (  );
00583 
00585   void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
00587   void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
00589   void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
00591   void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
00592 
00594   void updateColumnR ( CoinIndexedVector * region ) const;
00597   void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
00598 
00600   void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
00601 
00603   void updateColumnUSparse ( CoinIndexedVector * regionSparse, 
00604                              int * indexIn) const;
00606   void updateColumnUSparsish ( CoinIndexedVector * regionSparse, 
00607                                int * indexIn) const;
00609   int updateColumnUDensish ( double * COIN_RESTRICT region, 
00610                              int * COIN_RESTRICT regionIndex) const;
00612   void updateTwoColumnsUDensish (
00613                                  int & numberNonZero1,
00614                                  double * COIN_RESTRICT region1, 
00615                                  int * COIN_RESTRICT index1,
00616                                  int & numberNonZero2,
00617                                  double * COIN_RESTRICT region2, 
00618                                  int * COIN_RESTRICT index2) const;
00620   void updateColumnPFI ( CoinIndexedVector * regionSparse) const; 
00622   void permuteBack ( CoinIndexedVector * regionSparse, 
00623                      CoinIndexedVector * outVector) const;
00624 
00626   void updateColumnTransposePFI ( CoinIndexedVector * region) const;
00629   void updateColumnTransposeU ( CoinIndexedVector * region,
00630                                 int smallestIndex) const;
00633   void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
00634                                         int smallestIndex) const;
00637   void updateColumnTransposeUDensish ( CoinIndexedVector * region,
00638                                        int smallestIndex) const;
00641   void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
00642 
00644   void updateColumnTransposeR ( CoinIndexedVector * region ) const;
00646   void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
00648   void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
00649 
00651   void updateColumnTransposeL ( CoinIndexedVector * region ) const;
00653   void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
00655   void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
00657   void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
00659   void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
00660 public:
00665   int replaceColumnPFI ( CoinIndexedVector * regionSparse,
00666                          int pivotRow, double alpha);
00667 protected:
00670   int checkPivot(double saveFromU, double oldPivot) const;
00671   /********************************* START LARGE TEMPLATE ********/
00672 #ifdef INT_IS_8
00673 #define COINFACTORIZATION_BITS_PER_INT 64
00674 #define COINFACTORIZATION_SHIFT_PER_INT 6
00675 #define COINFACTORIZATION_MASK_PER_INT 0x3f
00676 #else
00677 #define COINFACTORIZATION_BITS_PER_INT 32
00678 #define COINFACTORIZATION_SHIFT_PER_INT 5
00679 #define COINFACTORIZATION_MASK_PER_INT 0x1f
00680 #endif
00681   template <class T>  inline bool
00682   pivot ( int pivotRow,
00683           int pivotColumn,
00684           CoinBigIndex pivotRowPosition,
00685           CoinBigIndex pivotColumnPosition,
00686           double work[],
00687           unsigned int workArea2[],
00688           int increment,
00689           int increment2,
00690           T markRow[] ,
00691           int largeInteger)
00692 {
00693   int *indexColumnU = indexColumnU_.array();
00694   CoinBigIndex *startColumnU = startColumnU_.array();
00695   int *numberInColumn = numberInColumn_.array();
00696   double *elementU = elementU_.array();
00697   int *indexRowU = indexRowU_.array();
00698   CoinBigIndex *startRowU = startRowU_.array();
00699   int *numberInRow = numberInRow_.array();
00700   double *elementL = elementL_.array();
00701   int *indexRowL = indexRowL_.array();
00702   int *saveColumn = saveColumn_.array();
00703   int *nextRow = nextRow_.array();
00704   int *lastRow = lastRow_.array() ;
00705 
00706   //store pivot columns (so can easily compress)
00707   int numberInPivotRow = numberInRow[pivotRow] - 1;
00708   CoinBigIndex startColumn = startColumnU[pivotColumn];
00709   int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
00710   CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
00711   int put = 0;
00712   CoinBigIndex startRow = startRowU[pivotRow];
00713   CoinBigIndex endRow = startRow + numberInPivotRow + 1;
00714 
00715   if ( pivotColumnPosition < 0 ) {
00716     for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00717       int iColumn = indexColumnU[pivotColumnPosition];
00718       if ( iColumn != pivotColumn ) {
00719         saveColumn[put++] = iColumn;
00720       } else {
00721         break;
00722       }
00723     }
00724   } else {
00725     for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
00726       saveColumn[put++] = indexColumnU[i];
00727     }
00728   }
00729   assert (pivotColumnPosition<endRow);
00730   assert (indexColumnU[pivotColumnPosition]==pivotColumn);
00731   pivotColumnPosition++;
00732   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
00733     saveColumn[put++] = indexColumnU[pivotColumnPosition];
00734   }
00735   //take out this bit of indexColumnU
00736   int next = nextRow[pivotRow];
00737   int last = lastRow[pivotRow];
00738 
00739   nextRow[last] = next;
00740   lastRow[next] = last;
00741   nextRow[pivotRow] = numberGoodU_;     //use for permute
00742   lastRow[pivotRow] = -2;
00743   numberInRow[pivotRow] = 0;
00744   //store column in L, compress in U and take column out
00745   CoinBigIndex l = lengthL_;
00746 
00747   if ( l + numberInPivotColumn > lengthAreaL_ ) {
00748     //need more memory
00749     printf("more memory needed in middle of invert\n");
00750     return false;
00751   }
00752   //l+=currentAreaL_->elementByColumn-elementL;
00753   CoinBigIndex lSave = l;
00754 
00755   pivotRowL_.array()[numberGoodL_] = pivotRow;
00756   CoinBigIndex * startColumnL = startColumnL_.array();
00757   startColumnL[numberGoodL_] = l;       //for luck and first time
00758   numberGoodL_++;
00759   startColumnL[numberGoodL_] = l + numberInPivotColumn;
00760   lengthL_ += numberInPivotColumn;
00761   if ( pivotRowPosition < 0 ) {
00762     for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00763       int iRow = indexRowU[pivotRowPosition];
00764       if ( iRow != pivotRow ) {
00765         indexRowL[l] = iRow;
00766         elementL[l] = elementU[pivotRowPosition];
00767         markRow[iRow] = static_cast<T>(l - lSave);
00768         l++;
00769         //take out of row list
00770         CoinBigIndex start = startRowU[iRow];
00771         CoinBigIndex end = start + numberInRow[iRow];
00772         CoinBigIndex where = start;
00773 
00774         while ( indexColumnU[where] != pivotColumn ) {
00775           where++;
00776         }                       /* endwhile */
00777 #if DEBUG_COIN
00778         if ( where >= end ) {
00779           abort (  );
00780         }
00781 #endif
00782         indexColumnU[where] = indexColumnU[end - 1];
00783         numberInRow[iRow]--;
00784       } else {
00785         break;
00786       }
00787     }
00788   } else {
00789     CoinBigIndex i;
00790 
00791     for ( i = startColumn; i < pivotRowPosition; i++ ) {
00792       int iRow = indexRowU[i];
00793 
00794       markRow[iRow] = static_cast<T>(l - lSave);
00795       indexRowL[l] = iRow;
00796       elementL[l] = elementU[i];
00797       l++;
00798       //take out of row list
00799       CoinBigIndex start = startRowU[iRow];
00800       CoinBigIndex end = start + numberInRow[iRow];
00801       CoinBigIndex where = start;
00802 
00803       while ( indexColumnU[where] != pivotColumn ) {
00804         where++;
00805       }                         /* endwhile */
00806 #if DEBUG_COIN
00807       if ( where >= end ) {
00808         abort (  );
00809       }
00810 #endif
00811       indexColumnU[where] = indexColumnU[end - 1];
00812       numberInRow[iRow]--;
00813       assert (numberInRow[iRow]>=0);
00814     }
00815   }
00816   assert (pivotRowPosition<endColumn);
00817   assert (indexRowU[pivotRowPosition]==pivotRow);
00818   double pivotElement = elementU[pivotRowPosition];
00819   double pivotMultiplier = 1.0 / pivotElement;
00820 
00821   pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
00822   pivotRowPosition++;
00823   for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
00824     int iRow = indexRowU[pivotRowPosition];
00825     
00826     markRow[iRow] = static_cast<T>(l - lSave);
00827     indexRowL[l] = iRow;
00828     elementL[l] = elementU[pivotRowPosition];
00829     l++;
00830     //take out of row list
00831     CoinBigIndex start = startRowU[iRow];
00832     CoinBigIndex end = start + numberInRow[iRow];
00833     CoinBigIndex where = start;
00834     
00835     while ( indexColumnU[where] != pivotColumn ) {
00836       where++;
00837     }                           /* endwhile */
00838 #if DEBUG_COIN
00839     if ( where >= end ) {
00840       abort (  );
00841     }
00842 #endif
00843     indexColumnU[where] = indexColumnU[end - 1];
00844     numberInRow[iRow]--;
00845     assert (numberInRow[iRow]>=0);
00846   }
00847   markRow[pivotRow] = static_cast<T>(largeInteger);
00848   //compress pivot column (move pivot to front including saved)
00849   numberInColumn[pivotColumn] = 0;
00850   //use end of L for temporary space
00851   int *indexL = &indexRowL[lSave];
00852   double *multipliersL = &elementL[lSave];
00853 
00854   //adjust
00855   int j;
00856 
00857   for ( j = 0; j < numberInPivotColumn; j++ ) {
00858     multipliersL[j] *= pivotMultiplier;
00859   }
00860   //zero out fill
00861   CoinBigIndex iErase;
00862   for ( iErase = 0; iErase < increment2 * numberInPivotRow;
00863         iErase++ ) {
00864     workArea2[iErase] = 0;
00865   }
00866   CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
00867   unsigned int *temp2 = workArea2;
00868   int * nextColumn = nextColumn_.array();
00869 
00870   //pack down and move to work
00871   int jColumn;
00872   for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
00873     int iColumn = saveColumn[jColumn];
00874     CoinBigIndex startColumn = startColumnU[iColumn];
00875     CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
00876     int iRow = indexRowU[startColumn];
00877     double value = elementU[startColumn];
00878     double largest;
00879     CoinBigIndex put = startColumn;
00880     CoinBigIndex positionLargest = -1;
00881     double thisPivotValue = 0.0;
00882 
00883     //compress column and find largest not updated
00884     bool checkLargest;
00885     int mark = markRow[iRow];
00886 
00887     if ( mark == largeInteger+1 ) {
00888       largest = fabs ( value );
00889       positionLargest = put;
00890       put++;
00891       checkLargest = false;
00892     } else {
00893       //need to find largest
00894       largest = 0.0;
00895       checkLargest = true;
00896       if ( mark != largeInteger ) {
00897         //will be updated
00898         work[mark] = value;
00899         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00900         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00901 
00902         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
00903         added--;
00904       } else {
00905         thisPivotValue = value;
00906       }
00907     }
00908     CoinBigIndex i;
00909     for ( i = startColumn + 1; i < endColumn; i++ ) {
00910       iRow = indexRowU[i];
00911       value = elementU[i];
00912       int mark = markRow[iRow];
00913 
00914       if ( mark == largeInteger+1 ) {
00915         //keep
00916         indexRowU[put] = iRow;
00917         elementU[put] = value;
00918         if ( checkLargest ) {
00919           double absValue = fabs ( value );
00920 
00921           if ( absValue > largest ) {
00922             largest = absValue;
00923             positionLargest = put;
00924           }
00925         }
00926         put++;
00927       } else if ( mark != largeInteger ) {
00928         //will be updated
00929         work[mark] = value;
00930         int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
00931         int bit = mark & COINFACTORIZATION_MASK_PER_INT;
00932 
00933         temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
00934         added--;
00935       } else {
00936         thisPivotValue = value;
00937       }
00938     }
00939     //slot in pivot
00940     elementU[put] = elementU[startColumn];
00941     indexRowU[put] = indexRowU[startColumn];
00942     if ( positionLargest == startColumn ) {
00943       positionLargest = put;    //follow if was largest
00944     }
00945     put++;
00946     elementU[startColumn] = thisPivotValue;
00947     indexRowU[startColumn] = pivotRow;
00948     //clean up counts
00949     startColumn++;
00950     numberInColumn[iColumn] = put - startColumn;
00951     int * numberInColumnPlus = numberInColumnPlus_.array();
00952     numberInColumnPlus[iColumn]++;
00953     startColumnU[iColumn]++;
00954     //how much space have we got
00955     int next = nextColumn[iColumn];
00956     CoinBigIndex space;
00957 
00958     space = startColumnU[next] - put - numberInColumnPlus[next];
00959     //assume no zero elements
00960     if ( numberInPivotColumn > space ) {
00961       //getColumnSpace also moves fixed part
00962       if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
00963         return false;
00964       }
00965       //redo starts
00966       positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
00967       startColumn = startColumnU[iColumn];
00968       put = startColumn + numberInColumn[iColumn];
00969     }
00970     double tolerance = zeroTolerance_;
00971 
00972     int *nextCount = nextCount_.array();
00973     for ( j = 0; j < numberInPivotColumn; j++ ) {
00974       value = work[j] - thisPivotValue * multipliersL[j];
00975       double absValue = fabs ( value );
00976 
00977       if ( absValue > tolerance ) {
00978         work[j] = 0.0;
00979         elementU[put] = value;
00980         indexRowU[put] = indexL[j];
00981         if ( absValue > largest ) {
00982           largest = absValue;
00983           positionLargest = put;
00984         }
00985         put++;
00986       } else {
00987         work[j] = 0.0;
00988         added--;
00989         int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
00990         int bit = j & COINFACTORIZATION_MASK_PER_INT;
00991 
00992         if ( temp2[word] & ( 1 << bit ) ) {
00993           //take out of row list
00994           iRow = indexL[j];
00995           CoinBigIndex start = startRowU[iRow];
00996           CoinBigIndex end = start + numberInRow[iRow];
00997           CoinBigIndex where = start;
00998 
00999           while ( indexColumnU[where] != iColumn ) {
01000             where++;
01001           }                     /* endwhile */
01002 #if DEBUG_COIN
01003           if ( where >= end ) {
01004             abort (  );
01005           }
01006 #endif
01007           indexColumnU[where] = indexColumnU[end - 1];
01008           numberInRow[iRow]--;
01009         } else {
01010           //make sure won't be added
01011           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01012           int bit = j & COINFACTORIZATION_MASK_PER_INT;
01013 
01014           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
01015         }
01016       }
01017     }
01018     numberInColumn[iColumn] = put - startColumn;
01019     //move largest
01020     if ( positionLargest >= 0 ) {
01021       value = elementU[positionLargest];
01022       iRow = indexRowU[positionLargest];
01023       elementU[positionLargest] = elementU[startColumn];
01024       indexRowU[positionLargest] = indexRowU[startColumn];
01025       elementU[startColumn] = value;
01026       indexRowU[startColumn] = iRow;
01027     }
01028     //linked list for column
01029     if ( nextCount[iColumn + numberRows_] != -2 ) {
01030       //modify linked list
01031       deleteLink ( iColumn + numberRows_ );
01032       addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01033     }
01034     temp2 += increment2;
01035   }
01036   //get space for row list
01037   unsigned int *putBase = workArea2;
01038   int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01039   int i = 0;
01040 
01041   // do linked lists and update counts
01042   while ( bigLoops ) {
01043     bigLoops--;
01044     int bit;
01045     for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01046       unsigned int *putThis = putBase;
01047       int iRow = indexL[i];
01048 
01049       //get space
01050       int number = 0;
01051       int jColumn;
01052 
01053       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01054         unsigned int test = *putThis;
01055 
01056         putThis += increment2;
01057         test = 1 - ( ( test >> bit ) & 1 );
01058         number += test;
01059       }
01060       int next = nextRow[iRow];
01061       CoinBigIndex space;
01062 
01063       space = startRowU[next] - startRowU[iRow];
01064       number += numberInRow[iRow];
01065       if ( space < number ) {
01066         if ( !getRowSpace ( iRow, number ) ) {
01067           return false;
01068         }
01069       }
01070       // now do
01071       putThis = putBase;
01072       next = nextRow[iRow];
01073       number = numberInRow[iRow];
01074       CoinBigIndex end = startRowU[iRow] + number;
01075       int saveIndex = indexColumnU[startRowU[next]];
01076 
01077       //add in
01078       for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01079         unsigned int test = *putThis;
01080 
01081         putThis += increment2;
01082         test = 1 - ( ( test >> bit ) & 1 );
01083         indexColumnU[end] = saveColumn[jColumn];
01084         end += test;
01085       }
01086       //put back next one in case zapped
01087       indexColumnU[startRowU[next]] = saveIndex;
01088       markRow[iRow] = static_cast<T>(largeInteger+1);
01089       number = end - startRowU[iRow];
01090       numberInRow[iRow] = number;
01091       deleteLink ( iRow );
01092       addLink ( iRow, number );
01093     }
01094     putBase++;
01095   }                             /* endwhile */
01096   int bit;
01097 
01098   for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
01099     unsigned int *putThis = putBase;
01100     int iRow = indexL[i];
01101 
01102     //get space
01103     int number = 0;
01104     int jColumn;
01105 
01106     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01107       unsigned int test = *putThis;
01108 
01109       putThis += increment2;
01110       test = 1 - ( ( test >> bit ) & 1 );
01111       number += test;
01112     }
01113     int next = nextRow[iRow];
01114     CoinBigIndex space;
01115 
01116     space = startRowU[next] - startRowU[iRow];
01117     number += numberInRow[iRow];
01118     if ( space < number ) {
01119       if ( !getRowSpace ( iRow, number ) ) {
01120         return false;
01121       }
01122     }
01123     // now do
01124     putThis = putBase;
01125     next = nextRow[iRow];
01126     number = numberInRow[iRow];
01127     CoinBigIndex end = startRowU[iRow] + number;
01128     int saveIndex;
01129 
01130     saveIndex = indexColumnU[startRowU[next]];
01131 
01132     //add in
01133     for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
01134       unsigned int test = *putThis;
01135 
01136       putThis += increment2;
01137       test = 1 - ( ( test >> bit ) & 1 );
01138 
01139       indexColumnU[end] = saveColumn[jColumn];
01140       end += test;
01141     }
01142     indexColumnU[startRowU[next]] = saveIndex;
01143     markRow[iRow] = static_cast<T>(largeInteger+1);
01144     number = end - startRowU[iRow];
01145     numberInRow[iRow] = number;
01146     deleteLink ( iRow );
01147     addLink ( iRow, number );
01148   }
01149   markRow[pivotRow] = static_cast<T>(largeInteger+1);
01150   //modify linked list for pivots
01151   deleteLink ( pivotRow );
01152   deleteLink ( pivotColumn + numberRows_ );
01153   totalElements_ += added;
01154   return true;
01155 }
01156 
01157   /********************************* END LARGE TEMPLATE ********/
01159 
01160 protected:
01161 
01164 
01165   double pivotTolerance_;
01167   double zeroTolerance_;
01168 #ifndef COIN_FAST_CODE
01170   double slackValue_;
01171 #else
01172 #ifndef slackValue_
01173 #define slackValue_ -1.0
01174 #endif
01175 #endif
01177   double areaFactor_;
01179   double relaxCheck_;
01181   int numberRows_;
01183   int numberRowsExtra_;
01185   int maximumRowsExtra_;
01187   int numberColumns_;
01189   int numberColumnsExtra_;
01191   int maximumColumnsExtra_;
01193   int numberGoodU_;
01195   int numberGoodL_;
01197   int maximumPivots_;
01199   int numberPivots_;
01202   CoinBigIndex totalElements_;
01204   CoinBigIndex factorElements_;
01206   CoinIntArrayWithLength pivotColumn_;
01208   CoinIntArrayWithLength permute_;
01210   CoinIntArrayWithLength permuteBack_;
01212   CoinIntArrayWithLength pivotColumnBack_;
01214   int status_;
01215 
01220   //int increasingRows_;
01221 
01223   int numberTrials_;
01225   CoinBigIndexArrayWithLength startRowU_;
01226 
01228   CoinIntArrayWithLength numberInRow_;
01229 
01231   CoinIntArrayWithLength numberInColumn_;
01232 
01234   CoinIntArrayWithLength numberInColumnPlus_;
01235 
01238   CoinIntArrayWithLength firstCount_;
01239 
01241   CoinIntArrayWithLength nextCount_;
01242 
01244   CoinIntArrayWithLength lastCount_;
01245 
01247   CoinIntArrayWithLength nextColumn_;
01248 
01250   CoinIntArrayWithLength lastColumn_;
01251 
01253   CoinIntArrayWithLength nextRow_;
01254 
01256   CoinIntArrayWithLength lastRow_;
01257 
01259   CoinIntArrayWithLength saveColumn_;
01260 
01262   CoinIntArrayWithLength markRow_;
01263 
01265   int messageLevel_;
01266 
01268   int biggerDimension_;
01269 
01271   CoinIntArrayWithLength indexColumnU_;
01272 
01274   CoinIntArrayWithLength pivotRowL_;
01275 
01277   CoinDoubleArrayWithLength pivotRegion_;
01278 
01280   int numberSlacks_;
01281 
01283   int numberU_;
01284 
01286   CoinBigIndex maximumU_;
01287 
01289   //int baseU_;
01290 
01292   CoinBigIndex lengthU_;
01293 
01295   CoinBigIndex lengthAreaU_;
01296 
01298   CoinDoubleArrayWithLength elementU_;
01299 
01301   CoinIntArrayWithLength indexRowU_;
01302 
01304   CoinBigIndexArrayWithLength startColumnU_;
01305 
01307   CoinBigIndexArrayWithLength convertRowToColumnU_;
01308 
01310   CoinBigIndex numberL_;
01311 
01313   CoinBigIndex baseL_;
01314 
01316   CoinBigIndex lengthL_;
01317 
01319   CoinBigIndex lengthAreaL_;
01320 
01322   CoinDoubleArrayWithLength elementL_;
01323 
01325   CoinIntArrayWithLength indexRowL_;
01326 
01328   CoinBigIndexArrayWithLength startColumnL_;
01329 
01331   bool doForrestTomlin_;
01332 
01334   int numberR_;
01335 
01337   CoinBigIndex lengthR_;
01338 
01340   CoinBigIndex lengthAreaR_;
01341 
01343   double *elementR_;
01344 
01346   int *indexRowR_;
01347 
01349   CoinBigIndexArrayWithLength startColumnR_;
01350 
01352   double  * denseArea_;
01353 
01355   int * densePermute_;
01356 
01358   int numberDense_;
01359 
01361   int denseThreshold_;
01362 
01364   CoinDoubleArrayWithLength workArea_;
01365 
01367   CoinUnsignedIntArrayWithLength workArea2_;
01368 
01370   CoinBigIndex numberCompressions_;
01371 
01373   mutable double ftranCountInput_;
01374   mutable double ftranCountAfterL_;
01375   mutable double ftranCountAfterR_;
01376   mutable double ftranCountAfterU_;
01377   mutable double btranCountInput_;
01378   mutable double btranCountAfterU_;
01379   mutable double btranCountAfterR_;
01380   mutable double btranCountAfterL_;
01381 
01383   mutable int numberFtranCounts_;
01384   mutable int numberBtranCounts_;
01385 
01387   double ftranAverageAfterL_;
01388   double ftranAverageAfterR_;
01389   double ftranAverageAfterU_;
01390   double btranAverageAfterU_;
01391   double btranAverageAfterR_;
01392   double btranAverageAfterL_;
01393 
01395   mutable bool collectStatistics_;
01396 
01398   int sparseThreshold_;
01399 
01401   int sparseThreshold2_;
01402 
01404   CoinBigIndexArrayWithLength startRowL_;
01405 
01407   CoinIntArrayWithLength indexColumnL_;
01408 
01410   CoinDoubleArrayWithLength elementByRowL_;
01411 
01413   mutable CoinIntArrayWithLength sparse_;
01417   int biasLU_;
01423   int persistenceFlag_;
01425 };
01426 // Dense coding
01427 #ifdef COIN_HAS_LAPACK
01428 #define DENSE_CODE 1
01429 /* Type of Fortran integer translated into C */
01430 #ifndef ipfint
01431 //typedef ipfint FORTRAN_INTEGER_TYPE ;
01432 typedef int ipfint;
01433 typedef const int cipfint;
01434 #endif
01435 #endif
01436 #endif
01437 // Extra for ugly include
01438 #ifdef UGLY_COIN_FACTOR_CODING
01439 #define FAC_UNSET (FAC_SET+1)
01440 {
01441   goodPivot=false;
01442   //store pivot columns (so can easily compress)
01443   CoinBigIndex startColumnThis = startColumn[iPivotColumn];
01444   CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
01445   int put = 0;
01446   CoinBigIndex startRowThis = startRow[iPivotRow];
01447   CoinBigIndex endRow = startRowThis + numberDoRow + 1;
01448   if ( pivotColumnPosition < 0 ) {
01449     for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01450       int iColumn = indexColumn[pivotColumnPosition];
01451       if ( iColumn != iPivotColumn ) {
01452         saveColumn[put++] = iColumn;
01453       } else {
01454         break;
01455       }
01456     }
01457   } else {
01458     for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
01459       saveColumn[put++] = indexColumn[i];
01460     }
01461   }
01462   assert (pivotColumnPosition<endRow);
01463   assert (indexColumn[pivotColumnPosition]==iPivotColumn);
01464   pivotColumnPosition++;
01465   for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
01466     saveColumn[put++] = indexColumn[pivotColumnPosition];
01467   }
01468   //take out this bit of indexColumn
01469   int next = nextRow[iPivotRow];
01470   int last = lastRow[iPivotRow];
01471   
01472   nextRow[last] = next;
01473   lastRow[next] = last;
01474   nextRow[iPivotRow] = numberGoodU_;    //use for permute
01475   lastRow[iPivotRow] = -2;
01476   numberInRow[iPivotRow] = 0;
01477   //store column in L, compress in U and take column out
01478   CoinBigIndex l = lengthL_;
01479   // **** HORRID coding coming up but a goto seems best!
01480   {
01481     if ( l + numberDoColumn > lengthAreaL_ ) {
01482       //need more memory
01483       printf("more memory needed in middle of invert\n");
01484       goto BAD_PIVOT;
01485     }
01486     //l+=currentAreaL_->elementByColumn-elementL;
01487     CoinBigIndex lSave = l;
01488     
01489     pivotRowL_.array()[numberGoodL_] = iPivotRow;
01490     CoinBigIndex * startColumnL = startColumnL_.array();
01491     startColumnL[numberGoodL_] = l;     //for luck and first time
01492     numberGoodL_++;
01493     startColumnL[numberGoodL_] = l + numberDoColumn;
01494     lengthL_ += numberDoColumn;
01495     if ( pivotRowPosition < 0 ) {
01496       for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01497         int iRow = indexRow[pivotRowPosition];
01498         if ( iRow != iPivotRow ) {
01499           indexRowL[l] = iRow;
01500           elementL[l] = element[pivotRowPosition];
01501           markRow[iRow] = l - lSave;
01502           l++;
01503           //take out of row list
01504           CoinBigIndex start = startRow[iRow];
01505           CoinBigIndex end = start + numberInRow[iRow];
01506           CoinBigIndex where = start;
01507           
01508           while ( indexColumn[where] != iPivotColumn ) {
01509             where++;
01510           }                     /* endwhile */
01511 #if DEBUG_COIN
01512           if ( where >= end ) {
01513             abort (  );
01514           }
01515 #endif
01516           indexColumn[where] = indexColumn[end - 1];
01517           numberInRow[iRow]--;
01518         } else {
01519           break;
01520         }
01521       }
01522     } else {
01523       CoinBigIndex i;
01524       
01525       for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
01526         int iRow = indexRow[i];
01527         
01528         markRow[iRow] = l - lSave;
01529         indexRowL[l] = iRow;
01530         elementL[l] = element[i];
01531         l++;
01532         //take out of row list
01533         CoinBigIndex start = startRow[iRow];
01534         CoinBigIndex end = start + numberInRow[iRow];
01535         CoinBigIndex where = start;
01536         
01537         while ( indexColumn[where] != iPivotColumn ) {
01538           where++;
01539         }                               /* endwhile */
01540 #if DEBUG_COIN
01541         if ( where >= end ) {
01542           abort (  );
01543         }
01544 #endif
01545         indexColumn[where] = indexColumn[end - 1];
01546         numberInRow[iRow]--;
01547         assert (numberInRow[iRow]>=0);
01548       }
01549     }
01550     assert (pivotRowPosition<endColumn);
01551     assert (indexRow[pivotRowPosition]==iPivotRow);
01552     double pivotElement = element[pivotRowPosition];
01553     double pivotMultiplier = 1.0 / pivotElement;
01554     
01555     pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
01556     pivotRowPosition++;
01557     for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
01558       int iRow = indexRow[pivotRowPosition];
01559       
01560       markRow[iRow] = l - lSave;
01561       indexRowL[l] = iRow;
01562       elementL[l] = element[pivotRowPosition];
01563       l++;
01564       //take out of row list
01565       CoinBigIndex start = startRow[iRow];
01566       CoinBigIndex end = start + numberInRow[iRow];
01567       CoinBigIndex where = start;
01568       
01569       while ( indexColumn[where] != iPivotColumn ) {
01570         where++;
01571       }                         /* endwhile */
01572 #if DEBUG_COIN
01573       if ( where >= end ) {
01574         abort (  );
01575       }
01576 #endif
01577       indexColumn[where] = indexColumn[end - 1];
01578       numberInRow[iRow]--;
01579       assert (numberInRow[iRow]>=0);
01580     }
01581     markRow[iPivotRow] = FAC_SET;
01582     //compress pivot column (move pivot to front including saved)
01583     numberInColumn[iPivotColumn] = 0;
01584     //use end of L for temporary space
01585     int *indexL = &indexRowL[lSave];
01586     double *multipliersL = &elementL[lSave];
01587     
01588     //adjust
01589     int j;
01590     
01591     for ( j = 0; j < numberDoColumn; j++ ) {
01592       multipliersL[j] *= pivotMultiplier;
01593     }
01594     //zero out fill
01595     CoinBigIndex iErase;
01596     for ( iErase = 0; iErase < increment2 * numberDoRow;
01597           iErase++ ) {
01598       workArea2[iErase] = 0;
01599     }
01600     CoinBigIndex added = numberDoRow * numberDoColumn;
01601     unsigned int *temp2 = workArea2;
01602     int * nextColumn = nextColumn_.array();
01603     
01604     //pack down and move to work
01605     int jColumn;
01606     for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01607       int iColumn = saveColumn[jColumn];
01608       CoinBigIndex startColumnThis = startColumn[iColumn];
01609       CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
01610       int iRow = indexRow[startColumnThis];
01611       double value = element[startColumnThis];
01612       double largest;
01613       CoinBigIndex put = startColumnThis;
01614       CoinBigIndex positionLargest = -1;
01615       double thisPivotValue = 0.0;
01616       
01617       //compress column and find largest not updated
01618       bool checkLargest;
01619       int mark = markRow[iRow];
01620       
01621       if ( mark == FAC_UNSET ) {
01622         largest = fabs ( value );
01623         positionLargest = put;
01624         put++;
01625         checkLargest = false;
01626       } else {
01627         //need to find largest
01628         largest = 0.0;
01629         checkLargest = true;
01630         if ( mark != FAC_SET ) {
01631           //will be updated
01632           workArea[mark] = value;
01633           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01634           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01635           
01636           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
01637           added--;
01638         } else {
01639           thisPivotValue = value;
01640         }
01641       }
01642       CoinBigIndex i;
01643       for ( i = startColumnThis + 1; i < endColumn; i++ ) {
01644         iRow = indexRow[i];
01645         value = element[i];
01646         int mark = markRow[iRow];
01647         
01648         if ( mark == FAC_UNSET ) {
01649           //keep
01650           indexRow[put] = iRow;
01651           element[put] = value;
01652           if ( checkLargest ) {
01653             double absValue = fabs ( value );
01654             
01655             if ( absValue > largest ) {
01656               largest = absValue;
01657               positionLargest = put;
01658             }
01659           }
01660           put++;
01661         } else if ( mark != FAC_SET ) {
01662           //will be updated
01663           workArea[mark] = value;
01664           int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
01665           int bit = mark & COINFACTORIZATION_MASK_PER_INT;
01666           
01667           temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
01668           added--;
01669         } else {
01670           thisPivotValue = value;
01671         }
01672       }
01673       //slot in pivot
01674       element[put] = element[startColumnThis];
01675       indexRow[put] = indexRow[startColumnThis];
01676       if ( positionLargest == startColumnThis ) {
01677         positionLargest = put;  //follow if was largest
01678       }
01679       put++;
01680       element[startColumnThis] = thisPivotValue;
01681       indexRow[startColumnThis] = iPivotRow;
01682       //clean up counts
01683       startColumnThis++;
01684       numberInColumn[iColumn] = put - startColumnThis;
01685       int * numberInColumnPlus = numberInColumnPlus_.array();
01686       numberInColumnPlus[iColumn]++;
01687       startColumn[iColumn]++;
01688       //how much space have we got
01689       int next = nextColumn[iColumn];
01690       CoinBigIndex space;
01691       
01692       space = startColumn[next] - put - numberInColumnPlus[next];
01693       //assume no zero elements
01694       if ( numberDoColumn > space ) {
01695         //getColumnSpace also moves fixed part
01696         if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
01697           goto BAD_PIVOT;
01698         }
01699         //redo starts
01700         positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
01701         startColumnThis = startColumn[iColumn];
01702         put = startColumnThis + numberInColumn[iColumn];
01703       }
01704       double tolerance = zeroTolerance_;
01705       
01706       int *nextCount = nextCount_.array();
01707       for ( j = 0; j < numberDoColumn; j++ ) {
01708         value = workArea[j] - thisPivotValue * multipliersL[j];
01709         double absValue = fabs ( value );
01710         
01711         if ( absValue > tolerance ) {
01712           workArea[j] = 0.0;
01713           element[put] = value;
01714           indexRow[put] = indexL[j];
01715           if ( absValue > largest ) {
01716             largest = absValue;
01717             positionLargest = put;
01718           }
01719           put++;
01720         } else {
01721           workArea[j] = 0.0;
01722           added--;
01723           int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01724           int bit = j & COINFACTORIZATION_MASK_PER_INT;
01725           
01726           if ( temp2[word] & ( 1 << bit ) ) {
01727             //take out of row list
01728             iRow = indexL[j];
01729             CoinBigIndex start = startRow[iRow];
01730             CoinBigIndex end = start + numberInRow[iRow];
01731             CoinBigIndex where = start;
01732             
01733             while ( indexColumn[where] != iColumn ) {
01734               where++;
01735             }                   /* endwhile */
01736 #if DEBUG_COIN
01737             if ( where >= end ) {
01738               abort (  );
01739             }
01740 #endif
01741             indexColumn[where] = indexColumn[end - 1];
01742             numberInRow[iRow]--;
01743           } else {
01744             //make sure won't be added
01745             int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
01746             int bit = j & COINFACTORIZATION_MASK_PER_INT;
01747             
01748             temp2[word] = temp2[word] | ( 1 << bit );   //say already in counts
01749           }
01750         }
01751       }
01752       numberInColumn[iColumn] = put - startColumnThis;
01753       //move largest
01754       if ( positionLargest >= 0 ) {
01755         value = element[positionLargest];
01756         iRow = indexRow[positionLargest];
01757         element[positionLargest] = element[startColumnThis];
01758         indexRow[positionLargest] = indexRow[startColumnThis];
01759         element[startColumnThis] = value;
01760         indexRow[startColumnThis] = iRow;
01761       }
01762       //linked list for column
01763       if ( nextCount[iColumn + numberRows_] != -2 ) {
01764         //modify linked list
01765         deleteLink ( iColumn + numberRows_ );
01766         addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
01767       }
01768       temp2 += increment2;
01769     }
01770     //get space for row list
01771     unsigned int *putBase = workArea2;
01772     int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
01773     int i = 0;
01774     
01775     // do linked lists and update counts
01776     while ( bigLoops ) {
01777       bigLoops--;
01778       int bit;
01779       for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
01780         unsigned int *putThis = putBase;
01781         int iRow = indexL[i];
01782         
01783         //get space
01784         int number = 0;
01785         int jColumn;
01786         
01787         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01788           unsigned int test = *putThis;
01789           
01790           putThis += increment2;
01791           test = 1 - ( ( test >> bit ) & 1 );
01792           number += test;
01793         }
01794         int next = nextRow[iRow];
01795         CoinBigIndex space;
01796         
01797         space = startRow[next] - startRow[iRow];
01798         number += numberInRow[iRow];
01799         if ( space < number ) {
01800           if ( !getRowSpace ( iRow, number ) ) {
01801             goto BAD_PIVOT;
01802           }
01803         }
01804         // now do
01805         putThis = putBase;
01806         next = nextRow[iRow];
01807         number = numberInRow[iRow];
01808         CoinBigIndex end = startRow[iRow] + number;
01809         int saveIndex = indexColumn[startRow[next]];
01810         
01811         //add in
01812         for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01813           unsigned int test = *putThis;
01814           
01815           putThis += increment2;
01816           test = 1 - ( ( test >> bit ) & 1 );
01817           indexColumn[end] = saveColumn[jColumn];
01818           end += test;
01819         }
01820         //put back next one in case zapped
01821         indexColumn[startRow[next]] = saveIndex;
01822         markRow[iRow] = FAC_UNSET;
01823         number = end - startRow[iRow];
01824         numberInRow[iRow] = number;
01825         deleteLink ( iRow );
01826         addLink ( iRow, number );
01827       }
01828       putBase++;
01829     }                           /* endwhile */
01830     int bit;
01831     
01832     for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
01833       unsigned int *putThis = putBase;
01834       int iRow = indexL[i];
01835       
01836       //get space
01837       int number = 0;
01838       int jColumn;
01839       
01840       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01841         unsigned int test = *putThis;
01842         
01843         putThis += increment2;
01844         test = 1 - ( ( test >> bit ) & 1 );
01845         number += test;
01846       }
01847       int next = nextRow[iRow];
01848       CoinBigIndex space;
01849       
01850       space = startRow[next] - startRow[iRow];
01851       number += numberInRow[iRow];
01852       if ( space < number ) {
01853         if ( !getRowSpace ( iRow, number ) ) {
01854           goto BAD_PIVOT;
01855         }
01856       }
01857       // now do
01858       putThis = putBase;
01859       next = nextRow[iRow];
01860       number = numberInRow[iRow];
01861       CoinBigIndex end = startRow[iRow] + number;
01862       int saveIndex;
01863       
01864       saveIndex = indexColumn[startRow[next]];
01865       
01866       //add in
01867       for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
01868         unsigned int test = *putThis;
01869         
01870         putThis += increment2;
01871         test = 1 - ( ( test >> bit ) & 1 );
01872         
01873         indexColumn[end] = saveColumn[jColumn];
01874         end += test;
01875       }
01876       indexColumn[startRow[next]] = saveIndex;
01877       markRow[iRow] = FAC_UNSET;
01878       number = end - startRow[iRow];
01879       numberInRow[iRow] = number;
01880       deleteLink ( iRow );
01881       addLink ( iRow, number );
01882     }
01883     markRow[iPivotRow] = FAC_UNSET;
01884     //modify linked list for pivots
01885     deleteLink ( iPivotRow );
01886     deleteLink ( iPivotColumn + numberRows_ );
01887     totalElements_ += added;
01888     goodPivot= true;
01889     // **** UGLY UGLY UGLY
01890   }
01891  BAD_PIVOT:
01892 
01893   ;
01894 }
01895 #undef FAC_UNSET
01896 #endif

Generated on Thu Jan 15 03:01:00 2009 for coin-Bcp by  doxygen 1.4.7