cbcintlinprog added

1 files changed, 2044 insertions, 0 deletions

diff --git a/thirdparty/linux/include/coin/CoinFactorization.hpp b/thirdparty/linux/include/coin/CoinFactorization.hpp
new file mode 100644
index 0000000..0a532bf
--- /dev/null
+++ b/thirdparty/linux/include/coin/CoinFactorization.hpp
@@ -0,0 +1,2044 @@
+/* $Id: CoinFactorization.hpp 1767 2015-01-05 12:36:13Z forrest $ */
+// Copyright (C) 2002, International Business Machines
+// Corporation and others.  All Rights Reserved.
+// This code is licensed under the terms of the Eclipse Public License (EPL).
+
+/* 
+   Authors
+   
+   John Forrest
+
+ */
+#ifndef CoinFactorization_H
+#define CoinFactorization_H
+//#define COIN_ONE_ETA_COPY 100
+
+#include <iostream>
+#include <string>
+#include <cassert>
+#include <cstdio>
+#include <cmath>
+#include "CoinTypes.hpp"
+#include "CoinIndexedVector.hpp"
+
+class CoinPackedMatrix;
+/** This deals with Factorization and Updates
+
+    This class started with a parallel simplex code I was writing in the
+    mid 90's.  The need for parallelism led to many complications and
+    I have simplified as much as I could to get back to this.
+
+    I was aiming at problems where I might get speed-up so I was looking at dense
+    problems or ones with structure.  This led to permuting input and output
+    vectors and to increasing the number of rows each rank-one update.  This is 
+    still in as a minor overhead.
+
+    I have also put in handling for hyper-sparsity.  I have taken out
+    all outer loop unrolling, dense matrix handling and most of the
+    book-keeping for slacks.  Also I always use FTRAN approach to updating
+    even if factorization fairly dense.  All these could improve performance.
+
+    I blame some of the coding peculiarities on the history of the code
+    but mostly it is just because I can't do elegant code (or useful
+    comments).
+
+    I am assuming that 32 bits is enough for number of rows or columns, but CoinBigIndex
+    may be redefined to get 64 bits.
+ */
+
+
+class CoinFactorization {
+   friend void CoinFactorizationUnitTest( const std::string & mpsDir );
+
+public:
+
+  /**@name Constructors and destructor and copy */
+  //@{
+  /// Default constructor
+    CoinFactorization (  );
+  /// Copy constructor 
+  CoinFactorization ( const CoinFactorization &other);
+
+  /// Destructor
+   ~CoinFactorization (  );
+  /// Delete all stuff (leaves as after CoinFactorization())
+  void almostDestructor();
+  /// Debug show object (shows one representation)
+  void show_self (  ) const;
+  /// Debug - save on file - 0 if no error
+  int saveFactorization (const char * file  ) const;
+  /** Debug - restore from file - 0 if no error on file.
+      If factor true then factorizes as if called from ClpFactorization
+  */
+  int restoreFactorization (const char * file  , bool factor=false) ;
+  /// Debug - sort so can compare
+  void sort (  ) const;
+  /// = copy
+    CoinFactorization & operator = ( const CoinFactorization & other );
+  //@}
+
+  /**@name Do factorization */
+  //@{
+  /** When part of LP - given by basic variables.
+  Actually does factorization.
+  Arrays passed in have non negative value to say basic.
+  If status is okay, basic variables have pivot row - this is only needed
+  If status is singular, then basic variables have pivot row
+  and ones thrown out have -1
+  returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
+  int factorize ( const CoinPackedMatrix & matrix, 
+		  int rowIsBasic[], int columnIsBasic[] , 
+		  double areaFactor = 0.0 );
+  /** When given as triplets.
+  Actually does factorization.  maximumL is guessed maximum size of L part of
+  final factorization, maximumU of U part.  These are multiplied by
+  areaFactor which can be computed by user or internally.  
+  Arrays are copied in.  I could add flag to delete arrays to save a 
+  bit of memory.
+  If status okay, permutation has pivot rows - this is only needed
+  If status is singular, then basic variables have pivot row
+  and ones thrown out have -1
+  returns 0 -okay, -1 singular, -99 memory */
+  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);
+  /** Two part version for maximum flexibility
+      This part creates arrays for user to fill.
+      estimateNumberElements is safe estimate of number
+      returns 0 -okay, -99 memory */
+  int factorizePart1 ( int numberRows,
+		       int numberColumns,
+		       CoinBigIndex estimateNumberElements,
+		       int * indicesRow[],
+		       int * indicesColumn[],
+		       CoinFactorizationDouble * elements[],
+		  double areaFactor = 0.0);
+  /** This is part two of factorization
+      Arrays belong to factorization and were returned by part 1
+      If status okay, permutation has pivot rows - this is only needed
+      If status is singular, then basic variables have pivot row
+      and ones thrown out have -1
+      returns 0 -okay, -1 singular, -99 memory */
+  int factorizePart2 (int permutation[],int exactNumberElements);
+  /// Condition number - product of pivots after factorization
+  double conditionNumber() const;
+  
+  //@}
+
+  /**@name general stuff such as permutation or status */
+  //@{ 
+  /// Returns status
+  inline int status (  ) const {
+    return status_;
+  }
+  /// Sets status
+  inline void setStatus (  int value)
+  {  status_=value;  }
+  /// Returns number of pivots since factorization
+  inline int pivots (  ) const {
+    return numberPivots_;
+  }
+  /// Sets number of pivots since factorization
+  inline void setPivots (  int value ) 
+  { numberPivots_=value; }
+  /// Returns address of permute region
+  inline int *permute (  ) const {
+    return permute_.array();
+  }
+  /// Returns address of pivotColumn region (also used for permuting)
+  inline int *pivotColumn (  ) const {
+    return pivotColumn_.array();
+  }
+  /// Returns address of pivot region
+  inline CoinFactorizationDouble *pivotRegion (  ) const {
+    return pivotRegion_.array();
+  }
+  /// Returns address of permuteBack region
+  inline int *permuteBack (  ) const {
+    return permuteBack_.array();
+  }
+  /// Returns address of lastRow region
+  inline int *lastRow (  ) const {
+    return lastRow_.array();
+  }
+  /** Returns address of pivotColumnBack region (also used for permuting)
+      Now uses firstCount to save memory allocation */
+  inline int *pivotColumnBack (  ) const {
+    //return firstCount_.array();
+    return pivotColumnBack_.array();
+  }
+  /// Start of each row in L
+  inline CoinBigIndex * startRowL() const
+  { return startRowL_.array();}
+
+  /// Start of each column in L
+  inline CoinBigIndex * startColumnL() const
+  { return startColumnL_.array();}
+
+  /// Index of column in row for L
+  inline int * indexColumnL() const
+  { return indexColumnL_.array();}
+
+  /// Row indices of L
+  inline int * indexRowL() const
+  { return indexRowL_.array();}
+
+  /// Elements in L (row copy)
+  inline CoinFactorizationDouble * elementByRowL() const
+  { return elementByRowL_.array();}
+
+  /// Number of Rows after iterating
+  inline int numberRowsExtra (  ) const {
+    return numberRowsExtra_;
+  }
+  /// Set number of Rows after factorization
+  inline void setNumberRows(int value)
+  { numberRows_ = value; }
+  /// Number of Rows after factorization
+  inline int numberRows (  ) const {
+    return numberRows_;
+  }
+  /// Number in L
+  inline CoinBigIndex numberL() const
+  { return numberL_;}
+
+  /// Base of L
+  inline CoinBigIndex baseL() const
+  { return baseL_;}
+  /// Maximum of Rows after iterating
+  inline int maximumRowsExtra (  ) const {
+    return maximumRowsExtra_;
+  }
+  /// Total number of columns in factorization
+  inline int numberColumns (  ) const {
+    return numberColumns_;
+  }
+  /// Total number of elements in factorization
+  inline int numberElements (  ) const {
+    return totalElements_;
+  }
+  /// Length of FT vector
+  inline int numberForrestTomlin (  ) const {
+    return numberInColumn_.array()[numberColumnsExtra_];
+  }
+  /// Number of good columns in factorization
+  inline int numberGoodColumns (  ) const {
+    return numberGoodU_;
+  }
+  /// Whether larger areas needed
+  inline double areaFactor (  ) const {
+    return areaFactor_;
+  }
+  inline void areaFactor ( double value ) {
+    areaFactor_=value;
+  }
+  /// Returns areaFactor but adjusted for dense
+  double adjustedAreaFactor() const;
+  /// Allows change of pivot accuracy check 1.0 == none >1.0 relaxed
+  inline void relaxAccuracyCheck(double value)
+  { relaxCheck_ = value;}
+  inline double getAccuracyCheck() const
+  { return relaxCheck_;}
+  /// Level of detail of messages
+  inline int messageLevel (  ) const {
+    return messageLevel_ ;
+  }
+  void messageLevel (  int value );
+  /// Maximum number of pivots between factorizations
+  inline int maximumPivots (  ) const {
+    return maximumPivots_ ;
+  }
+  void maximumPivots (  int value );
+
+  /// Gets dense threshold
+  inline int denseThreshold() const 
+    { return denseThreshold_;}
+  /// Sets dense threshold
+  inline void setDenseThreshold(int value)
+    { denseThreshold_ = value;}
+  /// Pivot tolerance
+  inline double pivotTolerance (  ) const {
+    return pivotTolerance_ ;
+  }
+  void pivotTolerance (  double value );
+  /// Zero tolerance
+  inline double zeroTolerance (  ) const {
+    return zeroTolerance_ ;
+  }
+  void zeroTolerance (  double value );
+#ifndef COIN_FAST_CODE
+  /// Whether slack value is +1 or -1
+  inline double slackValue (  ) const {
+    return slackValue_ ;
+  }
+  void slackValue (  double value );
+#endif
+  /// Returns maximum absolute value in factorization
+  double maximumCoefficient() const;
+  /// true if Forrest Tomlin update, false if PFI 
+  inline bool forrestTomlin() const
+  { return doForrestTomlin_;}
+  inline void setForrestTomlin(bool value)
+  { doForrestTomlin_=value;}
+  /// True if FT update and space
+  inline bool spaceForForrestTomlin() const
+  {
+    CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
+    CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
+    return (space>=0)&&doForrestTomlin_;
+  }
+  //@}
+
+  /**@name some simple stuff */
+  //@{
+
+  /// Returns number of dense rows
+  inline int numberDense() const
+  { return numberDense_;}
+
+  /// Returns number in U area
+  inline CoinBigIndex numberElementsU (  ) const {
+    return lengthU_;
+  }
+  /// Setss number in U area
+  inline void setNumberElementsU(CoinBigIndex value)
+  { lengthU_ = value; }
+  /// Returns length of U area
+  inline CoinBigIndex lengthAreaU (  ) const {
+    return lengthAreaU_;
+  }
+  /// Returns number in L area
+  inline CoinBigIndex numberElementsL (  ) const {
+    return lengthL_;
+  }
+  /// Returns length of L area
+  inline CoinBigIndex lengthAreaL (  ) const {
+    return lengthAreaL_;
+  }
+  /// Returns number in R area
+  inline CoinBigIndex numberElementsR (  ) const {
+    return lengthR_;
+  }
+  /// Number of compressions done
+  inline CoinBigIndex numberCompressions() const
+  { return numberCompressions_;}
+  /// Number of entries in each row
+  inline int * numberInRow() const
+  { return numberInRow_.array();}
+  /// Number of entries in each column
+  inline int * numberInColumn() const
+  { return numberInColumn_.array();}
+  /// Elements of U
+  inline CoinFactorizationDouble * elementU() const
+  { return elementU_.array();}
+  /// Row indices of U
+  inline int * indexRowU() const
+  { return indexRowU_.array();}
+  /// Start of each column in U
+  inline CoinBigIndex * startColumnU() const
+  { return startColumnU_.array();}
+  /// Maximum number of Columns after iterating
+  inline int maximumColumnsExtra()
+  { return maximumColumnsExtra_;}
+  /** L to U bias
+      0 - U bias, 1 - some U bias, 2 some L bias, 3 L bias
+  */
+  inline int biasLU() const
+  { return biasLU_;}
+  inline void setBiasLU(int value)
+  { biasLU_=value;}
+  /** Array persistence flag
+      If 0 then as now (delete/new)
+      1 then only do arrays if bigger needed
+      2 as 1 but give a bit extra if bigger needed
+  */
+  inline int persistenceFlag() const
+  { return persistenceFlag_;}
+  void setPersistenceFlag(int value);
+  //@}
+
+  /**@name rank one updates which do exist */
+  //@{
+
+  /** Replaces one Column to basis,
+   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
+      If checkBeforeModifying is true will do all accuracy checks
+      before modifying factorization.  Whether to set this depends on
+      speed considerations.  You could just do this on first iteration
+      after factorization and thereafter re-factorize
+   partial update already in U */
+  int replaceColumn ( CoinIndexedVector * regionSparse,
+		      int pivotRow,
+		      double pivotCheck ,
+		      bool checkBeforeModifying=false,
+		      double acceptablePivot=1.0e-8);
+  /** Combines BtranU and delete elements
+      If deleted is NULL then delete elements
+      otherwise store where elements are
+  */
+  void replaceColumnU ( CoinIndexedVector * regionSparse,
+			CoinBigIndex * deleted,
+			int internalPivotRow);
+#ifdef ABC_USE_COIN_FACTORIZATION
+  /** returns empty fake vector carved out of existing
+      later - maybe use associated arrays */
+  CoinIndexedVector * fakeVector(CoinIndexedVector * vector,
+				 int already=0) const;
+  void deleteFakeVector(CoinIndexedVector * vector,
+			CoinIndexedVector * fakeVector) const;
+  /** Checks if can replace one Column to basis,
+      returns update alpha
+      Fills in region for use later
+      partial update already in U */
+  double checkReplacePart1 (  CoinIndexedVector * regionSparse,
+				     int pivotRow);
+  /** Checks if can replace one Column to basis,
+      returns update alpha
+      Fills in region for use later
+      partial update in vector */
+  double checkReplacePart1 (  CoinIndexedVector * regionSparse,
+				      CoinIndexedVector * partialUpdate,
+				     int pivotRow);
+  /** Checks if can replace one Column in basis,
+      returns 0=OK, 1=Probably OK, 2=singular, 3=no room, 5 max pivots */
+  int checkReplacePart2 ( int pivotRow,
+				  double btranAlpha, 
+				  double ftranAlpha, 
+				  double ftAlpha,
+				  double acceptablePivot = 1.0e-8);
+  /** Replaces one Column to basis,
+      partial update already in U */
+  void replaceColumnPart3 ( CoinIndexedVector * regionSparse,
+			    int pivotRow,
+			    double alpha );
+  /** Replaces one Column to basis,
+      partial update in vector */
+  void replaceColumnPart3 ( CoinIndexedVector * regionSparse,
+			    CoinIndexedVector * partialUpdate,
+			    int pivotRow,
+			    double alpha );
+  /** Updates one column (FTRAN) from regionSparse2
+      Tries to do FT update
+      number returned is negative if no room
+      regionSparse starts as zero and is zero at end.
+      Note - if regionSparse2 packed on input - will be packed on output
+      long regions
+  */
+  int updateColumnFT ( CoinIndexedVector & regionSparse);
+  int updateColumnFTPart1 ( CoinIndexedVector & regionSparse) ;
+  void updateColumnFTPart2 ( CoinIndexedVector & regionSparse) ;
+  /** Updates one column (FTRAN) - long region
+      Tries to do FT update
+      puts partial update in vector */
+  void updateColumnFT ( CoinIndexedVector & regionSparseFT,
+			CoinIndexedVector & partialUpdate,
+			int which);
+  /** Updates one column (FTRAN) long region */
+  int updateColumn ( CoinIndexedVector & regionSparse) const;
+  /** Updates one column (FTRAN) from regionFT
+      Tries to do FT update
+      number returned is negative if no room.
+      Also updates regionOther - long region*/
+  int updateTwoColumnsFT ( CoinIndexedVector & regionSparseFT,
+			   CoinIndexedVector & regionSparseOther);
+  /** Updates one column (BTRAN) - long region*/
+  int updateColumnTranspose ( CoinIndexedVector & regionSparse) const;
+  /** Updates one column (FTRAN) - long region */
+  void updateColumnCpu ( CoinIndexedVector & regionSparse,int whichCpu) const;
+  /** Updates one column (BTRAN) - long region */
+  void updateColumnTransposeCpu ( CoinIndexedVector & regionSparse,int whichCpu) const;
+  /** Updates one full column (FTRAN) - long region */
+  void updateFullColumn ( CoinIndexedVector & regionSparse) const;
+  /** Updates one full column (BTRAN) - long region */
+  void updateFullColumnTranspose ( CoinIndexedVector & regionSparse) const;
+  /** Updates one column for dual steepest edge weights (FTRAN) - long region */
+  void updateWeights ( CoinIndexedVector & regionSparse) const;
+  /// Returns true if wants tableauColumn in replaceColumn
+  inline bool wantsTableauColumn() const
+  {return false;}
+  /// Pivot tolerance
+  inline double minimumPivotTolerance (  ) const {
+    return pivotTolerance_ ;
+  }
+  inline void minimumPivotTolerance (  double value )
+  { pivotTolerance(value);}
+  /// Says parallel
+  inline void setParallelMode(int value)
+  { parallelMode_=value;}
+  /// Sets solve mode
+  inline void setSolveMode(int value)
+  { parallelMode_ &= 3;parallelMode_ |= (value<<2);}
+  /// Sets solve mode
+  inline int solveMode() const
+  { return parallelMode_ >> 2;}
+  /// Update partial Ftran by R update
+  void updatePartialUpdate(CoinIndexedVector & partialUpdate);
+  /// Makes a non-singular basis by replacing variables
+  void makeNonSingular(int *  COIN_RESTRICT sequence);
+#endif
+  //@}
+
+  /**@name various uses of factorization (return code number elements) 
+   which user may want to know about */
+  //@{
+  /** Updates one column (FTRAN) from regionSparse2
+      Tries to do FT update
+      number returned is negative if no room
+      regionSparse starts as zero and is zero at end.
+      Note - if regionSparse2 packed on input - will be packed on output
+  */
+  int updateColumnFT ( CoinIndexedVector * regionSparse,
+		       CoinIndexedVector * regionSparse2);
+  /** This version has same effect as above with FTUpdate==false
+      so number returned is always >=0 */
+  int updateColumn ( CoinIndexedVector * regionSparse,
+		     CoinIndexedVector * regionSparse2,
+		     bool noPermute=false) const;
+  /** Updates one column (FTRAN) from region2
+      Tries to do FT update
+      number returned is negative if no room.
+      Also updates region3
+      region1 starts as zero and is zero at end */
+  int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
+			   CoinIndexedVector * regionSparse2,
+			   CoinIndexedVector * regionSparse3,
+			   bool noPermuteRegion3=false) ;
+  /** Updates one column (BTRAN) from regionSparse2
+      regionSparse starts as zero and is zero at end 
+      Note - if regionSparse2 packed on input - will be packed on output
+  */
+  int updateColumnTranspose ( CoinIndexedVector * regionSparse,
+			      CoinIndexedVector * regionSparse2) const;
+  /** makes a row copy of L for speed and to allow very sparse problems */
+  void goSparse();
+  /**  get sparse threshold */
+  inline int sparseThreshold ( ) const
+  { return sparseThreshold_;}
+  /**  set sparse threshold */
+  void sparseThreshold ( int value );
+  //@}
+  /// *** Below this user may not want to know about
+
+  /**@name various uses of factorization (return code number elements) 
+   which user may not want to know about (left over from my LP code) */
+  //@{
+  /// Get rid of all memory
+  inline void clearArrays()
+  { gutsOfDestructor();}
+  //@}
+
+  /**@name various updates - none of which have been written! */
+  //@{
+
+  /** Adds given elements to Basis and updates factorization,
+      can increase size of basis. Returns rank */
+  int add ( CoinBigIndex numberElements,
+	       int indicesRow[],
+	       int indicesColumn[], double elements[] );
+
+  /** Adds one Column to basis,
+      can increase size of basis. Returns rank */
+  int addColumn ( CoinBigIndex numberElements,
+		     int indicesRow[], double elements[] );
+
+  /** Adds one Row to basis,
+      can increase size of basis. Returns rank */
+  int addRow ( CoinBigIndex numberElements,
+		  int indicesColumn[], double elements[] );
+
+  /// Deletes one Column from basis, returns rank
+  int deleteColumn ( int Row );
+  /// Deletes one Row from basis, returns rank
+  int deleteRow ( int Row );
+
+  /** Replaces one Row in basis,
+      At present assumes just a singleton on row is in basis
+      returns 0=OK, 1=Probably OK, 2=singular, 3 no space */
+  int replaceRow ( int whichRow, int numberElements,
+		      const int indicesColumn[], const double elements[] );
+  /// Takes out all entries for given rows
+  void emptyRows(int numberToEmpty, const int which[]);
+  //@}
+  /**@name used by ClpFactorization */
+  /// See if worth going sparse
+  void checkSparse();
+  /// For statistics
+#if 0 //def CLP_FACTORIZATION_INSTRUMENT
+  inline bool collectStatistics() const
+  { return collectStatistics_;}
+  /// For statistics 
+  inline void setCollectStatistics(bool onOff) const
+  { collectStatistics_ = onOff;}
+#else
+  inline bool collectStatistics() const
+  { return true;}
+  /// For statistics 
+  inline void setCollectStatistics(bool onOff) const
+  { }
+#endif
+  /// The real work of constructors etc 0 just scalars, 1 bit normal 
+  void gutsOfDestructor(int type=1);
+  /// 1 bit - tolerances etc, 2 more, 4 dummy arrays
+  void gutsOfInitialize(int type);
+  void gutsOfCopy(const CoinFactorization &other);
+
+  /// Reset all sparsity etc statistics
+  void resetStatistics();
+
+
+  //@}
+
+  /**@name used by factorization */
+  /// Gets space for a factorization, called by constructors
+  void getAreas ( int numberRows,
+		  int numberColumns,
+		  CoinBigIndex maximumL,
+		  CoinBigIndex maximumU );
+
+  /** PreProcesses raw triplet data.
+      state is 0 - triplets, 1 - some counts etc , 2 - .. */
+  void preProcess ( int state,
+		    int possibleDuplicates = -1 );
+  /// Does most of factorization
+  int factor (  );
+protected:
+  /** Does sparse phase of factorization
+      return code is <0 error, 0= finished */
+  int factorSparse (  );
+  /** Does sparse phase of factorization (for smaller problems)
+      return code is <0 error, 0= finished */
+  int factorSparseSmall (  );
+  /** Does sparse phase of factorization (for larger problems)
+      return code is <0 error, 0= finished */
+  int factorSparseLarge (  );
+  /** Does dense phase of factorization
+      return code is <0 error, 0= finished */
+  int factorDense (  );
+
+  /// Pivots when just one other row so faster?
+  bool pivotOneOtherRow ( int pivotRow,
+			  int pivotColumn );
+  /// Does one pivot on Row Singleton in factorization
+  bool pivotRowSingleton ( int pivotRow,
+			   int pivotColumn );
+  /// Does one pivot on Column Singleton in factorization
+  bool pivotColumnSingleton ( int pivotRow,
+			      int pivotColumn );
+
+  /** Gets space for one Column with given length,
+   may have to do compression  (returns True if successful),
+   also moves existing vector,
+   extraNeeded is over and above present */
+  bool getColumnSpace ( int iColumn,
+			int extraNeeded );
+
+  /** Reorders U so contiguous and in order (if there is space)
+      Returns true if it could */
+  bool reorderU();
+  /**  getColumnSpaceIterateR.  Gets space for one extra R element in Column
+       may have to do compression  (returns true)
+       also moves existing vector */
+  bool getColumnSpaceIterateR ( int iColumn, double value,
+			       int iRow);
+  /**  getColumnSpaceIterate.  Gets space for one extra U element in Column
+       may have to do compression  (returns true)
+       also moves existing vector.
+       Returns -1 if no memory or where element was put
+       Used by replaceRow (turns off R version) */
+  CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
+			       int iRow);
+  /** Gets space for one Row with given length,
+  may have to do compression  (returns True if successful),
+  also moves existing vector */
+  bool getRowSpace ( int iRow, int extraNeeded );
+
+  /** Gets space for one Row with given length while iterating,
+  may have to do compression  (returns True if successful),
+  also moves existing vector */
+  bool getRowSpaceIterate ( int iRow,
+			    int extraNeeded );
+  /// Checks that row and column copies look OK
+  void checkConsistency (  );
+  /// Adds a link in chain of equal counts
+  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;
+  }}
+  /// Deletes a link in chain of equal counts
+  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;
+  }
+  /// Separate out links with same row/column count
+  void separateLinks(int count,bool rowsFirst);
+  /// Cleans up at end of factorization
+  void cleanup (  );
+
+  /// Updates part of column (FTRANL)
+  void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
+  /// Updates part of column (FTRANL) when densish
+  void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
+  /// Updates part of column (FTRANL) when sparse
+  void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
+  /// Updates part of column (FTRANL) when sparsish
+  void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
+
+  /// Updates part of column (FTRANR) without FT update
+  void updateColumnR ( CoinIndexedVector * region ) const;
+  /** Updates part of column (FTRANR) with FT update.
+      Also stores update after L and R */
+  void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
+
+  /// Updates part of column (FTRANU)
+  void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
+
+  /// Updates part of column (FTRANU) when sparse
+  void updateColumnUSparse ( CoinIndexedVector * regionSparse, 
+			     int * indexIn) const;
+  /// Updates part of column (FTRANU) when sparsish
+  void updateColumnUSparsish ( CoinIndexedVector * regionSparse, 
+			       int * indexIn) const;
+  /// Updates part of column (FTRANU)
+  int updateColumnUDensish ( double * COIN_RESTRICT region, 
+			     int * COIN_RESTRICT regionIndex) const;
+  /// Updates part of 2 columns (FTRANU) real work
+  void updateTwoColumnsUDensish (
+				 int & numberNonZero1,
+				 double * COIN_RESTRICT region1, 
+				 int * COIN_RESTRICT index1,
+				 int & numberNonZero2,
+				 double * COIN_RESTRICT region2, 
+				 int * COIN_RESTRICT index2) const;
+  /// Updates part of column PFI (FTRAN) (after rest)
+  void updateColumnPFI ( CoinIndexedVector * regionSparse) const; 
+  /// Permutes back at end of updateColumn
+  void permuteBack ( CoinIndexedVector * regionSparse, 
+		     CoinIndexedVector * outVector) const;
+
+  /// Updates part of column transpose PFI (BTRAN) (before rest)
+  void updateColumnTransposePFI ( CoinIndexedVector * region) const;
+  /** Updates part of column transpose (BTRANU),
+      assumes index is sorted i.e. region is correct */
+  void updateColumnTransposeU ( CoinIndexedVector * region,
+				int smallestIndex) const;
+  /** Updates part of column transpose (BTRANU) when sparsish,
+      assumes index is sorted i.e. region is correct */
+  void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
+					int smallestIndex) const;
+  /** Updates part of column transpose (BTRANU) when densish,
+      assumes index is sorted i.e. region is correct */
+  void updateColumnTransposeUDensish ( CoinIndexedVector * region,
+				       int smallestIndex) const;
+  /** Updates part of column transpose (BTRANU) when sparse,
+      assumes index is sorted i.e. region is correct */
+  void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
+  /** Updates part of column transpose (BTRANU) by column
+      assumes index is sorted i.e. region is correct */
+  void updateColumnTransposeUByColumn ( CoinIndexedVector * region,
+					int smallestIndex) const;
+
+  /// Updates part of column transpose (BTRANR)
+  void updateColumnTransposeR ( CoinIndexedVector * region ) const;
+  /// Updates part of column transpose (BTRANR) when dense
+  void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
+  /// Updates part of column transpose (BTRANR) when sparse
+  void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
+
+  /// Updates part of column transpose (BTRANL)
+  void updateColumnTransposeL ( CoinIndexedVector * region ) const;
+  /// Updates part of column transpose (BTRANL) when densish by column
+  void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
+  /// Updates part of column transpose (BTRANL) when densish by row
+  void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
+  /// Updates part of column transpose (BTRANL) when sparsish by row
+  void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
+  /// Updates part of column transpose (BTRANL) when sparse (by Row)
+  void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
+public:
+  /** Replaces one Column to basis for PFI
+   returns 0=OK, 1=Probably OK, 2=singular, 3=no room.
+   In this case region is not empty - it is incoming variable (updated)
+  */
+  int replaceColumnPFI ( CoinIndexedVector * regionSparse,
+			 int pivotRow, double alpha);
+protected:
+  /** Returns accuracy status of replaceColumn
+      returns 0=OK, 1=Probably OK, 2=singular */
+  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,
+	  CoinFactorizationDouble work[],
+	  unsigned int workArea2[],
+	  int increment2,
+	  T markRow[] ,
+	  int largeInteger)
+{
+  int *indexColumnU = indexColumnU_.array();
+  CoinBigIndex *startColumnU = startColumnU_.array();
+  int *numberInColumn = numberInColumn_.array();
+  CoinFactorizationDouble *elementU = elementU_.array();
+  int *indexRowU = indexRowU_.array();
+  CoinBigIndex *startRowU = startRowU_.array();
+  int *numberInRow = numberInRow_.array();
+  CoinFactorizationDouble *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
+    if ((messageLevel_&4)!=0) 
+      printf("more memory needed in middle of invert\n");
+    return false;
+  }
+  //l+=currentAreaL_->elementByColumn-elementL;
+  CoinBigIndex lSave = l;
+
+  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);
+  CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
+  CoinFactorizationDouble 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];
+  CoinFactorizationDouble *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];
+    CoinFactorizationDouble value = elementU[startColumn];
+    double largest;
+    CoinBigIndex put = startColumn;
+    CoinBigIndex positionLargest = -1;
+    CoinFactorizationDouble 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
+      if (positionLargest >= 0)
+         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;
+	assert (put<lengthAreaU_); 
+	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 ********/
+  //@}
+////////////////// data //////////////////
+protected:
+
+  /**@name data */
+  //@{
+  /// Pivot tolerance
+  double pivotTolerance_;
+  /// Zero tolerance
+  double zeroTolerance_;
+#ifndef COIN_FAST_CODE
+  /// Whether slack value is  +1 or -1
+  double slackValue_;
+#else
+#ifndef slackValue_
+#define slackValue_ -1.0
+#endif
+#endif
+  /// How much to multiply areas by
+  double areaFactor_;
+  /// Relax check on accuracy in replaceColumn
+  double relaxCheck_;
+  /// Number of Rows in factorization
+  int numberRows_;
+  /// Number of Rows after iterating
+  int numberRowsExtra_;
+  /// Maximum number of Rows after iterating
+  int maximumRowsExtra_;
+  /// Number of Columns in factorization
+  int numberColumns_;
+  /// Number of Columns after iterating
+  int numberColumnsExtra_;
+  /// Maximum number of Columns after iterating
+  int maximumColumnsExtra_;
+  /// Number factorized in U (not row singletons)
+  int numberGoodU_;
+  /// Number factorized in L
+  int numberGoodL_;
+  /// Maximum number of pivots before factorization
+  int maximumPivots_;
+  /// Number pivots since last factorization
+  int numberPivots_;
+  /// Number of elements in U (to go)
+  ///       or while iterating total overall
+  CoinBigIndex totalElements_;
+  /// Number of elements after factorization
+  CoinBigIndex factorElements_;
+  /// Pivot order for each Column
+  CoinIntArrayWithLength pivotColumn_;
+  /// Permutation vector for pivot row order
+  CoinIntArrayWithLength permute_;
+  /// DePermutation vector for pivot row order
+  CoinIntArrayWithLength permuteBack_;
+  /// Inverse Pivot order for each Column
+  CoinIntArrayWithLength pivotColumnBack_;
+  /// Status of factorization
+  int status_;
+
+  /** 0 - no increasing rows - no permutations,
+   1 - no increasing rows but permutations 
+   2 - increasing rows 
+     - taken out as always 2 */
+  //int increasingRows_;
+
+  /// Number of trials before rejection
+  int numberTrials_;
+  /// Start of each Row as pointer
+  CoinBigIndexArrayWithLength startRowU_;
+
+  /// Number in each Row
+  CoinIntArrayWithLength numberInRow_;
+
+  /// Number in each Column
+  CoinIntArrayWithLength numberInColumn_;
+
+  /// Number in each Column including pivoted
+  CoinIntArrayWithLength numberInColumnPlus_;
+
+  /** First Row/Column with count of k,
+      can tell which by offset - Rows then Columns */
+  CoinIntArrayWithLength firstCount_;
+
+  /// Next Row/Column with count
+  CoinIntArrayWithLength nextCount_;
+
+  /// Previous Row/Column with count
+  CoinIntArrayWithLength lastCount_;
+
+  /// Next Column in memory order
+  CoinIntArrayWithLength nextColumn_;
+
+  /// Previous Column in memory order
+  CoinIntArrayWithLength lastColumn_;
+
+  /// Next Row in memory order
+  CoinIntArrayWithLength nextRow_;
+
+  /// Previous Row in memory order
+  CoinIntArrayWithLength lastRow_;
+
+  /// Columns left to do in a single pivot
+  CoinIntArrayWithLength saveColumn_;
+
+  /// Marks rows to be updated
+  CoinIntArrayWithLength markRow_;
+
+  /// Detail in messages
+  int messageLevel_;
+
+  /// Larger of row and column size
+  int biggerDimension_;
+
+  /// Base address for U (may change)
+  CoinIntArrayWithLength indexColumnU_;
+
+  /// Pivots for L
+  CoinIntArrayWithLength pivotRowL_;
+
+  /// Inverses of pivot values
+  CoinFactorizationDoubleArrayWithLength pivotRegion_;
+
+  /// Number of slacks at beginning of U
+  int numberSlacks_;
+
+  /// Number in U
+  int numberU_;
+
+  /// Maximum space used in U
+  CoinBigIndex maximumU_;
+
+  /// Base of U is always 0
+  //int baseU_;
+
+  /// Length of U
+  CoinBigIndex lengthU_;
+
+  /// Length of area reserved for U
+  CoinBigIndex lengthAreaU_;
+
+/// Elements of U
+  CoinFactorizationDoubleArrayWithLength elementU_;
+
+/// Row indices of U
+  CoinIntArrayWithLength indexRowU_;
+
+/// Start of each column in U
+  CoinBigIndexArrayWithLength startColumnU_;
+
+/// Converts rows to columns in U 
+  CoinBigIndexArrayWithLength convertRowToColumnU_;
+
+  /// Number in L
+  CoinBigIndex numberL_;
+
+/// Base of L
+  CoinBigIndex baseL_;
+
+  /// Length of L
+  CoinBigIndex lengthL_;
+
+  /// Length of area reserved for L
+  CoinBigIndex lengthAreaL_;
+
+  /// Elements of L
+  CoinFactorizationDoubleArrayWithLength elementL_;
+
+  /// Row indices of L
+  CoinIntArrayWithLength indexRowL_;
+
+  /// Start of each column in L
+  CoinBigIndexArrayWithLength startColumnL_;
+
+  /// true if Forrest Tomlin update, false if PFI 
+  bool doForrestTomlin_;
+
+  /// Number in R
+  int numberR_;
+
+  /// Length of R stuff
+  CoinBigIndex lengthR_;
+
+  /// length of area reserved for R
+  CoinBigIndex lengthAreaR_;
+
+  /// Elements of R
+  CoinFactorizationDouble *elementR_;
+
+  /// Row indices for R
+  int *indexRowR_;
+
+  /// Start of columns for R
+  CoinBigIndexArrayWithLength startColumnR_;
+
+  /// Dense area
+  double  * denseArea_;
+
+  /// Dense area - actually used (for alignment etc)
+  double  * denseAreaAddress_;
+
+  /// Dense permutation
+  int * densePermute_;
+
+  /// Number of dense rows
+  int numberDense_;
+
+  /// Dense threshold
+  int denseThreshold_;
+
+  /// First work area
+  CoinFactorizationDoubleArrayWithLength workArea_;
+
+  /// Second work area
+  CoinUnsignedIntArrayWithLength workArea2_;
+
+  /// Number of compressions done
+  CoinBigIndex numberCompressions_;
+
+public:
+  /// Below are all to collect
+  mutable double ftranCountInput_;
+  mutable double ftranCountAfterL_;
+  mutable double ftranCountAfterR_;
+  mutable double ftranCountAfterU_;
+  mutable double btranCountInput_;
+  mutable double btranCountAfterU_;
+  mutable double btranCountAfterR_;
+  mutable double btranCountAfterL_;
+
+  /// We can roll over factorizations
+  mutable int numberFtranCounts_;
+  mutable int numberBtranCounts_;
+
+  /// While these are average ratios collected over last period
+  double ftranAverageAfterL_;
+  double ftranAverageAfterR_;
+  double ftranAverageAfterU_;
+  double btranAverageAfterU_;
+  double btranAverageAfterR_;
+  double btranAverageAfterL_;
+protected:
+
+  /// For statistics 
+#if 0
+  mutable bool collectStatistics_;
+#else
+#define collectStatistics_ 1
+#endif
+
+  /// Below this use sparse technology - if 0 then no L row copy
+  int sparseThreshold_;
+
+  /// And one for "sparsish"
+  int sparseThreshold2_;
+
+  /// Start of each row in L
+  CoinBigIndexArrayWithLength startRowL_;
+
+  /// Index of column in row for L
+  CoinIntArrayWithLength indexColumnL_;
+
+  /// Elements in L (row copy)
+  CoinFactorizationDoubleArrayWithLength elementByRowL_;
+
+  /// Sparse regions
+  mutable CoinIntArrayWithLength sparse_;
+  /** L to U bias
+      0 - U bias, 1 - some U bias, 2 some L bias, 3 L bias
+  */
+  int biasLU_;
+  /** Array persistence flag
+      If 0 then as now (delete/new)
+      1 then only do arrays if bigger needed
+      2 as 1 but give a bit extra if bigger needed
+  */
+  int persistenceFlag_;
+#ifdef ABC_USE_COIN_FACTORIZATION
+  /// Says if parallel
+  int parallelMode_;
+#endif
+  //@}
+};
+// Dense coding
+#ifdef COIN_HAS_LAPACK
+#ifndef COIN_FACTORIZATION_DENSE_CODE
+#define COIN_FACTORIZATION_DENSE_CODE 1
+#endif
+#endif
+#ifdef COIN_FACTORIZATION_DENSE_CODE
+/* 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
+      if ((messageLevel_&4)!=0) 
+	printf("more memory needed in middle of invert\n");
+      goto BAD_PIVOT;
+    }
+    //l+=currentAreaL_->elementByColumn-elementL;
+    CoinBigIndex lSave = l;
+    
+    CoinBigIndex * startColumnL = startColumnL_.array();
+    startColumnL[numberGoodL_] = l;	//for luck and first time
+    numberGoodL_++;
+    startColumnL[numberGoodL_] = l + numberDoColumn;
+    lengthL_ += numberDoColumn;
+    if ( pivotRowPosition < 0 ) {
+      for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
+	int iRow = indexRow[pivotRowPosition];
+	if ( iRow != iPivotRow ) {
+	  indexRowL[l] = iRow;
+	  elementL[l] = element[pivotRowPosition];
+	  markRow[iRow] = l - lSave;
+	  l++;
+	  //take out of row list
+	  CoinBigIndex start = startRow[iRow];
+	  CoinBigIndex end = start + numberInRow[iRow];
+	  CoinBigIndex where = start;
+	  
+	  while ( indexColumn[where] != iPivotColumn ) {
+	    where++;
+	  }			/* endwhile */
+#if DEBUG_COIN
+	  if ( where >= end ) {
+	    abort (  );
+	  }
+#endif
+	  indexColumn[where] = indexColumn[end - 1];
+	  numberInRow[iRow]--;
+	} else {
+	  break;
+	}
+      }
+    } else {
+      CoinBigIndex i;
+      
+      for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
+	int iRow = indexRow[i];
+	
+	markRow[iRow] = l - lSave;
+	indexRowL[l] = iRow;
+	elementL[l] = element[i];
+	l++;
+	//take out of row list
+	CoinBigIndex start = startRow[iRow];
+	CoinBigIndex end = start + numberInRow[iRow];
+	CoinBigIndex where = start;
+	
+	while ( indexColumn[where] != iPivotColumn ) {
+	  where++;
+	}				/* endwhile */
+#if DEBUG_COIN
+	if ( where >= end ) {
+	  abort (  );
+	}
+#endif
+	indexColumn[where] = indexColumn[end - 1];
+	numberInRow[iRow]--;
+	assert (numberInRow[iRow]>=0);
+      }
+    }
+    assert (pivotRowPosition<endColumn);
+    assert (indexRow[pivotRowPosition]==iPivotRow);
+    CoinFactorizationDouble pivotElement = element[pivotRowPosition];
+    CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
+    
+    pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
+    pivotRowPosition++;
+    for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
+      int iRow = indexRow[pivotRowPosition];
+      
+      markRow[iRow] = l - lSave;
+      indexRowL[l] = iRow;
+      elementL[l] = element[pivotRowPosition];
+      l++;
+      //take out of row list
+      CoinBigIndex start = startRow[iRow];
+      CoinBigIndex end = start + numberInRow[iRow];
+      CoinBigIndex where = start;
+      
+      while ( indexColumn[where] != iPivotColumn ) {
+	where++;
+      }				/* endwhile */
+#if DEBUG_COIN
+      if ( where >= end ) {
+	abort (  );
+      }
+#endif
+      indexColumn[where] = indexColumn[end - 1];
+      numberInRow[iRow]--;
+      assert (numberInRow[iRow]>=0);
+    }
+    markRow[iPivotRow] = FAC_SET;
+    //compress pivot column (move pivot to front including saved)
+    numberInColumn[iPivotColumn] = 0;
+    //use end of L for temporary space
+    int *indexL = &indexRowL[lSave];
+    CoinFactorizationDouble *multipliersL = &elementL[lSave];
+    
+    //adjust
+    int j;
+    
+    for ( j = 0; j < numberDoColumn; j++ ) {
+      multipliersL[j] *= pivotMultiplier;
+    }
+    //zero out fill
+    CoinBigIndex iErase;
+    for ( iErase = 0; iErase < increment2 * numberDoRow;
+	  iErase++ ) {
+      workArea2[iErase] = 0;
+    }
+    CoinBigIndex added = numberDoRow * numberDoColumn;
+    unsigned int *temp2 = workArea2;
+    int * nextColumn = nextColumn_.array();
+    
+    //pack down and move to work
+    int jColumn;
+    for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
+      int iColumn = saveColumn[jColumn];
+      CoinBigIndex startColumnThis = startColumn[iColumn];
+      CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
+      int iRow = indexRow[startColumnThis];
+      CoinFactorizationDouble value = element[startColumnThis];
+      double largest;
+      CoinBigIndex put = startColumnThis;
+      CoinBigIndex positionLargest = -1;
+      CoinFactorizationDouble thisPivotValue = 0.0;
+      
+      //compress column and find largest not updated
+      bool checkLargest;
+      int mark = markRow[iRow];
+      
+      if ( mark == FAC_UNSET ) {
+	largest = fabs ( value );
+	positionLargest = put;
+	put++;
+	checkLargest = false;
+      } else {
+	//need to find largest
+	largest = 0.0;
+	checkLargest = true;
+	if ( mark != FAC_SET ) {
+	  //will be updated
+	  workArea[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 = startColumnThis + 1; i < endColumn; i++ ) {
+	iRow = indexRow[i];
+	value = element[i];
+	int mark = markRow[iRow];
+	
+	if ( mark == FAC_UNSET ) {
+	  //keep
+	  indexRow[put] = iRow;
+	  element[put] = value;
+	  if ( checkLargest ) {
+	    double absValue = fabs ( value );
+	    
+	    if ( absValue > largest ) {
+	      largest = absValue;
+	      positionLargest = put;
+	    }
+	  }
+	  put++;
+	} else if ( mark != FAC_SET ) {
+	  //will be updated
+	  workArea[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
+      element[put] = element[startColumnThis];
+      indexRow[put] = indexRow[startColumnThis];
+      if ( positionLargest == startColumnThis ) {
+	positionLargest = put;	//follow if was largest
+      }
+      put++;
+      element[startColumnThis] = thisPivotValue;
+      indexRow[startColumnThis] = iPivotRow;
+      //clean up counts
+      startColumnThis++;
+      numberInColumn[iColumn] = put - startColumnThis;
+      int * numberInColumnPlus = numberInColumnPlus_.array();
+      numberInColumnPlus[iColumn]++;
+      startColumn[iColumn]++;
+      //how much space have we got
+      int next = nextColumn[iColumn];
+      CoinBigIndex space;
+      
+      space = startColumn[next] - put - numberInColumnPlus[next];
+      //assume no zero elements
+      if ( numberDoColumn > space ) {
+	//getColumnSpace also moves fixed part
+	if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
+	  goto BAD_PIVOT;
+	}
+	//redo starts
+	positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
+	startColumnThis = startColumn[iColumn];
+	put = startColumnThis + numberInColumn[iColumn];
+      }
+      double tolerance = zeroTolerance_;
+      
+      int *nextCount = nextCount_.array();
+      for ( j = 0; j < numberDoColumn; j++ ) {
+	value = workArea[j] - thisPivotValue * multipliersL[j];
+	double absValue = fabs ( value );
+	
+	if ( absValue > tolerance ) {
+	  workArea[j] = 0.0;
+	  element[put] = value;
+	  indexRow[put] = indexL[j];
+	  if ( absValue > largest ) {
+	    largest = absValue;
+	    positionLargest = put;
+	  }
+	  put++;
+	} else {
+	  workArea[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 = startRow[iRow];
+	    CoinBigIndex end = start + numberInRow[iRow];
+	    CoinBigIndex where = start;
+	    
+	    while ( indexColumn[where] != iColumn ) {
+	      where++;
+	    }			/* endwhile */
+#if DEBUG_COIN
+	    if ( where >= end ) {
+	      abort (  );
+	    }
+#endif
+	    indexColumn[where] = indexColumn[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 - startColumnThis;
+      //move largest
+      if ( positionLargest >= 0 ) {
+	value = element[positionLargest];
+	iRow = indexRow[positionLargest];
+	element[positionLargest] = element[startColumnThis];
+	indexRow[positionLargest] = indexRow[startColumnThis];
+	element[startColumnThis] = value;
+	indexRow[startColumnThis] = 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 = numberDoColumn >> 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 < numberDoRow; jColumn++ ) {
+	  unsigned int test = *putThis;
+	  
+	  putThis += increment2;
+	  test = 1 - ( ( test >> bit ) & 1 );
+	  number += test;
+	}
+	int next = nextRow[iRow];
+	CoinBigIndex space;
+	
+	space = startRow[next] - startRow[iRow];
+	number += numberInRow[iRow];
+	if ( space < number ) {
+	  if ( !getRowSpace ( iRow, number ) ) {
+	    goto BAD_PIVOT;
+	  }
+	}
+	// now do
+	putThis = putBase;
+	next = nextRow[iRow];
+	number = numberInRow[iRow];
+	CoinBigIndex end = startRow[iRow] + number;
+	int saveIndex = indexColumn[startRow[next]];
+	
+	//add in
+	for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
+	  unsigned int test = *putThis;
+	  
+	  putThis += increment2;
+	  test = 1 - ( ( test >> bit ) & 1 );
+	  indexColumn[end] = saveColumn[jColumn];
+	  end += test;
+	}
+	//put back next one in case zapped
+	indexColumn[startRow[next]] = saveIndex;
+	markRow[iRow] = FAC_UNSET;
+	number = end - startRow[iRow];
+	numberInRow[iRow] = number;
+	deleteLink ( iRow );
+	addLink ( iRow, number );
+      }
+      putBase++;
+    }				/* endwhile */
+    int bit;
+    
+    for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
+      unsigned int *putThis = putBase;
+      int iRow = indexL[i];
+      
+      //get space
+      int number = 0;
+      int jColumn;
+      
+      for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
+	unsigned int test = *putThis;
+	
+	putThis += increment2;
+	test = 1 - ( ( test >> bit ) & 1 );
+	number += test;
+      }
+      int next = nextRow[iRow];
+      CoinBigIndex space;
+      
+      space = startRow[next] - startRow[iRow];
+      number += numberInRow[iRow];
+      if ( space < number ) {
+	if ( !getRowSpace ( iRow, number ) ) {
+	  goto BAD_PIVOT;
+	}
+      }
+      // now do
+      putThis = putBase;
+      next = nextRow[iRow];
+      number = numberInRow[iRow];
+      CoinBigIndex end = startRow[iRow] + number;
+      int saveIndex;
+      
+      saveIndex = indexColumn[startRow[next]];
+      
+      //add in
+      for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
+	unsigned int test = *putThis;
+	
+	putThis += increment2;
+	test = 1 - ( ( test >> bit ) & 1 );
+	
+	indexColumn[end] = saveColumn[jColumn];
+	end += test;
+      }
+      indexColumn[startRow[next]] = saveIndex;
+      markRow[iRow] = FAC_UNSET;
+      number = end - startRow[iRow];
+      numberInRow[iRow] = number;
+      deleteLink ( iRow );
+      addLink ( iRow, number );
+    }
+    markRow[iPivotRow] = FAC_UNSET;
+    //modify linked list for pivots
+    deleteLink ( iPivotRow );
+    deleteLink ( iPivotColumn + numberRows_ );
+    totalElements_ += added;
+    goodPivot= true;
+    // **** UGLY UGLY UGLY
+  }
+ BAD_PIVOT:
+
+  ;
+}
+#undef FAC_UNSET
+#endif