/home/coin/SVN-release/Clp-1.8.2/CoinUtils/src/CoinFactorization.hpp

Go to the documentation of this file.

// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.

/* 
   Authors
   
   John Forrest

 */
#ifndef CoinFactorization_H
#define CoinFactorization_H

#include <iostream>
#include <string>
#include <cassert>
#include <cstdio>
#include "CoinFinite.hpp"
#include "CoinIndexedVector.hpp"
class CoinPackedMatrix;
class CoinFactorization {
   friend void CoinFactorizationUnitTest( const std::string & mpsDir );

public:


    CoinFactorization (  );
  CoinFactorization ( const CoinFactorization &other);

   ~CoinFactorization (  );
  void almostDestructor();
  void show_self (  ) const;
  int saveFactorization (const char * file  ) const;
  int restoreFactorization (const char * file  , bool factor=false) ;
  void sort (  ) const;
    CoinFactorization & operator = ( const CoinFactorization & other );

  int factorize ( const CoinPackedMatrix & matrix, 
                  int rowIsBasic[], int columnIsBasic[] , 
                  double areaFactor = 0.0 );
  int factorize ( int numberRows,
                  int numberColumns,
                  CoinBigIndex numberElements,
                  CoinBigIndex maximumL,
                  CoinBigIndex maximumU,
                  const int indicesRow[],
                  const int indicesColumn[], const double elements[] ,
                  int permutation[],
                  double areaFactor = 0.0);
  int factorizePart1 ( int numberRows,
                       int numberColumns,
                       CoinBigIndex estimateNumberElements,
                       int * indicesRow[],
                       int * indicesColumn[],
                       double * elements[],
                  double areaFactor = 0.0);
  int factorizePart2 (int permutation[],int exactNumberElements);
  double conditionNumber() const;
  


  inline int status (  ) const {
    return status_;
  }
  inline void setStatus (  int value)
  {  status_=value;  }
  inline int pivots (  ) const {
    return numberPivots_;
  }
  inline void setPivots (  int value ) 
  { numberPivots_=value; }
  inline int *permute (  ) const {
    return permute_.array();
  }
  inline int *pivotColumn (  ) const {
    return pivotColumn_.array();
  }
  inline double *pivotRegion (  ) const {
    return pivotRegion_.array();
  }
  inline int *permuteBack (  ) const {
    return permuteBack_.array();
  }
  inline int *pivotColumnBack (  ) const {
    return pivotColumnBack_.array();
  }
  inline CoinBigIndex * startRowL() const
  { return startRowL_.array();}

  inline CoinBigIndex * startColumnL() const
  { return startColumnL_.array();}

  inline int * indexColumnL() const
  { return indexColumnL_.array();}

  inline int * indexRowL() const
  { return indexRowL_.array();}

  inline double * elementByRowL() const
  { return elementByRowL_.array();}

  inline int numberRowsExtra (  ) const {
    return numberRowsExtra_;
  }
  inline void setNumberRows(int value)
  { numberRows_ = value; }
  inline int numberRows (  ) const {
    return numberRows_;
  }
  inline CoinBigIndex numberL() const
  { return numberL_;}

  inline CoinBigIndex baseL() const
  { return baseL_;}
  inline int maximumRowsExtra (  ) const {
    return maximumRowsExtra_;
  }
  inline int numberColumns (  ) const {
    return numberColumns_;
  }
  inline int numberElements (  ) const {
    return totalElements_;
  }
  inline int numberForrestTomlin (  ) const {
    return numberInColumn_.array()[numberColumnsExtra_];
  }
  inline int numberGoodColumns (  ) const {
    return numberGoodU_;
  }
  inline double areaFactor (  ) const {
    return areaFactor_;
  }
  inline void areaFactor ( double value ) {
    areaFactor_=value;
  }
  double adjustedAreaFactor() const;
  inline void relaxAccuracyCheck(double value)
  { relaxCheck_ = value;}
  inline double getAccuracyCheck() const
  { return relaxCheck_;}
  inline int messageLevel (  ) const {
    return messageLevel_ ;
  }
  void messageLevel (  int value );
  inline int maximumPivots (  ) const {
    return maximumPivots_ ;
  }
  void maximumPivots (  int value );

