summaryrefslogtreecommitdiff
path: root/thirdparty/windows/include/coin/CoinFactorization.hpp
diff options
context:
space:
mode:
authorharpreet2016-07-01 00:11:55 +0530
committerharpreet2016-07-01 00:11:55 +0530
commit89b96cb79cbe75ba8c1afea61b3caca37a083f62 (patch)
tree1b879fb0dce2406c1c49bf5bdabcedfb12ff020c /thirdparty/windows/include/coin/CoinFactorization.hpp
parent10e2e4d8b4a7592a8631ddac8e8d1664d6f0b9e3 (diff)
downloadFOSSEE-Optimization-toolbox-89b96cb79cbe75ba8c1afea61b3caca37a083f62.tar.gz
FOSSEE-Optimization-toolbox-89b96cb79cbe75ba8c1afea61b3caca37a083f62.tar.bz2
FOSSEE-Optimization-toolbox-89b96cb79cbe75ba8c1afea61b3caca37a083f62.zip
Windows Included
Diffstat (limited to 'thirdparty/windows/include/coin/CoinFactorization.hpp')
-rw-r--r--thirdparty/windows/include/coin/CoinFactorization.hpp2044
1 files changed, 2044 insertions, 0 deletions
diff --git a/thirdparty/windows/include/coin/CoinFactorization.hpp b/thirdparty/windows/include/coin/CoinFactorization.hpp
new file mode 100644
index 0000000..0a532bf
--- /dev/null
+++ b/thirdparty/windows/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