Skip to content

Commit

Permalink
Optimizations to LU and FT basis factorizations
Browse files Browse the repository at this point in the history 
* a more efficient computation of the explicit basis:
1. Compute from left to write instead of from right to left. This
makes multipliciation by eta matrices easier
2. Compute invLP once-and-for-all, when the basis is factorized

* pivoting stats

* bug fix in store/restore

* supprt for a "basis column oracle", which allows the basis
factorization classes to query the tableau for constraint matrix columns

* cleanup in FT factorization

* grab a fresh basis from the oracle instead of condensing th etas
  • Loading branch information
@guykatzz
guykatzz authored May 9, 2018
1 parent 8b6f1ce commit 6b8125f
Show file tree
Hide file tree
Showing 15 changed files with 296 additions and 322 deletions.
6 changes: 3 additions & 3 deletions src/basis_factorization/BasisFactorizationFactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,14 +16,14 @@
#include "GlobalConfiguration.h"
#include "LUFactorization.h"

IBasisFactorization *BasisFactorizationFactory::createBasisFactorization( unsigned basisSize )
IBasisFactorization *BasisFactorizationFactory::createBasisFactorization( unsigned basisSize, const IBasisFactorization::BasisColumnOracle &basisColumnOracle )
{
if ( GlobalConfiguration::BASIS_FACTORIZATION_TYPE ==
GlobalConfiguration::LU_FACTORIZATION )
return new LUFactorization( basisSize );
return new LUFactorization( basisSize, basisColumnOracle );
else if ( GlobalConfiguration::BASIS_FACTORIZATION_TYPE ==
GlobalConfiguration::FORREST_TOMLIN_FACTORIZATION )
return new ForrestTomlinFactorization( basisSize );
return new ForrestTomlinFactorization( basisSize, basisColumnOracle );

throw BasisFactorizationError( BasisFactorizationError::UNKNOWN_BASIS_FACTORIZATION_TYPE );
}
Expand Down
2 changes: 1 addition & 1 deletion src/basis_factorization/BasisFactorizationFactory.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
class BasisFactorizationFactory
{
public:
static IBasisFactorization *createBasisFactorization( unsigned basisSize );
static IBasisFactorization *createBasisFactorization( unsigned basisSize, const IBasisFactorization::BasisColumnOracle &basisColumnOracle );
};

#endif // __BasisFactorizationFactory_h__
Expand Down
157 changes: 17 additions & 140 deletions src/basis_factorization/ForrestTomlinFactorization.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,13 +18,13 @@
#include <cstdlib>
#include <cstring>

ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m )
: _m( m )
ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m, const BasisColumnOracle &basisColumnOracle )
: IBasisFactorization( basisColumnOracle )
, _m( m )
, _B( NULL )
, _Q( m )
, _invQ( m )
, _U( NULL )
, _invLP( NULL )
, _explicitBasisAvailable( false )
, _workMatrix( NULL )
, _workVector( NULL )
Expand All @@ -43,11 +43,6 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m )
throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED,
"ForrestTomlinFactorization::U" );

_invLP = new double[m * m];
if ( !_invLP )
throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED,
"ForrestTomlinFactorization::invLPx" );

for ( unsigned i = 0; i < _m; ++i )
{
_U[i] = new EtaMatrix( _m, i );
Expand Down Expand Up @@ -81,11 +76,6 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m )
for ( unsigned row = 0; row < _m; ++row )
_B[row * _m + row] = 1.0;

// Also, invLP = I
std::fill_n( _invLP, _m * _m, 0.0 );
for ( unsigned row = 0; row < _m; ++row )
_invLP[row * _m + row] = 1.0;

// Us, Q and invQ are all initialized to the identity matrix anyway.
}

Expand All @@ -112,12 +102,6 @@ ForrestTomlinFactorization::~ForrestTomlinFactorization()
_U = NULL;
}

if ( _invLP )
{
delete[] _invLP;
_invLP = NULL;
}

