00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #ifndef __EVT_MATRIX_HH__
00016 #define __EVT_MATRIX_HH__
00017
00018 #include <vector>
00019 #include <sstream>
00020
00021
00022 template <class T> class EvtMatrix
00023 {
00024 private:
00025 T** _mat;
00026 int _range;
00027 public:
00028 EvtMatrix() : _range( 0 ) {};
00029 ~EvtMatrix();
00030 inline void setRange( int range );
00031
00032 T& operator()( int row, int col ) { return _mat[ row ][ col ]; }
00033 T* operator[]( int row ) { return _mat[ row ]; }
00034 T det();
00035 EvtMatrix* min( int row, int col );
00036 EvtMatrix* inverse();
00037 std::string dump();
00038
00039 template <class M> friend EvtMatrix< M >* operator*( const EvtMatrix< M >& left, const EvtMatrix< M >& right );
00040 };
00041
00042
00043
00044
00045 template <class T> inline void EvtMatrix< T >::setRange( int range )
00046 {
00047
00048
00049 if ( _range != range )
00050 {
00051 if ( _range )
00052 {
00053 for ( int row = 0; row < _range; row++ )
00054 delete[] _mat[ row ];
00055 delete[] _mat;
00056 }
00057
00058 _mat = new T*[ range ];
00059 for ( int row = 0; row < range; row++ )
00060 _mat[ row ] = new T[ range ];
00061
00062
00063 _range = range;
00064 }
00065
00066
00067 for ( int row = 0; row < _range; row++ )
00068 for ( int col = 0; col < _range; col++ )
00069 _mat[ row ][ col ] = 0.;
00070 }
00071
00072
00073 template <class T> EvtMatrix< T >::~EvtMatrix()
00074 {
00075 for( int row = 0; row < _range; row++ )
00076 delete[] _mat[ row ];
00077 delete[] _mat;
00078 }
00079
00080
00081 template <class T> std::string EvtMatrix< T >::dump()
00082 {
00083 std::ostringstream str;
00084
00085 for ( int row = 0; row < _range; row++ )
00086 {
00087 str << "|";
00088 for ( int col = 0; col < _range; col++ )
00089 str << "\t" << _mat[ row ][ col ];
00090 str << "\t|" << std::endl;
00091 }
00092
00093 return str.str();
00094 }
00095
00096
00097 template <class T> T EvtMatrix< T >::det()
00098 {
00099 if ( _range == 1 )
00100 return _mat[ 0 ][ 0 ];
00101
00102
00103
00104 if ( _range == 2 )
00105 return _mat[ 0 ][ 0 ] * _mat[ 1 ][ 1 ] - _mat[ 0 ][ 1 ] * _mat[ 1 ][ 0 ];
00106
00107 T sum = 0.;
00108
00109 for ( int col = 0; col < _range; col++ )
00110 {
00111 EvtMatrix< T >* minor = min( 0, col );
00112 sum += std::pow( -1., col ) * _mat[ 0 ][ col ] * minor->det();
00113 delete minor;
00114 }
00115
00116 return sum;
00117 }
00118
00119
00120
00121 template <class T> EvtMatrix< T >* EvtMatrix< T >::min( int row, int col )
00122 {
00123 EvtMatrix< T >* minor = new EvtMatrix< T >();
00124 minor->setRange( _range - 1 );
00125
00126 int minIndex = 0;
00127
00128 for ( int r = 0; r < _range; r++ )
00129 for ( int c = 0; c < _range; c++ )
00130 if ( ( r != row ) && ( c != col ) )
00131 {
00132 (*minor)( minIndex / ( _range - 1 ), minIndex % ( _range - 1 ) ) = _mat[ r ][ c ];
00133 minIndex++;
00134 }
00135
00136 return minor;
00137 }
00138
00139
00140 template <class T> EvtMatrix< T >* EvtMatrix< T >::inverse()
00141 {
00142 EvtMatrix< T >* inv = new EvtMatrix< T >();
00143 inv->setRange( _range );
00144
00145 if ( det() == 0 )
00146 {
00147 std::cerr << "This matrix has a null determinant and cannot be inverted. Returning zero matrix." << std::endl;
00148 for ( int row = 0; row < _range; row++ )
00149 for ( int col = 0; col < _range; col++ )
00150 (*inv)( row, col ) = 0.;
00151 return inv;
00152 }
00153
00154 T determinant = det();
00155
00156 for ( int row = 0; row < _range; row++ )
00157 for ( int col = 0; col < _range; col++ )
00158 {
00159 EvtMatrix< T >* minor = min( row, col );
00160 inv->_mat[col][row] = std::pow( -1., row + col ) * minor->det() / determinant;
00161 delete minor;
00162 }
00163
00164 return inv;
00165 }
00166
00167
00168 template <class T>
00169 EvtMatrix< T >* operator*( const EvtMatrix< T >& left, const EvtMatrix< T >& right )
00170 {
00171
00172 if ( left._range != right._range )
00173 {
00174 std::cerr << "These matrices cannot be multiplied." << std::endl;
00175 return new EvtMatrix< T >();
00176 }
00177
00178 EvtMatrix< T >* mat = new EvtMatrix< T >();
00179 mat->setRange( left._range );
00180
00181
00182 for ( int row = 0; row < left._range; row++ )
00183 for ( int col = 0; col < right._range; col++ )
00184 (*mat)[ row ][ col ] = 0;
00185
00186 for ( int row = 0; row < left._range; row++ )
00187 for ( int col = 0; col < right._range; col++ )
00188 for ( int line = 0; line < right._range; line++ )
00189 (*mat)[ row ][ col ] += left._mat[ row ][ line ] * right._mat[ line ][ col ];
00190
00191 return mat;
00192 }
00193
00194
00195 #endif