  inline int denseThreshold() const 
    { return denseThreshold_;}
  inline void setDenseThreshold(int value)
    { denseThreshold_ = value;}
  inline double pivotTolerance (  ) const {
    return pivotTolerance_ ;
  }
  void pivotTolerance (  double value );
  inline double zeroTolerance (  ) const {
    return zeroTolerance_ ;
  }
  void zeroTolerance (  double value );
#ifndef COIN_FAST_CODE
  inline double slackValue (  ) const {
    return slackValue_ ;
  }
  void slackValue (  double value );
#endif
  double maximumCoefficient() const;
  inline bool forrestTomlin() const
  { return doForrestTomlin_;}
  inline void setForrestTomlin(bool value)
  { doForrestTomlin_=value;}
  inline bool spaceForForrestTomlin() const
  {
    CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
    CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
    return (space>=0)&&doForrestTomlin_;
  }


  inline int numberDense() const
  { return numberDense_;}

  inline CoinBigIndex numberElementsU (  ) const {
    return lengthU_;
  }
  inline void setNumberElementsU(CoinBigIndex value)
  { lengthU_ = value; }
  inline CoinBigIndex lengthAreaU (  ) const {
    return lengthAreaU_;
  }
  inline CoinBigIndex numberElementsL (  ) const {
    return lengthL_;
  }
  inline CoinBigIndex lengthAreaL (  ) const {
    return lengthAreaL_;
  }
  inline CoinBigIndex numberElementsR (  ) const {
    return lengthR_;
  }
  inline CoinBigIndex numberCompressions() const
  { return numberCompressions_;}
  inline int * numberInRow() const
  { return numberInRow_.array();}
  inline int * numberInColumn() const
  { return numberInColumn_.array();}
  inline double * elementU() const
  { return elementU_.array();}
  inline int * indexRowU() const
  { return indexRowU_.array();}
  inline CoinBigIndex * startColumnU() const
  { return startColumnU_.array();}
  inline int maximumColumnsExtra()
  { return maximumColumnsExtra_;}
  inline int biasLU() const
  { return biasLU_;}
  inline void setBiasLU(int value)
  { biasLU_=value;}
  inline int persistenceFlag() const
  { return persistenceFlag_;}
  void setPersistenceFlag(int value);


  int replaceColumn ( CoinIndexedVector * regionSparse,
                      int pivotRow,
                      double pivotCheck ,
                      bool checkBeforeModifying=false);

  int updateColumnFT ( CoinIndexedVector * regionSparse,
                       CoinIndexedVector * regionSparse2);
  int updateColumn ( CoinIndexedVector * regionSparse,
                     CoinIndexedVector * regionSparse2,
                     bool noPermute=false) const;
  int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
                           CoinIndexedVector * regionSparse2,
                           CoinIndexedVector * regionSparse3,
                           bool noPermuteRegion3=false) ;
  int updateColumnTranspose ( CoinIndexedVector * regionSparse,
                              CoinIndexedVector * regionSparse2) const;
  void goSparse();
  inline int sparseThreshold ( ) const
  { return sparseThreshold_;}
  void sparseThreshold ( int value );



  inline void clearArrays()
  { gutsOfDestructor();}


  int add ( CoinBigIndex numberElements,
               int indicesRow[],
               int indicesColumn[], double elements[] );

  int addColumn ( CoinBigIndex numberElements,
                     int indicesRow[], double elements[] );

  int addRow ( CoinBigIndex numberElements,
                  int indicesColumn[], double elements[] );

  int deleteColumn ( int Row );
  int deleteRow ( int Row );

  int replaceRow ( int whichRow, int numberElements,
                      const int indicesColumn[], const double elements[] );
  void emptyRows(int numberToEmpty, const int which[]);


  void checkSparse();
  inline bool collectStatistics() const
  { return collectStatistics_;}
  inline void setCollectStatistics(bool onOff) const
  { collectStatistics_ = onOff;}
  void gutsOfDestructor(int type=1);
  void gutsOfInitialize(int type);
  void gutsOfCopy(const CoinFactorization &other);

  void resetStatistics();




  void getAreas ( int numberRows,
                  int numberColumns,
                  CoinBigIndex maximumL,
                  CoinBigIndex maximumU );

  void preProcess ( int state,
                    int possibleDuplicates = -1 );
  int factor (  );