if ( _workMatrix )
{
delete[] _workMatrix;
Expand Down Expand Up @@ -309,7 +293,7 @@ void ForrestTomlinFactorization::pushEtaMatrix( unsigned columnIndex, const doub

// If the number of A matrices is too great, condense them.
if ( _A.size() > GlobalConfiguration::REFACTORIZATION_THRESHOLD )
makeExplicitBasisAvailable();
refactorizeBasis();
}

void ForrestTomlinFactorization::forwardTransformation( const double *y, double *x ) const
Expand Down Expand Up @@ -533,8 +517,6 @@ void ForrestTomlinFactorization::storeFactorization( IBasisFactorization *other
for ( const auto &lp : _LP )
otherFTFactorization->_LP.append( lp->duplicate() );

memcpy( otherFTFactorization->_invLP, _invLP, sizeof(double) * _m * _m );

List<AlmostIdentityMatrix *>::iterator aIt;
for ( aIt = otherFTFactorization->_A.begin(); aIt != otherFTFactorization->_A.end(); ++aIt )
delete *aIt;
Expand Down Expand Up @@ -572,8 +554,6 @@ void ForrestTomlinFactorization::restoreFactorization( const IBasisFactorization
for ( const auto &lp : otherFTFactorization->_LP )
_LP.append( lp->duplicate() );

memcpy( _invLP, otherFTFactorization->_invLP, sizeof(double) * _m * _m );

List<AlmostIdentityMatrix *>::iterator aIt;
for ( aIt = _A.begin(); aIt != _A.end(); ++aIt )
delete *aIt;
Expand Down Expand Up @@ -691,56 +671,6 @@ void ForrestTomlinFactorization::initialLUFactorization()
for ( unsigned j = 0; j <= i; ++j )
u->_column[j] = _workMatrix[j * _m + i];
}

// Compute inv(P1)inv(L1) ... inv(Ps)inv(LS) for later.

// Initialize to the identity matrix
std::fill_n( _invLP, _m * _m, 0 );
for ( unsigned i = 0; i < _m; ++i )
_invLP[i*_m + i] = 1.0;

// Multiply by inverted Ps and Ls
for ( auto lp = _LP.rbegin(); lp != _LP.rend(); ++lp )
{
if ( (*lp)->_pair )
{
// The inverse of a general permutation matrix is its transpose.
// However, since these are permutations of just two rows, the
// transpose is equal to the original - so no need to invert.
// Multiplying on the right permutes columns.
double temp;
unsigned first = (*lp)->_pair->first;
unsigned second = (*lp)->_pair->second;

for ( unsigned i = 0; i < _m; ++i )
{
temp = _invLP[i * _m + first];
_invLP[i * _m + first] = _invLP[i * _m + second];
_invLP[i * _m + second] = temp;
}
}
else
{
// The inverse of an eta matrix is an eta matrix with the special
// column negated and divided by the diagonal entry. The only
// exception is the diagonal entry itself, which is just the
// inverse of the original diagonal entry.
EtaMatrix *eta = (*lp)->_eta;
unsigned columnIndex = eta->_columnIndex;
double etaDiagonalEntry = 1 / eta->_column[columnIndex];

for ( unsigned row = 0; row < _m; ++row )
{
for ( unsigned col = columnIndex + 1; col < _m; ++col )
{
_invLP[row * _m + columnIndex] -=
( eta->_column[col] * _invLP[row * _m + col] );
}

_invLP[row * _m + columnIndex] *= etaDiagonalEntry;
}
}
}
}

bool ForrestTomlinFactorization::explicitBasisAvailable() const
Expand All @@ -752,70 +682,8 @@ void ForrestTomlinFactorization::makeExplicitBasisAvailable()
{
if ( explicitBasisAvailable() )
return;
/*
We know that the following equation holds:
Ak...A1 * LsPs...L1P1 * B = Q * Um...U1 * invQ
Or:

B = inv( Ak...A1 * LsPs...L1P1 ) * Q * Um...U1 * invQ
= inv(P1)inv(L1) ... inv(Ps)inv(LS) * inv(A1)...inv(Ak)
* Q * Um...U1 * invQ
*/

// Start with the already-computed inv(P1)inv(L1) ... inv(Ps)inv(LS)
memcpy( _B, _invLP, sizeof(double) * _m * _m );

// Multiply by the inverse As. An inverse of an almost identity matrix
// is just that matrix with its special value negated (unless that
// matrix is diagonal.
for ( const auto &a : _A )
{
if ( a->_row == a->_column )
{
for ( unsigned j = 0; j < _m; ++j )
_B[j * _m + a->_column] *= ( 1 / a->_value );
}
else
{
for ( unsigned j = 0; j < _m; ++j )
_B[j * _m + a->_column] -= _B[j * _m + a->_row] * a->_value;
}
}

// Permute columns according to Q
for ( unsigned i = 0; i < _m; ++i )
{
unsigned columnToCopy = _invQ._ordering[i];
for ( unsigned j = 0; j < _m; ++j )
_workMatrix[j * _m + i] = _B[j * _m + columnToCopy];
}

// Multiply by Um...U1
for ( int i = _m - 1; i >= 0; --i )
{
for ( unsigned row = 0; row < _m; ++row )
{
for ( unsigned j = 0; j < _U[i]->_columnIndex; ++j )
{
_workMatrix[row * _m + _U[i]->_columnIndex] +=
( _U[i]->_column[j] * _workMatrix[row * _m + j] );
}
}
}

// Permute columns according to invQ
for ( unsigned i = 0; i < _m; ++i )
{
unsigned columnToCopy = _Q._ordering[i];
for ( unsigned j = 0; j < _m; ++j )
_B[j * _m + i] = _workMatrix[j * _m + columnToCopy];
}

clearFactorization();
initialLUFactorization();
_explicitBasisAvailable = true;
refactorizeBasis();
}

const double *ForrestTomlinFactorization::getBasis() const
Expand Down Expand Up @@ -988,12 +856,21 @@ void ForrestTomlinFactorization::dump() const
printf( "\nDumping invQ:\n" );
_invQ.dump();

printf( "*** Done dumping FT factorization ***\n\no" );
printf( "*** Done dumping FT factorization ***\n\n" );
}

