]> git.tdb.fi Git - libs/math.git/blob - source/linal/dynamicmatrix.h
Rename the low-level matrix inversion function to gauss_jordan
[libs/math.git] / source / linal / dynamicmatrix.h
1 #ifndef MSP_LINAL_DYNAMICMATRIX_H_
2 #define MSP_LINAL_DYNAMICMATRIX_H_
3
4 #include <algorithm>
5 #include "dynamicvector.h"
6 #include "matrixops.h"
7
8 namespace Msp {
9 namespace LinAl {
10
11 /**
12 A general mathematical matrix.  Size is specified at runtime.  May be slower
13 than a fixed-size matrix.  There are no compile-time diagnostics for mismatched
14 matrix sizes.
15 */
16 template<typename T>
17 class DynamicMatrix
18 {
19 public:
20         typedef T ElementType;
21
22 private:
23         unsigned rows_;
24         unsigned columns_;
25         T *data;
26
27 public:
28         DynamicMatrix(unsigned, unsigned);
29         DynamicMatrix(unsigned, unsigned, const T *);
30         DynamicMatrix(const DynamicMatrix &);
31         DynamicMatrix &operator=(const DynamicMatrix &);
32         ~DynamicMatrix();
33
34         unsigned rows() const { return rows_; }
35         unsigned columns() const { return columns_; }
36
37         T &element(unsigned, unsigned);
38         const T &element(unsigned, unsigned) const;
39         T &operator()(unsigned i, unsigned j) { return element(i, j); }
40         const T &operator()(unsigned i, unsigned j) const { return element(i, j); }
41
42         DynamicMatrix &operator*=(T);
43         DynamicMatrix &operator/=(T);
44         DynamicMatrix &operator+=(const DynamicMatrix &);
45         DynamicMatrix &operator-=(const DynamicMatrix &);
46
47         DynamicMatrix &exchange_rows(unsigned, unsigned);
48         DynamicMatrix &multiply_row(unsigned, T);
49         DynamicMatrix &add_row(unsigned, unsigned, T);
50
51         DynamicMatrix &invert();
52 };
53
54
55 template<typename T>
56 inline DynamicMatrix<T>::DynamicMatrix(unsigned r, unsigned c):
57         rows_(r),
58         columns_(c),
59         data(new T[rows_*columns_])
60 {
61         std::fill(data, data+rows_*columns_, T());
62 }
63
64 template<typename T>
65 inline DynamicMatrix<T>::DynamicMatrix(unsigned r, unsigned c, const T *d):
66         rows_(r),
67         columns_(c),
68         data(new T[rows_*columns_])
69 {
70         std::copy(d, d+rows_*columns_, data);
71 }
72
73 template<typename T>
74 inline DynamicMatrix<T>::DynamicMatrix(const DynamicMatrix &other):
75         rows_(other.rows()),
76         columns_(other.columns()),
77         data(new T[rows_*columns_])
78 {
79         std::copy(other.data, other.data+rows_*columns_, data);
80 }
81
82 template<typename T>
83 inline DynamicMatrix<T> &DynamicMatrix<T>::operator=(const DynamicMatrix &other)
84 {
85         if(rows_!=other.rows() || columns_!=other.columns())
86         {
87                 delete[] data;
88                 rows_ = other.rows();
89                 columns_ = other.columns();
90                 data = new T[rows_*columns_];
91         }
92
93         std::copy(other.data, other.data+rows_*columns_, data);
94
95         return *this;
96 }
97
98 template<typename T>
99 inline DynamicMatrix<T>::~DynamicMatrix()
100 {
101         delete[] data;
102 }
103
104 template<typename T>
105 inline T &DynamicMatrix<T>::element(unsigned i, unsigned j)
106 {
107         if(i>=rows_ || j>=columns_)
108                 throw std::out_of_range("DynamicMatrix::element");
109
110         return data[i+rows_*j];
111 }
112
113 template<typename T>
114 inline const T &DynamicMatrix<T>::element(unsigned i, unsigned j) const
115 {
116         if(i>=rows_ || j>=columns_)
117                 throw std::out_of_range("DynamicMatrix::element");
118
119         return data[i+rows_*j];
120 }
121
122 template<typename T>
123 inline DynamicMatrix<T> &DynamicMatrix<T>::operator*=(T s)
124 {
125         for(unsigned i=0; i<rows_*columns_; ++i)
126                 data[i] *= s;
127         return *this;
128 }
129
130 template<typename T>
131 inline DynamicMatrix<T> operator*(const DynamicMatrix<T> &m, T s)
132 {
133         DynamicMatrix<T> r(m);
134         return r *= s;
135 }
136
137 template<typename T>
138 inline DynamicMatrix<T> operator*(T s, const DynamicMatrix<T> &m)
139 {
140         return m*s;
141 }
142
143 template<typename T>
144 inline DynamicMatrix<T> operator*(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
145 {
146         if(m1.columns()!=m2.rows())
147                 throw size_mismatch("matrix*matrix");
148
149         DynamicMatrix<T> r(m1.rows(), m2.columns());
150         for(unsigned i=0; i<m1.rows(); ++i)
151                 for(unsigned j=0; j<m2.columns(); ++j)
152                         for(unsigned k=0; k<m1.columns(); ++k)
153                                 r(i, j) += m1(i, k)*m2(k, j);
154
155         return r;
156 }
157
158 template<typename T>
159 inline DynamicVector<T> operator*(const DynamicMatrix<T> &m, const DynamicVector<T> &v)
160 {
161         if(m.columns()!=v.size())
162                 throw size_mismatch("matrix*vector");
163
164         DynamicVector<T> r(m.rows());
165         for(unsigned i=0; i<m.rows(); ++i)
166                 for(unsigned j=0; j<m.columns(); ++j)
167                         r[i] += m(i, j)*v[j];
168
169         return r;
170 }
171
172 template<typename T>
173 inline DynamicVector<T> operator*(const DynamicVector<T> &v, const DynamicMatrix<T> &m)
174 {
175         if(v.size()!=m.rows())
176                 throw size_mismatch("vector*matrix");
177
178         DynamicVector<T> r(m.columns());
179         for(unsigned j=0; j<m.columns(); ++j)
180                 for(unsigned i=0; i<m.rows(); ++i)
181                         r[i] += v[i]*m(i, j);
182
183         return r;
184 }
185
186 template<typename T>
187 inline DynamicMatrix<T> &DynamicMatrix<T>::operator/=(T s)
188 {
189         for(unsigned i=0; i<rows_*columns_; ++i)
190                 data[i] /= s;
191         return *this;
192 }
193
194 template<typename T>
195 inline DynamicMatrix<T> operator/(const DynamicMatrix<T> &m, T s)
196 {
197         DynamicMatrix<T> r(m);
198         return r /= s;
199 }
200
201 template<typename T>
202 inline DynamicMatrix<T> &DynamicMatrix<T>::operator+=(const DynamicMatrix<T> &m)
203 {
204         if(rows_!=m.rows_ || columns_!=m.columns_)
205                 throw size_mismatch("matrix+matrix");
206
207         for(unsigned i=0; i<rows_*columns_; ++i)
208                 data[i] += m.data[i];
209
210         return *this;
211 }
212
213 template<typename T>
214 inline DynamicMatrix<T> operator+(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
215 {
216         DynamicMatrix<T> r(m1);
217         return r += m2;
218 }
219
220 template<typename T>
221 inline DynamicMatrix<T> &DynamicMatrix<T>::operator-=(const DynamicMatrix<T> &m)
222 {
223         if(rows_!=m.rows_ || columns_!=m.columns_)
224                 throw size_mismatch("matrix-matrix");
225
226         for(unsigned i=0; i<rows_*columns_; ++i)
227                 data[i] += m.data[i];
228
229         return *this;
230 }
231
232 template<typename T>
233 inline DynamicMatrix<T> operator-(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
234 {
235         DynamicMatrix<T> r(m1);
236         return r += m2;
237 }
238
239 template<typename T>
240 inline bool operator==(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
241 {
242         if(m1.rows()!=m2.rows() || m1.columns()!=m2.columns())
243                 throw size_mismatch("matrix==matrix");
244
245         for(unsigned i=0; i<m1.rows(); ++i)
246                 for(unsigned j=0; j<m1.columns(); ++j)
247                         if(m1(i, j)!=m2(i, j))
248                                 return false;
249
250         return true;
251 }
252
253 template<typename T>
254 inline DynamicMatrix<T> &DynamicMatrix<T>::exchange_rows(unsigned i, unsigned j)
255 {
256         if(i>=rows_ || j>=rows_)
257                 throw std::out_of_range("DynamicMatrix::exchange_rows");
258
259         using std::swap;
260         for(unsigned k=0; k<columns_; ++k)
261                 swap(element(i, k), element(j, k));
262
263         return *this;
264 }
265
266 template<typename T>
267 inline DynamicMatrix<T> &DynamicMatrix<T>::multiply_row(unsigned i, T s)
268 {
269         if(i>=rows_)
270                 throw std::out_of_range("DynamicMatrix::multiply_row");
271
272         for(unsigned k=0; k<columns_; ++k)
273                 element(i, k) *= s;
274
275         return *this;
276 }
277
278 template<typename T>
279 inline DynamicMatrix<T> &DynamicMatrix<T>::add_row(unsigned i, unsigned j, T s)
280 {
281         if(i>=rows_ || j>=rows_)
282                 throw std::out_of_range("DynamicMatrix::exchange_rows");
283
284         for(unsigned k=0; k<columns_; ++k)
285                 element(j, k) += element(i, k)*s;
286
287         return *this;
288 }
289
290 template<typename T>
291 inline DynamicMatrix<T> transpose(const DynamicMatrix<T> &m)
292 {
293         DynamicMatrix<T> r(m.columns(), m.rows());
294         for(unsigned i=0; i<m.rows(); ++i)
295                 for(unsigned j=0; j<m.columns(); ++j)
296                         r(j, i) = m(i, j);
297         return r;
298 }
299
300 template<typename T>
301 inline DynamicMatrix<T> &DynamicMatrix<T>::invert()
302 {
303         if(rows_!=columns_)
304                 throw size_mismatch("matrix invert");
305
306         DynamicMatrix<T> r(rows_, columns_);
307         for(unsigned i=0; i<rows_; ++i)
308                 r(i, i) = T(1);
309
310         return invert_matrix(*this, r);
311 }
312
313 template<typename T>
314 inline DynamicMatrix<T> invert(const DynamicMatrix<T> &m)
315 {
316         DynamicMatrix<T> r = m;
317         return r.invert();
318 }
319
320 } // namespace LinAl
321 } // namespace Msp
322
323 #endif