protected:
  int factorSparse (  );
  int factorSparseSmall (  );
  int factorSparseLarge (  );
  int factorDense (  );

  bool pivotOneOtherRow ( int pivotRow,
                          int pivotColumn );
  bool pivotRowSingleton ( int pivotRow,
                           int pivotColumn );
  bool pivotColumnSingleton ( int pivotRow,
                              int pivotColumn );

  bool getColumnSpace ( int iColumn,
                        int extraNeeded );

  bool getColumnSpaceIterateR ( int iColumn, double value,
                               int iRow);
  CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
                               int iRow);
  bool getRowSpace ( int iRow, int extraNeeded );

  bool getRowSpaceIterate ( int iRow,
                            int extraNeeded );
  void checkConsistency (  );
  inline void addLink ( int index, int count ) {
    int *nextCount = nextCount_.array();
    int *firstCount = firstCount_.array();
    int *lastCount = lastCount_.array();
    int next = firstCount[count];
      lastCount[index] = -2 - count;
    if ( next < 0 ) {
      //first with that count
      firstCount[count] = index;
      nextCount[index] = -1;
    } else {
      firstCount[count] = index;
      nextCount[index] = next;
      lastCount[next] = index;
  }}
  inline void deleteLink ( int index ) {
    int *nextCount = nextCount_.array();
    int *firstCount = firstCount_.array();
    int *lastCount = lastCount_.array();
    int next = nextCount[index];
    int last = lastCount[index];
    if ( last >= 0 ) {
      nextCount[last] = next;
    } else {
      int count = -last - 2;

      firstCount[count] = next;
    }
    if ( next >= 0 ) {
      lastCount[next] = last;
    }
    nextCount[index] = -2;
    lastCount[index] = -2;
    return;
  }
  void separateLinks(int count,bool rowsFirst);
  void cleanup (  );

  void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
  void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
  void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
  void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;

  void updateColumnR ( CoinIndexedVector * region ) const;
  void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );

  void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;

  void updateColumnUSparse ( CoinIndexedVector * regionSparse, 
                             int * indexIn) const;
  void updateColumnUSparsish ( CoinIndexedVector * regionSparse, 
                               int * indexIn) const;
  int updateColumnUDensish ( double * COIN_RESTRICT region, 
                             int * COIN_RESTRICT regionIndex) const;
  void updateTwoColumnsUDensish (
                                 int & numberNonZero1,
                                 double * COIN_RESTRICT region1, 
                                 int * COIN_RESTRICT index1,
                                 int & numberNonZero2,
                                 double * COIN_RESTRICT region2, 
                                 int * COIN_RESTRICT index2) const;
  void updateColumnPFI ( CoinIndexedVector * regionSparse) const; 
  void permuteBack ( CoinIndexedVector * regionSparse, 
                     CoinIndexedVector * outVector) const;

  void updateColumnTransposePFI ( CoinIndexedVector * region) const;
  void updateColumnTransposeU ( CoinIndexedVector * region,
                                int smallestIndex) const;
  void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
                                        int smallestIndex) const;
  void updateColumnTransposeUDensish ( CoinIndexedVector * region,
                                       int smallestIndex) const;
  void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;

  void updateColumnTransposeR ( CoinIndexedVector * region ) const;
  void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
  void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;

  void updateColumnTransposeL ( CoinIndexedVector * region ) const;
  void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
  void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
  void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
  void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
public:
  int replaceColumnPFI ( CoinIndexedVector * regionSparse,
                         int pivotRow, double alpha);
protected:
  int checkPivot(double saveFromU, double oldPivot) const;
  /********************************* START LARGE TEMPLATE ********/
