Line data Source code
1 : /*****************************************************************************
2 : * Project: BaBar detector at the SLAC PEP-II B-factory
3 : * Package: EvtGenBase
4 : * File: $Id: EvtMatrix.hh,v 1.1 2009-03-16 16:50:49 robbep Exp $
5 : *
6 : * Description:
7 : * Class to make simple computations with matrices: assignment, product,
8 : * determinant, inverse... Still many functions could be implemented.
9 : *
10 : * Modification history:
11 : * Jordi Garra Ticó 2008/07/03 File created
12 : *****************************************************************************/
13 :
14 :
15 : #ifndef __EVT_MATRIX_HH__
16 : #define __EVT_MATRIX_HH__
17 :
18 : #include <vector>
19 : #include <sstream>
20 :
21 :
22 : template <class T> class EvtMatrix
23 : {
24 : private:
25 : T** _mat;
26 : int _range;
27 : public:
28 0 : EvtMatrix() : _range( 0 ) {};
29 : ~EvtMatrix();
30 : inline void setRange( int range );
31 :
32 0 : T& operator()( int row, int col ) { return _mat[ row ][ col ]; }
33 0 : T* operator[]( int row ) { return _mat[ row ]; }
34 : T det();
35 : EvtMatrix* min( int row, int col );
36 : EvtMatrix* inverse();
37 : std::string dump();
38 :
39 : template <class M> friend EvtMatrix< M >* operator*( const EvtMatrix< M >& left, const EvtMatrix< M >& right );
40 : };
41 :
42 :
43 :
44 :
45 : template <class T> inline void EvtMatrix< T >::setRange( int range )
46 : {
47 : // If the range is changed, delete any previous matrix stored
48 : // and allocate elements with the newly specified range.
49 0 : if ( _range != range )
50 : {
51 0 : if ( _range )
52 : {
53 0 : for ( int row = 0; row < _range; row++ )
54 0 : delete[] _mat[ row ];
55 0 : delete[] _mat;
56 : }
57 :
58 0 : _mat = new T*[ range ];
59 0 : for ( int row = 0; row < range; row++ )
60 0 : _mat[ row ] = new T[ range ];
61 :
62 : // Set the new range.
63 0 : _range = range;
64 0 : }
65 :
66 : // Since user is willing to change the range, reset the matrix elements.
67 0 : for ( int row = 0; row < _range; row++ )
68 0 : for ( int col = 0; col < _range; col++ )
69 0 : _mat[ row ][ col ] = 0.;
70 0 : }
71 :
72 :
73 : template <class T> EvtMatrix< T >::~EvtMatrix()
74 0 : {
75 0 : for( int row = 0; row < _range; row++ )
76 0 : delete[] _mat[ row ];
77 0 : delete[] _mat;
78 0 : }
79 :
80 :
81 : template <class T> std::string EvtMatrix< T >::dump()
82 : {
83 : std::ostringstream str;
84 :
85 : for ( int row = 0; row < _range; row++ )
86 : {
87 : str << "|";
88 : for ( int col = 0; col < _range; col++ )
89 : str << "\t" << _mat[ row ][ col ];
90 : str << "\t|" << std::endl;
91 : }
92 :
93 : return str.str();
94 : }
95 :
96 :
97 : template <class T> T EvtMatrix< T >::det()
98 : {
99 0 : if ( _range == 1 )
100 0 : return _mat[ 0 ][ 0 ];
101 :
102 : // There's no need to define the range 2 determinant manually, but it may
103 : // speed up the calculation.
104 0 : if ( _range == 2 )
105 0 : return _mat[ 0 ][ 0 ] * _mat[ 1 ][ 1 ] - _mat[ 0 ][ 1 ] * _mat[ 1 ][ 0 ];
106 :
107 0 : T sum = 0.;
108 :
109 0 : for ( int col = 0; col < _range; col++ )
110 : {
111 0 : EvtMatrix< T >* minor = min( 0, col );
112 0 : sum += std::pow( -1., col ) * _mat[ 0 ][ col ] * minor->det();
113 0 : delete minor;
114 : }
115 :
116 0 : return sum;
117 0 : }
118 :
119 :
120 : // Returns the minor at (i, j).
121 : template <class T> EvtMatrix< T >* EvtMatrix< T >::min( int row, int col )
122 : {
123 0 : EvtMatrix< T >* minor = new EvtMatrix< T >();
124 0 : minor->setRange( _range - 1 );
125 :
126 : int minIndex = 0;
127 :
128 0 : for ( int r = 0; r < _range; r++ )
129 0 : for ( int c = 0; c < _range; c++ )
130 0 : if ( ( r != row ) && ( c != col ) )
131 : {
132 0 : (*minor)( minIndex / ( _range - 1 ), minIndex % ( _range - 1 ) ) = _mat[ r ][ c ];
133 0 : minIndex++;
134 0 : }
135 :
136 0 : return minor;
137 0 : }
138 :
139 :
140 : template <class T> EvtMatrix< T >* EvtMatrix< T >::inverse()
141 : {
142 0 : EvtMatrix< T >* inv = new EvtMatrix< T >();
143 0 : inv->setRange( _range );
144 :
145 0 : if ( det() == 0 )
146 : {
147 0 : std::cerr << "This matrix has a null determinant and cannot be inverted. Returning zero matrix." << std::endl;
148 0 : for ( int row = 0; row < _range; row++ )
149 0 : for ( int col = 0; col < _range; col++ )
150 0 : (*inv)( row, col ) = 0.;
151 0 : return inv;
152 : }
153 :
154 0 : T determinant = det();
155 :
156 0 : for ( int row = 0; row < _range; row++ )
157 0 : for ( int col = 0; col < _range; col++ )
158 : {
159 0 : EvtMatrix< T >* minor = min( row, col );
160 0 : inv->_mat[col][row] = std::pow( -1., row + col ) * minor->det() / determinant;
161 0 : delete minor;
162 : }
163 :
164 : return inv;
165 0 : }
166 :
167 :
168 : template <class T>
169 : EvtMatrix< T >* operator*( const EvtMatrix< T >& left, const EvtMatrix< T >& right )
170 : {
171 : // Chech that the matrices have the correct range.
172 : if ( left._range != right._range )
173 : {
174 : std::cerr << "These matrices cannot be multiplied." << std::endl;
175 : return new EvtMatrix< T >();
176 : }
177 :
178 : EvtMatrix< T >* mat = new EvtMatrix< T >();
179 : mat->setRange( left._range );
180 :
181 : // Initialize the elements of the matrix.
182 : for ( int row = 0; row < left._range; row++ )
183 : for ( int col = 0; col < right._range; col++ )
184 : (*mat)[ row ][ col ] = 0;
185 :
186 : for ( int row = 0; row < left._range; row++ )
187 : for ( int col = 0; col < right._range; col++ )
188 : for ( int line = 0; line < right._range; line++ )
189 : (*mat)[ row ][ col ] += left._mat[ row ][ line ] * right._mat[ line ][ col ];
190 :
191 : return mat;
192 : }
193 :
194 :
195 : #endif
|
