Changeset View
Changeset View
Standalone View
Standalone View
EvtGenBase/EvtMatrix.hh
Show All 22 Lines | |||||
#include <cmath> | #include <cmath> | ||||
#include <sstream> | #include <sstream> | ||||
#include <vector> | #include <vector> | ||||
template <class T> | template <class T> | ||||
class EvtMatrix { | class EvtMatrix { | ||||
private: | private: | ||||
T** _mat; | T** m_mat; | ||||
int _range; | int m_range; | ||||
public: | public: | ||||
EvtMatrix() : _range( 0 ){}; | EvtMatrix() : m_range( 0 ){}; | ||||
~EvtMatrix(); | ~EvtMatrix(); | ||||
inline void setRange( int range ); | inline void setRange( int range ); | ||||
T& operator()( int row, int col ) { return _mat[row][col]; } | T& operator()( int row, int col ) { return m_mat[row][col]; } | ||||
T* operator[]( int row ) { return _mat[row]; } | T* operator[]( int row ) { return m_mat[row]; } | ||||
T det(); | T det(); | ||||
EvtMatrix* min( int row, int col ); | EvtMatrix* min( int row, int col ); | ||||
EvtMatrix* inverse(); | EvtMatrix* inverse(); | ||||
std::string dump(); | std::string dump(); | ||||
template <class M> | template <class M> | ||||
friend EvtMatrix<M>* operator*( const EvtMatrix<M>& left, | friend EvtMatrix<M>* operator*( const EvtMatrix<M>& left, | ||||
const EvtMatrix<M>& right ); | const EvtMatrix<M>& right ); | ||||
}; | }; | ||||
template <class T> | template <class T> | ||||
inline void EvtMatrix<T>::setRange( int range ) | inline void EvtMatrix<T>::setRange( int range ) | ||||
{ | { | ||||
// If the range is changed, delete any previous matrix stored | // If the range is changed, delete any previous matrix stored | ||||
// and allocate elements with the newly specified range. | // and allocate elements with the newly specified range. | ||||
if ( _range != range ) { | if ( m_range != range ) { | ||||
if ( _range ) { | if ( m_range ) { | ||||
for ( int row = 0; row < _range; row++ ) | for ( int row = 0; row < m_range; row++ ) | ||||
delete[] _mat[row]; | delete[] m_mat[row]; | ||||
delete[] _mat; | delete[] m_mat; | ||||
} | } | ||||
_mat = new T*[range]; | m_mat = new T*[range]; | ||||
for ( int row = 0; row < range; row++ ) | for ( int row = 0; row < range; row++ ) | ||||
_mat[row] = new T[range]; | m_mat[row] = new T[range]; | ||||
// Set the new range. | // Set the new range. | ||||
_range = range; | m_range = range; | ||||
} | } | ||||
// Since user is willing to change the range, reset the matrix elements. | // Since user is willing to change the range, reset the matrix elements. | ||||
for ( int row = 0; row < _range; row++ ) | for ( int row = 0; row < m_range; row++ ) | ||||
for ( int col = 0; col < _range; col++ ) | for ( int col = 0; col < m_range; col++ ) | ||||
_mat[row][col] = 0.; | m_mat[row][col] = 0.; | ||||
} | } | ||||
template <class T> | template <class T> | ||||
EvtMatrix<T>::~EvtMatrix() | EvtMatrix<T>::~EvtMatrix() | ||||
{ | { | ||||
for ( int row = 0; row < _range; row++ ) | for ( int row = 0; row < m_range; row++ ) | ||||
delete[] _mat[row]; | delete[] m_mat[row]; | ||||
delete[] _mat; | delete[] m_mat; | ||||
} | } | ||||
template <class T> | template <class T> | ||||
std::string EvtMatrix<T>::dump() | std::string EvtMatrix<T>::dump() | ||||
{ | { | ||||
std::ostringstream str; | std::ostringstream str; | ||||
for ( int row = 0; row < _range; row++ ) { | for ( int row = 0; row < m_range; row++ ) { | ||||
str << "|"; | str << "|"; | ||||
for ( int col = 0; col < _range; col++ ) | for ( int col = 0; col < m_range; col++ ) | ||||
str << "\t" << _mat[row][col]; | str << "\t" << m_mat[row][col]; | ||||
str << "\t|" << std::endl; | str << "\t|" << std::endl; | ||||
} | } | ||||
return str.str(); | return str.str(); | ||||
} | } | ||||
template <class T> | template <class T> | ||||
T EvtMatrix<T>::det() | T EvtMatrix<T>::det() | ||||
{ | { | ||||
if ( _range == 1 ) | if ( m_range == 1 ) | ||||
return _mat[0][0]; | return m_mat[0][0]; | ||||
// There's no need to define the range 2 determinant manually, but it may | // There's no need to define the range 2 determinant manually, but it may | ||||
// speed up the calculation. | // speed up the calculation. | ||||
if ( _range == 2 ) | if ( m_range == 2 ) | ||||
return _mat[0][0] * _mat[1][1] - _mat[0][1] * _mat[1][0]; | return m_mat[0][0] * m_mat[1][1] - m_mat[0][1] * m_mat[1][0]; | ||||
T sum = 0.; | T sum = 0.; | ||||
for ( int col = 0; col < _range; col++ ) { | for ( int col = 0; col < m_range; col++ ) { | ||||
EvtMatrix<T>* minor = min( 0, col ); | EvtMatrix<T>* minor = min( 0, col ); | ||||
sum += std::pow( -1., col ) * _mat[0][col] * minor->det(); | sum += std::pow( -1., col ) * m_mat[0][col] * minor->det(); | ||||
delete minor; | delete minor; | ||||
} | } | ||||
return sum; | return sum; | ||||
} | } | ||||
// Returns the minor at (i, j). | // Returns the minor at (i, j). | ||||
template <class T> | template <class T> | ||||
EvtMatrix<T>* EvtMatrix<T>::min( int row, int col ) | EvtMatrix<T>* EvtMatrix<T>::min( int row, int col ) | ||||
{ | { | ||||
EvtMatrix<T>* minor = new EvtMatrix<T>(); | EvtMatrix<T>* minor = new EvtMatrix<T>(); | ||||
minor->setRange( _range - 1 ); | minor->setRange( m_range - 1 ); | ||||
int minIndex = 0; | int minIndex = 0; | ||||
for ( int r = 0; r < _range; r++ ) | for ( int r = 0; r < m_range; r++ ) | ||||
for ( int c = 0; c < _range; c++ ) | for ( int c = 0; c < m_range; c++ ) | ||||
if ( ( r != row ) && ( c != col ) ) { | if ( ( r != row ) && ( c != col ) ) { | ||||
( *minor )( minIndex / ( _range - 1 ), | ( *minor )( minIndex / ( m_range - 1 ), | ||||
minIndex % ( _range - 1 ) ) = _mat[r][c]; | minIndex % ( m_range - 1 ) ) = m_mat[r][c]; | ||||
minIndex++; | minIndex++; | ||||
} | } | ||||
return minor; | return minor; | ||||
} | } | ||||
template <class T> | template <class T> | ||||
EvtMatrix<T>* EvtMatrix<T>::inverse() | EvtMatrix<T>* EvtMatrix<T>::inverse() | ||||
{ | { | ||||
EvtMatrix<T>* inv = new EvtMatrix<T>(); | EvtMatrix<T>* inv = new EvtMatrix<T>(); | ||||
inv->setRange( _range ); | inv->setRange( m_range ); | ||||
if ( det() == 0 ) { | if ( det() == 0 ) { | ||||
std::cerr << "This matrix has a null determinant and cannot be inverted. Returning zero matrix." | std::cerr << "This matrix has a null determinant and cannot be inverted. Returning zero matrix." | ||||
<< std::endl; | << std::endl; | ||||
for ( int row = 0; row < _range; row++ ) | for ( int row = 0; row < m_range; row++ ) | ||||
for ( int col = 0; col < _range; col++ ) | for ( int col = 0; col < m_range; col++ ) | ||||
( *inv )( row, col ) = 0.; | ( *inv )( row, col ) = 0.; | ||||
return inv; | return inv; | ||||
} | } | ||||
T determinant = det(); | T determinant = det(); | ||||
for ( int row = 0; row < _range; row++ ) | for ( int row = 0; row < m_range; row++ ) | ||||
for ( int col = 0; col < _range; col++ ) { | for ( int col = 0; col < m_range; col++ ) { | ||||
EvtMatrix<T>* minor = min( row, col ); | EvtMatrix<T>* minor = min( row, col ); | ||||
inv->_mat[col][row] = std::pow( -1., row + col ) * minor->det() / | inv->m_mat[col][row] = std::pow( -1., row + col ) * minor->det() / | ||||
determinant; | determinant; | ||||
delete minor; | delete minor; | ||||
} | } | ||||
return inv; | return inv; | ||||
} | } | ||||
template <class T> | template <class T> | ||||
EvtMatrix<T>* operator*( const EvtMatrix<T>& left, const EvtMatrix<T>& right ) | EvtMatrix<T>* operator*( const EvtMatrix<T>& left, const EvtMatrix<T>& right ) | ||||
{ | { | ||||
// Chech that the matrices have the correct range. | // Chech that the matrices have the correct range. | ||||
if ( left._range != right._range ) { | if ( left.m_range != right.m_range ) { | ||||
std::cerr << "These matrices cannot be multiplied." << std::endl; | std::cerr << "These matrices cannot be multiplied." << std::endl; | ||||
return new EvtMatrix<T>(); | return new EvtMatrix<T>(); | ||||
} | } | ||||
EvtMatrix<T>* mat = new EvtMatrix<T>(); | EvtMatrix<T>* mat = new EvtMatrix<T>(); | ||||
mat->setRange( left._range ); | mat->setRange( left.m_range ); | ||||
// Initialize the elements of the matrix. | // Initialize the elements of the matrix. | ||||
for ( int row = 0; row < left._range; row++ ) | for ( int row = 0; row < left.m_range; row++ ) | ||||
for ( int col = 0; col < right._range; col++ ) | for ( int col = 0; col < right.m_range; col++ ) | ||||
( *mat )[row][col] = 0; | ( *mat )[row][col] = 0; | ||||
for ( int row = 0; row < left._range; row++ ) | for ( int row = 0; row < left.m_range; row++ ) | ||||
for ( int col = 0; col < right._range; col++ ) | for ( int col = 0; col < right.m_range; col++ ) | ||||
for ( int line = 0; line < right._range; line++ ) | for ( int line = 0; line < right.m_range; line++ ) | ||||
( *mat )[row][col] += left._mat[row][line] * | ( *mat )[row][col] += left.m_mat[row][line] * | ||||
right._mat[line][col]; | right.m_mat[line][col]; | ||||
return mat; | return mat; | ||||
} | } | ||||
#endif | #endif |