#ifdef INT_IS_8
#define COINFACTORIZATION_BITS_PER_INT 64
#define COINFACTORIZATION_SHIFT_PER_INT 6
#define COINFACTORIZATION_MASK_PER_INT 0x3f
#else
#define COINFACTORIZATION_BITS_PER_INT 32
#define COINFACTORIZATION_SHIFT_PER_INT 5
#define COINFACTORIZATION_MASK_PER_INT 0x1f
#endif
  template <class T>  inline bool
  pivot ( int pivotRow,
          int pivotColumn,
          CoinBigIndex pivotRowPosition,
          CoinBigIndex pivotColumnPosition,
          double work[],
          unsigned int workArea2[],
          int increment,
          int increment2,
          T markRow[] ,
          int largeInteger)
{
  int *indexColumnU = indexColumnU_.array();
  CoinBigIndex *startColumnU = startColumnU_.array();
  int *numberInColumn = numberInColumn_.array();
  double *elementU = elementU_.array();
  int *indexRowU = indexRowU_.array();
  CoinBigIndex *startRowU = startRowU_.array();
  int *numberInRow = numberInRow_.array();
  double *elementL = elementL_.array();
  int *indexRowL = indexRowL_.array();
  int *saveColumn = saveColumn_.array();
  int *nextRow = nextRow_.array();
  int *lastRow = lastRow_.array() ;

  //store pivot columns (so can easily compress)
  int numberInPivotRow = numberInRow[pivotRow] - 1;
  CoinBigIndex startColumn = startColumnU[pivotColumn];
  int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
  CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
  int put = 0;
  CoinBigIndex startRow = startRowU[pivotRow];
  CoinBigIndex endRow = startRow + numberInPivotRow + 1;

  if ( pivotColumnPosition < 0 ) {
    for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
      int iColumn = indexColumnU[pivotColumnPosition];
      if ( iColumn != pivotColumn ) {
        saveColumn[put++] = iColumn;
      } else {
        break;
      }
    }
  } else {
    for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
      saveColumn[put++] = indexColumnU[i];
    }
  }
  assert (pivotColumnPosition<endRow);
  assert (indexColumnU[pivotColumnPosition]==pivotColumn);
  pivotColumnPosition++;
  for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    saveColumn[put++] = indexColumnU[pivotColumnPosition];
  }
  //take out this bit of indexColumnU
  int next = nextRow[pivotRow];
  int last = lastRow[pivotRow];

  nextRow[last] = next;
  lastRow[next] = last;
  nextRow[pivotRow] = numberGoodU_;     //use for permute
  lastRow[pivotRow] = -2;
  numberInRow[pivotRow] = 0;
  //store column in L, compress in U and take column out
  CoinBigIndex l = lengthL_;

  if ( l + numberInPivotColumn > lengthAreaL_ ) {
    //need more memory
    printf("more memory needed in middle of invert\n");
    return false;
  }
  //l+=currentAreaL_->elementByColumn-elementL;
  CoinBigIndex lSave = l;

  pivotRowL_.array()[numberGoodL_] = pivotRow;
  CoinBigIndex * startColumnL = startColumnL_.array();
  startColumnL[numberGoodL_] = l;       //for luck and first time
  numberGoodL_++;
  startColumnL[numberGoodL_] = l + numberInPivotColumn;
  lengthL_ += numberInPivotColumn;
  if ( pivotRowPosition < 0 ) {
    for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
      int iRow = indexRowU[pivotRowPosition];
      if ( iRow != pivotRow ) {
        indexRowL[l] = iRow;
        elementL[l] = elementU[pivotRowPosition];
        markRow[iRow] = static_cast<T>(l - lSave);
        l++;
        //take out of row list
        CoinBigIndex start = startRowU[iRow];
        CoinBigIndex end = start + numberInRow[iRow];
        CoinBigIndex where = start;

        while ( indexColumnU[where] != pivotColumn ) {
          where++;
        }                       /* endwhile */
#if DEBUG_COIN
        if ( where >= end ) {
          abort (  );
        }
#endif
        indexColumnU[where] = indexColumnU[end - 1];
        numberInRow[iRow]--;
      } else {
        break;
      }
    }
  } else {
    CoinBigIndex i;

    for ( i = startColumn; i < pivotRowPosition; i++ ) {
      int iRow = indexRowU[i];

      markRow[iRow] = static_cast<T>(l - lSave);
      indexRowL[l] = iRow;
      elementL[l] = elementU[i];
      l++;
      //take out of row list
      CoinBigIndex start = startRowU[iRow];
      CoinBigIndex end = start + numberInRow[iRow];
      CoinBigIndex where = start;

      while ( indexColumnU[where] != pivotColumn ) {
        where++;
      }                         /* endwhile */
#if DEBUG_COIN
      if ( where >= end ) {
        abort (  );
      }
#endif
      indexColumnU[where] = indexColumnU[end - 1];
      numberInRow[iRow]--;
      assert (numberInRow[iRow]>=0);
    }
  }
  assert (pivotRowPosition<endColumn);
  assert (indexRowU[pivotRowPosition]==pivotRow);
  double pivotElement = elementU[pivotRowPosition];
  double pivotMultiplier = 1.0 / pivotElement;

  pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
  pivotRowPosition++;
  for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
    int iRow = indexRowU[pivotRowPosition];
    
    markRow[iRow] = static_cast<T>(l - lSave);
    indexRowL[l] = iRow;
    elementL[l] = elementU[pivotRowPosition];
    l++;
    //take out of row list
    CoinBigIndex start = startRowU[iRow];
    CoinBigIndex end = start + numberInRow[iRow];
    CoinBigIndex where = start;
    
    while ( indexColumnU[where] != pivotColumn ) {
      where++;
    }                           /* endwhile */