const double *ForrestTomlinFactorization::getInvLP() const
void ForrestTomlinFactorization::refactorizeBasis()
{
return _invLP;
for ( unsigned column = 0; column < _m; ++column )
{
const double *basisColumn = _basisColumnOracle->getColumnOfBasis( column );
for ( unsigned row = 0; row < _m; ++row )
_B[row * _m + column] = basisColumn[row];
}

clearFactorization();
initialLUFactorization();
_explicitBasisAvailable = true;
}

//
Expand Down
15 changes: 6 additions & 9 deletions src/basis_factorization/ForrestTomlinFactorization.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@
class ForrestTomlinFactorization : public IBasisFactorization
{
public:
ForrestTomlinFactorization( unsigned m );
ForrestTomlinFactorization( unsigned m, const BasisColumnOracle &basisColumnOracle );
~ForrestTomlinFactorization();

/*
Expand Down Expand Up @@ -88,6 +88,11 @@ class ForrestTomlinFactorization : public IBasisFactorization
*/
void invertBasis( double *result );

/*
Obtain the basis matrix from the oracle and compute a fresh factorization
*/
void refactorizeBasis();

public:
/*
For testing purposes only
Expand All @@ -97,7 +102,6 @@ class ForrestTomlinFactorization : public IBasisFactorization
const EtaMatrix **getU() const;
const List<LPElement *> *getLP() const;
const List<AlmostIdentityMatrix *> *getA() const;
const double *getInvLP() const;

void pushA( const AlmostIdentityMatrix &matrix );
void setQ( const PermutationMatrix &Q );
Expand All @@ -124,13 +128,6 @@ class ForrestTomlinFactorization : public IBasisFactorization
// A1 is the first element of the list
List<AlmostIdentityMatrix *>_A;

/*
A helper matrix to facilitate computing the explicit basis.
This matrix stores:
inv(P1)inv(L1) ... inv(Ps)inv(LS)
*/
double *_invLP;

/*
Indicates whether the explicit basis matrix is available.
*/
Expand Down
16 changes: 15 additions & 1 deletion src/basis_factorization/IBasisFactorization.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,19 @@ class IBasisFactorization
This is the interfact class for a basis factorization.
*/
public:
IBasisFactorization()
/*
A callback for obtaining columns of the basis matrix
*/
class BasisColumnOracle
{
public:
virtual ~BasisColumnOracle() {}
virtual const double *getColumnOfBasis( unsigned column ) const = 0;
};

IBasisFactorization( const BasisColumnOracle &basisColumnOracle )
: _factorizationEnabled( true )
, _basisColumnOracle( &basisColumnOracle )
{
}

Expand Down Expand Up @@ -96,6 +107,9 @@ class IBasisFactorization
disabled.
*/
bool _factorizationEnabled;

protected:
const BasisColumnOracle *_basisColumnOracle;
};

#endif // __IBasisFactorization_h__
Expand Down
Loading

0 comments on commit 6b8125f

Please sign in to comment.