« back to all changes in this revision
Viewing changes to bin/user/poly_patch/math/MMatrix.hh
- browse files at revision 813
- compare with another revision
- download diff
- download tarball
- view history from revision 813
- Committer: Chris Rogers
- Date: 2015-03-06 09:44:06 UTC
- Revision ID: chris.rogers@stfc.ac.uk-20150306094406-872e34ejwevtdh7a
PolyPatch changes
- files added:
- bin/user/poly_patch
- bin/user/poly_patch/math
- bin/user/poly_patch/object
- bin/user/poly_patch/test
- files modified:
- .bzrignore
added
removed
bin/user/poly_patch/math/MMatrix.hh
1
// MAUS WARNING: THIS IS LEGACY CODE.
2
//This file is a part of MAUS
3
//
4
//MAUS is free software: you can redistribute it and/or modify
5
//it under the terms of the GNU General Public License as published by
6
//the Free Software Foundation, either version 3 of the License, or
7
//(at your option) any later version.
8
//
9
//MAUS is distributed in the hope that it will be useful,
10
//but WITHOUT ANY WARRANTY; without even the implied warranty of
11
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12
//GNU General Public License for more details.
13
//
14
//You should have received a copy of the GNU General Public License
15
//along with MAUS in the doc folder. If not, see
16
//<http://www.gnu.org/licenses/>.
17
18
#ifndef _SRC_COMMON_CPP_MATHS_MMATRIX_HH_
19
#define _SRC_COMMON_CPP_MATHS_MMATRIX_HH_
20
21
#include <iostream>
22
#include <vector>
23
24
#include "gsl/gsl_matrix.h"
25
#include "gsl/gsl_blas.h"
26
27
#include "Utils/Exception.hh"
28
typedef MAUS::Exception Squeal;
29
30
#include "math/MVector.hh"
31
32
33
///////////////// MMatrix /////////////////
34
35
/** @class MMatrix C++ wrapper for GSL matrix
36
* MMatrix class handles matrix algebra, maths operators and some
37
* higher level calculation like matrix inversion, eigenvector
38
* analysis etc
39
*
40
* Use template to define two types:\n
41
* (i) MMatrix<double> is a matrix of doubles\n
42
* (ii) MMatrix<m_complex> is a matrix of m_complex\n
43
*
44
* Maths operators and a few others are defined, but maths operators
45
* don't allow operations between types - i.e. you can't multiply
46
* a complex matrix by a double matrix. Instead use interface methods
47
* like real() and complex() to convert between types first
48
*
49
* Nb: CLHEP::HepMatrix can also be used but doesn't have any eigenvalue
50
* analysis. Conversion functions are provided in MMatrixToCLHEP
51
*/
52
template <class Tmplt>
53
class MMatrix
54
{
55
public:
56
/** @brief default constructor makes an empty MMatrix of size (0,0)
57
*/
58
MMatrix();
59
60
/** @brief Copy constructor makes a deep copy of mv
61
*/
62
MMatrix(const MMatrix<Tmplt>& mv );
63
64
/** @brief Construct a matrix and fill with data from memory data_beginning
65
*
66
* @params nrows number of rows
67
* @params ncols number of columns
68
* @params data_beginning pointer to the start of a memory block of size
69
* nrows*ncols with data laid out <row> <row> <row>. Note MMatrix
70
* does not take ownership of memory in data_beginning.
71
*/
72
MMatrix(size_t nrows, size_t ncols, Tmplt* data_beginning );
73
74
/** @brief Construct a matrix and fill with identical data
75
*
76
* @params nrows number of rows
77
* @params ncols number of columns
78
* @params data variable to be copied into all items in the matrix
79
*/
80
MMatrix(size_t nrows, size_t ncols, Tmplt value);
81
82
/** @brief Construct a matrix and fill all fields with 0
83
*
84
* @params nrows number of rows
85
* @params ncols number of columns
86
*/
87
MMatrix(size_t nrows, size_t ncols);
88
89
/** @brief Construct a square matrix filling on and off diagonal values
90
*
91
* @params i number of rows and number of columns
92
* @params diag_value fill values with row == column (i.e. on the diagonal)
93
* with this value
94
* @params off_diag_value fill values with row != column (i.e. off the
95
* diagonal) with this value
96
*/
97
static MMatrix<Tmplt> Diagonal(size_t i, Tmplt diag_value, Tmplt off_diag_value);
98
99
/** @brief destructor
100
*/
101
~MMatrix();
102
103
/** @brief returns number of rows in the matrix
104
*/
105
size_t num_row() const;
106
107
/** @brief returns number of columns in the matrix
108
*/
109
size_t num_col() const;
110
111
/** @brief returns sum of terms with row == column, even if matrix is not
112
* square
113
*/
114
Tmplt trace() const;
115
116
/** @brief returns matrix determinant; throws an exception if matrix is not
117
* square
118
*/
119
Tmplt determinant() const;
120
121
/** @brief returns matrix inverse leaving this matrix unchanged
122
*/
123
MMatrix<Tmplt> inverse() const;
124
125
/** @brief turns this matrix into its inverse
126
*/
127
void invert();
128
129
/** @brief returns matrix transpose T (such that M(i,j) = T(j,i))
130
*/
131
MMatrix<Tmplt> T() const;
132
133
/** @brief returns a vector of eigenvalues. Throws an exception if either this
134
* matrix is not square or the eigenvalues could not be found (e.g.
135
* singular matrix or whatever).
136
*/
137
MVector<m_complex> eigenvalues() const; //return vector of eigenvalues
138
std::pair<MVector<m_complex>, MMatrix<m_complex> >
139
eigenvectors() const; //return pair of eigenvalues, eigenvectors (access using my_pair.first or my_pair.second)
140
141
//return a submatrix, with data from min_row to max_row and min_row to max_row inclusive; indexing starts at 1
142
MMatrix<Tmplt> sub(size_t min_row, size_t max_row, size_t min_col, size_t max_col) const;
143
//create a column vector from the column^th column
144
MVector<Tmplt> get_mvector(size_t column) const;
145
//return a reference to the datum; indexing starts at 1, goes to num_row (inclusive)
146
const Tmplt& operator()(size_t row, size_t column) const;
147
Tmplt& operator()(size_t row, size_t column);
148
149
MMatrix<Tmplt>& operator= (const MMatrix<Tmplt>& mm);
150
151
//TODO - implement iterator
152
//class iterator
153
//{
154
//}
155
156
157
friend MMatrix<m_complex>& operator *=(MMatrix<m_complex>& m, m_complex c);
158
friend MMatrix<double>& operator *=(MMatrix<double>& m, double d);
159
friend MMatrix<m_complex>& operator *=(MMatrix<m_complex>& m1, MMatrix<m_complex> m2);
160
friend MMatrix<double>& operator *=(MMatrix<double>& m1, MMatrix<double> m2);
161
friend MVector<m_complex> operator * (MMatrix<m_complex> m, MVector<m_complex> v);
162
friend MVector<double> operator * (MMatrix<double> m, MVector<double> v);
163
friend MMatrix<m_complex>& operator +=(MMatrix<m_complex>& m1, const MMatrix<m_complex>& m2);
164
friend MMatrix<double>& operator +=(MMatrix<double>& m1, const MMatrix<double>& m2);
165
template <class Tmplt2> friend MMatrix<Tmplt2> operator + (MMatrix<Tmplt2> m1, const MMatrix<Tmplt2> m2);
166
167
friend const gsl_matrix* MMatrix_to_gsl(const MMatrix<double>& m);
168
friend const gsl_matrix_complex* MMatrix_to_gsl(const MMatrix<gsl_complex>& m);
169
170
friend class MMatrix<double>; //To do the eigenvector problem, MMatrix<double> needs to see MMatrix<complex>'s _matrix
171
172
173
private:
174
//build the matrix with size i,j, data not initialised
175
void build_matrix(size_t i, size_t j);
176
//build the matrix with size i,j, data initialised to data in temp
177
void build_matrix(size_t i, size_t j, Tmplt* temp);
178
//delete the matrix and set it to null
179
void delete_matrix();
180
//return the matrix overloaded to the correct type or throw
181
//if _matrix == NULL
182
static gsl_matrix* get_matrix(const MMatrix<double>& m);
183
static gsl_matrix_complex* get_matrix(const MMatrix<m_complex>& m);
184
//gsl object
185
void* _matrix;
186
};
187
188
//Operators do what you would expect - i.e. normal mathematical functions
189
MMatrix<m_complex>& operator *=(MMatrix<m_complex>& m, m_complex c);
190
MMatrix<double>& operator *=(MMatrix<double>& m, double d);
191
MMatrix<m_complex>& operator *=(MMatrix<m_complex>& m1, MMatrix<m_complex> m2); //M1 *= M2 returns M1 = M1*M2
192
MMatrix<double>& operator *=(MMatrix<double>& m1, MMatrix<double> m2); //M1 *= M2 returns M1 = M1*M2
193
MVector<m_complex> operator * (MMatrix<m_complex> m, MVector<m_complex> v);
194
MVector<double> operator * (MMatrix<double> m, MVector<double> v);
195
MMatrix<m_complex>& operator +=(MMatrix<m_complex>& m1, const MMatrix<m_complex>& m2);
196
MMatrix<double>& operator +=(MMatrix<double>& m1, const MMatrix<double>& m2);
197
198
template <class Tmplt> MMatrix<Tmplt>& operator -=(MMatrix<double>& m1, MMatrix<double> m2);
199
200
template <class Tmplt> MMatrix<Tmplt> operator*(MMatrix<Tmplt>, MMatrix<Tmplt>);
201
template <class Tmplt> MMatrix<Tmplt> operator*(MMatrix<Tmplt>, MVector<Tmplt>);
202
template <class Tmplt> MMatrix<Tmplt> operator*(MMatrix<Tmplt>, Tmplt);
203
template <class Tmplt> MMatrix<Tmplt> operator/(MMatrix<Tmplt>, Tmplt);
204
205
template <class Tmplt> MMatrix<Tmplt> operator+(MMatrix<Tmplt>, MMatrix<Tmplt>);
206
template <class Tmplt> MMatrix<Tmplt> operator-(MMatrix<Tmplt>, MMatrix<Tmplt>);
207
template <class Tmplt> MMatrix<Tmplt> operator-(MMatrix<Tmplt>); //unary minus returns matrix*(-1)
208
209
template <class Tmplt> bool operator==(const MMatrix<Tmplt>& c1, const MMatrix<Tmplt>& c2);
210
template <class Tmplt> bool operator!=(const MMatrix<Tmplt>& c1, const MMatrix<Tmplt>& c2);
211
212
//format is
213
// num_row num_col
214
// m11 m12 m13 ...
215
// m21 m22 m23 ...
216
// ...
217
template <class Tmplt> std::ostream& operator<<(std::ostream& out, MMatrix<Tmplt> mat);
218
template <class Tmplt> std::istream& operator>>(std::istream& in, MMatrix<Tmplt>& mat);
219
220
//return matrix of doubles filled with either real or imaginary part of m
221
MMatrix<double> re(MMatrix<m_complex> m);
222
MMatrix<double> im(MMatrix<m_complex> m);
223
//return matrix of m_complex filled with real and imaginary parts
224
MMatrix<m_complex> complex(MMatrix<double> real);
225
MMatrix<m_complex> complex(MMatrix<double> real, MMatrix<double> imaginary);
226
227
//return pointer to gsl_matrix objects that store matrix data in m
228
const gsl_matrix* MMatrix_to_gsl(const MMatrix<double>& m);
229
const gsl_matrix_complex* MMatrix_to_gsl(const MMatrix<gsl_complex>& m);
230
231
//////////////////////////// MMatrix declaration end ///////////////
232
233
234
235
//////////////////////////// MMatrix Inlined Functions //////////////
236
237
template <>
238
inline size_t MMatrix<double>::num_row() const
239
{
240
if(_matrix) return ((gsl_matrix*)_matrix)->size1;
241
return 0;
242
}
243
template <>
244
inline size_t MMatrix<m_complex>::num_row() const
245
{
246
if(_matrix) return ((gsl_matrix_complex*)_matrix)->size1;
247
return 0;
248
}
249
template <>
250
inline size_t MMatrix<double>::num_col() const
251
{
252
if(_matrix) return ((gsl_matrix*)_matrix)->size2;
253
return 0;
254
}
255
256
template <>
257
inline size_t MMatrix<m_complex>::num_col() const
258
{
259
if(_matrix) return ((gsl_matrix_complex*)_matrix)->size2;
260
return 0;
261
}
262
263
MMatrix<m_complex> inline & operator *=(MMatrix<m_complex>& m, m_complex c)
264
{gsl_matrix_complex_scale((gsl_matrix_complex*)m._matrix, c); return m;}
265
266
MMatrix<double> inline & operator *=(MMatrix<double>& m, double d)
267
{gsl_matrix_scale((gsl_matrix*)m._matrix, d); return m;}
268
269
MVector<m_complex> inline operator * (MMatrix<m_complex> m, MVector<m_complex> v)
270
{
271
MVector<m_complex> v0(m.num_row());
272
gsl_blas_zgemv(CblasNoTrans, m_complex_build(1.), (gsl_matrix_complex*)m._matrix, (gsl_vector_complex*)v._vector, m_complex_build(0.), (gsl_vector_complex*)v0._vector);
273
return v0;
274
}
275
276
MVector<double> inline operator * (MMatrix<double> m, MVector<double> v)
277
{
278
MVector<double> v0(m.num_row());
279
gsl_blas_dgemv(CblasNoTrans, 1., (gsl_matrix*)m._matrix, (gsl_vector*)v._vector, 0., (gsl_vector*)v0._vector);
280
return v0;
281
}
282
283
MMatrix<m_complex> inline & operator +=(MMatrix<m_complex>& m1, const MMatrix<m_complex>& m2)
284
{gsl_matrix_complex_add((gsl_matrix_complex*)m1._matrix, (gsl_matrix_complex*)m2._matrix); return m1;}
285
286
MMatrix<double> inline & operator +=(MMatrix<double>& m1, const MMatrix<double>& m2)
287
{gsl_matrix_add((gsl_matrix*)m1._matrix, (gsl_matrix*)m2._matrix); return m1;}
288
289
template <class Tmplt> MMatrix<Tmplt> inline operator - (MMatrix<Tmplt> m1)
290
{for(size_t i=1; i<=m1.num_row(); i++) for(size_t j=1; j<=m1.num_col(); j++) m1(i,j) = -1.*m1(i,j); return m1; }
291
292
template <class Tmplt> MMatrix<Tmplt> inline operator*(MMatrix<Tmplt> m1, MMatrix<Tmplt> m2)
293
{return m1*=m2;}
294
295
template <class Tmplt> MMatrix<Tmplt> inline operator*(MMatrix<Tmplt> m, MVector<Tmplt> v)
296
{return m*=v;}
297
298
template <class Tmplt> MMatrix<Tmplt> inline operator*(MMatrix<Tmplt> m, Tmplt t)
299
{return m*=t;}
300
301
template <class Tmplt> MMatrix<Tmplt> inline & operator/=(MMatrix<Tmplt>& m, Tmplt t)
302
{return m*=1./t;}
303
template <class Tmplt> MMatrix<Tmplt> inline operator/(MMatrix<Tmplt> m, Tmplt t)
304
{return m/=t;}
305
306
template <class Tmplt> MMatrix<Tmplt> inline operator+(MMatrix<Tmplt> m1, MMatrix<Tmplt> m2)
307
{ MMatrix<Tmplt> m3 = m1+=m2; return m3; }
308
309
template <class Tmplt> MMatrix<Tmplt> inline & operator-=(MMatrix<Tmplt>& m1, MMatrix<Tmplt> m2)
310
{return m1 += -m2;}
311
template <class Tmplt> MMatrix<Tmplt> inline operator-(MMatrix<Tmplt> m1, MMatrix<Tmplt> m2)
312
{return m1 -= m2;}
313
314
315
316
template <>
317
const double inline & MMatrix<double>::operator()(size_t i, size_t j) const
318
{ return *gsl_matrix_ptr((gsl_matrix*)_matrix, i-1, j-1); }
319
template <>
320
const m_complex inline & MMatrix<m_complex>::operator()(size_t i, size_t j) const
321
{ return *gsl_matrix_complex_ptr((gsl_matrix_complex*)_matrix, i-1, j-1); }
322
323
template <>
324
double inline & MMatrix<double>::operator()(size_t i, size_t j)
325
{ return *gsl_matrix_ptr((gsl_matrix*)_matrix, i-1, j-1); }
326
template <>
327
m_complex inline & MMatrix<m_complex>::operator()(size_t i, size_t j)
328
{ return *gsl_matrix_complex_ptr((gsl_matrix_complex*)_matrix, i-1, j-1); }
329
330
template <class Tmplt>
331
bool operator==(const MMatrix<Tmplt>& lhs, const MMatrix<Tmplt>& rhs)
332
{
333
if( lhs.num_row() != rhs.num_row() || lhs.num_col() != rhs.num_col() ) return false;
334
for(size_t i=1; i<=lhs.num_row(); i++)
335
for(size_t j=1; j<=lhs.num_col(); j++)
336
if(lhs(i,j) != rhs(i,j)) return false;
337
return true;
338
}
339
340
template <class Tmplt>
341
bool inline operator!=(const MMatrix<Tmplt>& lhs, const MMatrix<Tmplt>& rhs) {return ! (lhs == rhs);}
342
343
//iostream
344
template <class Tmplt> std::ostream& operator<<(std::ostream& out, MMatrix<Tmplt>);
345
//template <class Tmplt> std::istream& operator>>(std::istream& in, MMatrix<Tmplt>);
346
347
template <class Tmplt>
348
gsl_matrix inline* MMatrix<Tmplt>::get_matrix(const MMatrix<double>& m)
349
{
350
if(m._matrix == NULL)
351
throw(Squeal(Squeal::recoverable, "Attempt to access uninitialised matrix", "MMatrix::get_matrix"));
352
return (gsl_matrix*)m._matrix;
353
}
354
355
template <class Tmplt>
356
gsl_matrix_complex inline* MMatrix<Tmplt>::get_matrix(const MMatrix<m_complex>& m)
357
{
358
if(m._matrix == NULL)
359
throw(Squeal(Squeal::recoverable, "Attempt to access uninitialised matrix", "MMatrix::get_matrix"));
360
return (gsl_matrix_complex*)m._matrix;
361
362
}
363
364
////////////////////////// MMatrix inlined functions end //////////////////////
365
366
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1