#if DEBUG_COIN
    if ( where >= end ) {
      abort (  );
    }
#endif
    indexColumnU[where] = indexColumnU[end - 1];
    numberInRow[iRow]--;
    assert (numberInRow[iRow]>=0);
  }
  markRow[pivotRow] = static_cast<T>(largeInteger);
  //compress pivot column (move pivot to front including saved)
  numberInColumn[pivotColumn] = 0;
  //use end of L for temporary space
  int *indexL = &indexRowL[lSave];
  double *multipliersL = &elementL[lSave];

  //adjust
  int j;

  for ( j = 0; j < numberInPivotColumn; j++ ) {
    multipliersL[j] *= pivotMultiplier;
  }
  //zero out fill
  CoinBigIndex iErase;
  for ( iErase = 0; iErase < increment2 * numberInPivotRow;
        iErase++ ) {
    workArea2[iErase] = 0;
  }
  CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
  unsigned int *temp2 = workArea2;
  int * nextColumn = nextColumn_.array();

  //pack down and move to work
  int jColumn;
  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
    int iColumn = saveColumn[jColumn];
    CoinBigIndex startColumn = startColumnU[iColumn];
    CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
    int iRow = indexRowU[startColumn];
    double value = elementU[startColumn];
    double largest;
    CoinBigIndex put = startColumn;
    CoinBigIndex positionLargest = -1;
    double thisPivotValue = 0.0;

    //compress column and find largest not updated
    bool checkLargest;
    int mark = markRow[iRow];

    if ( mark == largeInteger+1 ) {
      largest = fabs ( value );
      positionLargest = put;
      put++;
      checkLargest = false;
    } else {
      //need to find largest
      largest = 0.0;
      checkLargest = true;
      if ( mark != largeInteger ) {
        //will be updated
        work[mark] = value;
        int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
        int bit = mark & COINFACTORIZATION_MASK_PER_INT;

        temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
        added--;
      } else {
        thisPivotValue = value;
      }
    }
    CoinBigIndex i;
    for ( i = startColumn + 1; i < endColumn; i++ ) {
      iRow = indexRowU[i];
      value = elementU[i];
      int mark = markRow[iRow];

      if ( mark == largeInteger+1 ) {
        //keep
        indexRowU[put] = iRow;
        elementU[put] = value;
        if ( checkLargest ) {
          double absValue = fabs ( value );

          if ( absValue > largest ) {
            largest = absValue;
            positionLargest = put;
          }
        }
        put++;
      } else if ( mark != largeInteger ) {
        //will be updated
        work[mark] = value;
        int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
        int bit = mark & COINFACTORIZATION_MASK_PER_INT;

        temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
        added--;
      } else {
        thisPivotValue = value;
      }
    }
    //slot in pivot
    elementU[put] = elementU[startColumn];
    indexRowU[put] = indexRowU[startColumn];
    if ( positionLargest == startColumn ) {
      positionLargest = put;    //follow if was largest
    }
    put++;
    elementU[startColumn] = thisPivotValue;
    indexRowU[startColumn] = pivotRow;
    //clean up counts
    startColumn++;
    numberInColumn[iColumn] = put - startColumn;
    int * numberInColumnPlus = numberInColumnPlus_.array();
    numberInColumnPlus[iColumn]++;
    startColumnU[iColumn]++;
    //how much space have we got
    int next = nextColumn[iColumn];
    CoinBigIndex space;

    space = startColumnU[next] - put - numberInColumnPlus[next];
    //assume no zero elements
    if ( numberInPivotColumn > space ) {
      //getColumnSpace also moves fixed part
      if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
        return false;
      }
      //redo starts
      positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
      startColumn = startColumnU[iColumn];
      put = startColumn + numberInColumn[iColumn];
    }
    double tolerance = zeroTolerance_;

    int *nextCount = nextCount_.array();
    for ( j = 0; j < numberInPivotColumn; j++ ) {
      value = work[j] - thisPivotValue * multipliersL[j];
      double absValue = fabs ( value );

      if ( absValue > tolerance ) {
        work[j] = 0.0;
        elementU[put] = value;
        indexRowU[put] = indexL[j];
        if ( absValue > largest ) {
          largest = absValue;
          positionLargest = put;
        }
        put++;
      } else {
        work[j] = 0.0;
        added--;
        int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
        int bit = j & COINFACTORIZATION_MASK_PER_INT;

        if ( temp2[word] & ( 1 << bit ) ) {
          //take out of row list
          iRow = indexL[j];
          CoinBigIndex start = startRowU[iRow];
          CoinBigIndex end = start + numberInRow[iRow];
          CoinBigIndex where = start;

          while ( indexColumnU[where] != iColumn ) {
            where++;
          }                     /* endwhile */
#if DEBUG_COIN
          if ( where >= end ) {
            abort (  );
          }
#endif
          indexColumnU[where] = indexColumnU[end - 1];
          numberInRow[iRow]--;
        } else {
          //make sure won't be added
          int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
          int bit = j & COINFACTORIZATION_MASK_PER_INT;

          temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
        }
      }
    }
    numberInColumn[iColumn] = put - startColumn;
    //move largest
    if ( positionLargest >= 0 ) {
      value = elementU[positionLargest];
      iRow = indexRowU[positionLargest];
      elementU[positionLargest] = elementU[startColumn];
      indexRowU[positionLargest] = indexRowU[startColumn];
      elementU[startColumn] = value;
      indexRowU[startColumn] = iRow;
    }
    //linked list for column
    if ( nextCount[iColumn + numberRows_] != -2 ) {
      //modify linked list
      deleteLink ( iColumn + numberRows_ );
      addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
    }
    temp2 += increment2;
  }
  //get space for row list
  unsigned int *putBase = workArea2;
  int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
  int i = 0;

  // do linked lists and update counts
  while ( bigLoops ) {
    bigLoops--;
    int bit;
    for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
      unsigned int *putThis = putBase;
      int iRow = indexL[i];

      //get space
      int number = 0;
      int jColumn;

      for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
        unsigned int test = *putThis;

        putThis += increment2;
        test = 1 - ( ( test >> bit ) & 1 );
        number += test;
      }
      int next = nextRow[iRow];
      CoinBigIndex space;

      space = startRowU[next] - startRowU[iRow];
      number += numberInRow[iRow];
      if ( space < number ) {
        if ( !getRowSpace ( iRow, number ) ) {
          return false;
        }
      }
      // now do
      putThis = putBase;
      next = nextRow[iRow];
      number = numberInRow[iRow];
      CoinBigIndex end = startRowU[iRow] + number;
      int saveIndex = indexColumnU[startRowU[next]];

      //add in
      for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
        unsigned int test = *putThis;

        putThis += increment2;
        test = 1 - ( ( test >> bit ) & 1 );
        indexColumnU[end] = saveColumn[jColumn];
        end += test;
      }
      //put back next one in case zapped
      indexColumnU[startRowU[next]] = saveIndex;
      markRow[iRow] = static_cast<T>(largeInteger+1);
      number = end - startRowU[iRow];
      numberInRow[iRow] = number;
      deleteLink ( iRow );
      addLink ( iRow, number );
    }
    putBase++;
  }                             /* endwhile */
  int bit;

  for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
    unsigned int *putThis = putBase;
    int iRow = indexL[i];

    //get space
    int number = 0;
    int jColumn;

    for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
      unsigned int test = *putThis;

      putThis += increment2;
      test = 1 - ( ( test >> bit ) & 1 );
      number += test;
    }
    int next = nextRow[iRow];
    CoinBigIndex space;

    space = startRowU[next] - startRowU[iRow];
    number += numberInRow[iRow];
    if ( space < number ) {
      if ( !getRowSpace ( iRow, number ) ) {
        return false;
      }
    }
    // now do
    putThis = putBase;
    next = nextRow[iRow];
    number = numberInRow[iRow];
    CoinBigIndex end = startRowU[iRow] + number;
    int saveIndex;

    saveIndex = indexColumnU[startRowU[next]];

    //add in
    for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
      unsigned int test = *putThis;

      putThis += increment2;
      test = 1 - ( ( test >> bit ) & 1 );

      indexColumnU[end] = saveColumn[jColumn];
      end += test;
    }
    indexColumnU[startRowU[next]] = saveIndex;
    markRow[iRow] = static_cast<T>(largeInteger+1);
    number = end - startRowU[iRow];
    numberInRow[iRow] = number;
    deleteLink ( iRow );
    addLink ( iRow, number );
  }
  markRow[pivotRow] = static_cast<T>(largeInteger+1);
  //modify linked list for pivots
  deleteLink ( pivotRow );
  deleteLink ( pivotColumn + numberRows_ );
  totalElements_ += added;
  return true;
}

  /********************************* END LARGE TEMPLATE ********/

protected:


  double pivotTolerance_;
  double zeroTolerance_;
#ifndef COIN_FAST_CODE
  double slackValue_;
#else
#ifndef slackValue_
#define slackValue_ -1.0
#endif
#endif
  double areaFactor_;
  double relaxCheck_;
  int numberRows_;
  int numberRowsExtra_;
  int maximumRowsExtra_;
  int numberColumns_;
  int numberColumnsExtra_;
  int maximumColumnsExtra_;
  int numberGoodU_;
  int numberGoodL_;
  int maximumPivots_;
  int numberPivots_;
  CoinBigIndex totalElements_;
  CoinBigIndex factorElements_;
  CoinIntArrayWithLength pivotColumn_;
  CoinIntArrayWithLength permute_;
  CoinIntArrayWithLength permuteBack_;
  CoinIntArrayWithLength pivotColumnBack_;
  int status_;

  //int increasingRows_;

  int numberTrials_;
  CoinBigIndexArrayWithLength startRowU_;

  CoinIntArrayWithLength numberInRow_;

  CoinIntArrayWithLength numberInColumn_;

  CoinIntArrayWithLength numberInColumnPlus_;

  CoinIntArrayWithLength firstCount_;

  CoinIntArrayWithLength nextCount_;

  CoinIntArrayWithLength lastCount_;

  CoinIntArrayWithLength nextColumn_;

  CoinIntArrayWithLength lastColumn_;

  CoinIntArrayWithLength nextRow_;

  CoinIntArrayWithLength lastRow_;

  CoinIntArrayWithLength saveColumn_;

  CoinIntArrayWithLength markRow_;

  int messageLevel_;

  int biggerDimension_;

  CoinIntArrayWithLength indexColumnU_;

  CoinIntArrayWithLength pivotRowL_;

  CoinDoubleArrayWithLength pivotRegion_;

  int numberSlacks_;

  int numberU_;

  CoinBigIndex maximumU_;

  //int baseU_;

  CoinBigIndex lengthU_;

  CoinBigIndex lengthAreaU_;

  CoinDoubleArrayWithLength elementU_;

  CoinIntArrayWithLength indexRowU_;

  CoinBigIndexArrayWithLength startColumnU_;

  CoinBigIndexArrayWithLength convertRowToColumnU_;

  CoinBigIndex numberL_;

  CoinBigIndex baseL_;

  CoinBigIndex lengthL_;

  CoinBigIndex lengthAreaL_;

  CoinDoubleArrayWithLength elementL_;

  CoinIntArrayWithLength indexRowL_;

  CoinBigIndexArrayWithLength startColumnL_;

  bool doForrestTomlin_;

  int numberR_;

  CoinBigIndex lengthR_;

  CoinBigIndex lengthAreaR_;

  double *elementR_;

  int *indexRowR_;

  CoinBigIndexArrayWithLength startColumnR_;

  double  * denseArea_;

  int * densePermute_;

  int numberDense_;

  int denseThreshold_;

  CoinDoubleArrayWithLength workArea_;

  CoinUnsignedIntArrayWithLength workArea2_;

  CoinBigIndex numberCompressions_;

  mutable double ftranCountInput_;
  mutable double ftranCountAfterL_;
  mutable double ftranCountAfterR_;
  mutable double ftranCountAfterU_;
  mutable double btranCountInput_;
  mutable double btranCountAfterU_;
  mutable double btranCountAfterR_;
  mutable double btranCountAfterL_;

  mutable int numberFtranCounts_;
  mutable int numberBtranCounts_;

  double ftranAverageAfterL_;
  double ftranAverageAfterR_;
  double ftranAverageAfterU_;
  double btranAverageAfterU_;
  double btranAverageAfterR_;
  double btranAverageAfterL_;

  mutable bool collectStatistics_;

  int sparseThreshold_;

  int sparseThreshold2_;

  CoinBigIndexArrayWithLength startRowL_;

  CoinIntArrayWithLength indexColumnL_;

  CoinDoubleArrayWithLength elementByRowL_;

  mutable CoinIntArrayWithLength sparse_;
  int biasLU_;
  int persistenceFlag_;
};
// Dense coding
#ifdef COIN_HAS_LAPACK
#define DENSE_CODE 1
/* Type of Fortran integer translated into C */
#ifndef ipfint
//typedef ipfint FORTRAN_INTEGER_TYPE ;
typedef int ipfint;
typedef const int cipfint;
#endif
#endif
#endif
// Extra for ugly include
#ifdef UGLY_COIN_FACTOR_CODING
#define FAC_UNSET (FAC_SET+1)
{
  goodPivot=false;
  //store pivot columns (so can easily compress)
  CoinBigIndex startColumnThis = startColumn[iPivotColumn];
  CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
  int put = 0;
  CoinBigIndex startRowThis = startRow[iPivotRow];
  CoinBigIndex endRow = startRowThis + numberDoRow + 1;
  if ( pivotColumnPosition < 0 ) {
    for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
      int iColumn = indexColumn[pivotColumnPosition];
      if ( iColumn != iPivotColumn ) {
        saveColumn[put++] = iColumn;
      } else {
        break;
      }
    }
  } else {
    for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
      saveColumn[put++] = indexColumn[i];
    }
  }
  assert (pivotColumnPosition<endRow);
  assert (indexColumn[pivotColumnPosition]==iPivotColumn);
  pivotColumnPosition++;
  for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
    saveColumn[put++] = indexColumn[pivotColumnPosition];
  }
  //take out this bit of indexColumn
  int next = nextRow[iPivotRow];
  int last = lastRow[iPivotRow];
  
  nextRow[last] = next;
  lastRow[next] = last;
  nextRow[iPivotRow] = numberGoodU_;    //use for permute
  lastRow[iPivotRow] = -2;
  numberInRow[iPivotRow] = 0;
  //store column in L, compress in U and take column out
  CoinBigIndex l = lengthL_;
  // **** HORRID coding coming up but a goto seems best!
  {
    if ( l + numberDoColumn > lengthAreaL_ ) {
      //need more memory
      printf("more memory needed in middle of invert\n");
      goto BAD_PIVOT;
    }
    //l+=currentAreaL_->elementByColumn-elementL;
    CoinBigIndex lSave = l;
    
    pivotRowL_.array()[numberGoodL_] = iPivotRow;