1 #ifndef MSP_LINAL_DYNAMICMATRIX_H_
2 #define MSP_LINAL_DYNAMICMATRIX_H_
6 #include "dynamicvector.h"
13 A general mathematical matrix. Size is specified at runtime. May be slower
14 than a fixed-size matrix. There are no compile-time diagnostics for mismatched
21 typedef T ElementType;
29 DynamicMatrix(unsigned, unsigned);
30 DynamicMatrix(unsigned, unsigned, const T *);
31 DynamicMatrix(const DynamicMatrix &);
32 DynamicMatrix &operator=(const DynamicMatrix &);
35 unsigned rows() const { return rows_; }
36 unsigned columns() const { return columns_; }
38 T &element(unsigned, unsigned);
39 const T &element(unsigned, unsigned) const;
40 T &operator()(unsigned i, unsigned j) { return element(i, j); }
41 const T &operator()(unsigned i, unsigned j) const { return element(i, j); }
43 DynamicMatrix &operator*=(T);
44 DynamicMatrix &operator/=(T);
45 DynamicMatrix &operator+=(const DynamicMatrix &);
46 DynamicMatrix &operator-=(const DynamicMatrix &);
48 DynamicMatrix &exchange_columns(unsigned, unsigned);
49 DynamicMatrix &multiply_column(unsigned, T);
50 DynamicMatrix &add_column(unsigned, unsigned, T);
52 DynamicMatrix &invert();
57 inline DynamicMatrix<T>::DynamicMatrix(unsigned r, unsigned c):
60 data(new T[rows_*columns_])
62 std::fill(data, data+rows_*columns_, T());
66 inline DynamicMatrix<T>::DynamicMatrix(unsigned r, unsigned c, const T *d):
69 data(new T[rows_*columns_])
71 std::copy(d, d+rows_*columns_, data);
75 inline DynamicMatrix<T>::DynamicMatrix(const DynamicMatrix &other):
77 columns_(other.columns()),
78 data(new T[rows_*columns_])
80 std::copy(other.data, other.data+rows_*columns_, data);
84 inline DynamicMatrix<T> &DynamicMatrix<T>::operator=(const DynamicMatrix &other)
86 if(rows_!=other.rows() || columns_!=other.columns())
90 columns_ = other.columns();
91 data = new T[rows_*columns_];
94 std::copy(other.data, other.data+rows_*columns_, data);
100 inline DynamicMatrix<T>::~DynamicMatrix()
106 inline T &DynamicMatrix<T>::element(unsigned i, unsigned j)
108 if(i>=rows_ || j>=columns_)
109 throw std::out_of_range("DynamicMatrix::element");
111 return data[i+rows_*j];
115 inline const T &DynamicMatrix<T>::element(unsigned i, unsigned j) const
117 if(i>=rows_ || j>=columns_)
118 throw std::out_of_range("DynamicMatrix::element");
120 return data[i+rows_*j];
124 inline DynamicMatrix<T> &DynamicMatrix<T>::operator*=(T s)
126 for(unsigned i=0; i<rows_*columns_; ++i)
132 inline DynamicMatrix<T> operator*(const DynamicMatrix<T> &m, T s)
134 DynamicMatrix<T> r(m);
139 inline DynamicMatrix<T> operator*(T s, const DynamicMatrix<T> &m)
145 inline DynamicMatrix<T> operator*(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
147 if(m1.columns()!=m2.rows())
148 throw size_mismatch("matrix*matrix");
150 DynamicMatrix<T> r(m1.rows(), m2.columns());
151 for(unsigned i=0; i<m1.rows(); ++i)
152 for(unsigned j=0; j<m2.columns(); ++j)
153 for(unsigned k=0; k<m1.columns(); ++k)
154 r(i, j) += m1(i, k)*m2(k, j);
160 inline DynamicVector<T> operator*(const DynamicMatrix<T> &m, const DynamicVector<T> &v)
162 if(m.columns()!=v.size())
163 throw size_mismatch("matrix*vector");
165 DynamicVector<T> r(m.rows());
166 for(unsigned i=0; i<m.rows(); ++i)
167 for(unsigned j=0; j<m.columns(); ++j)
168 r[i] += m(i, j)*v[j];
174 inline DynamicVector<T> operator*(const DynamicVector<T> &v, const DynamicMatrix<T> &m)
176 if(v.size()!=m.rows())
177 throw size_mismatch("vector*matrix");
179 DynamicVector<T> r(m.columns());
180 for(unsigned j=0; j<m.columns(); ++j)
181 for(unsigned i=0; i<m.rows(); ++i)
182 r[i] += v[i]*m(i, j);
188 inline DynamicMatrix<T> &DynamicMatrix<T>::operator/=(T s)
190 for(unsigned i=0; i<rows_*columns_; ++i)
196 inline DynamicMatrix<T> operator/(const DynamicMatrix<T> &m, T s)
198 DynamicMatrix<T> r(m);
203 inline DynamicMatrix<T> &DynamicMatrix<T>::operator+=(const DynamicMatrix<T> &m)
205 if(rows_!=m.rows_ || columns_!=m.columns_)
206 throw size_mismatch("matrix+matrix");
208 for(unsigned i=0; i<rows_*columns_; ++i)
209 data[i] += m.data[i];
215 inline DynamicMatrix<T> operator+(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
217 DynamicMatrix<T> r(m1);
222 inline DynamicMatrix<T> &DynamicMatrix<T>::operator-=(const DynamicMatrix<T> &m)
224 if(rows_!=m.rows_ || columns_!=m.columns_)
225 throw size_mismatch("matrix-matrix");
227 for(unsigned i=0; i<rows_*columns_; ++i)
228 data[i] += m.data[i];
234 inline DynamicMatrix<T> operator-(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
236 DynamicMatrix<T> r(m1);
241 inline bool operator==(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
243 if(m1.rows()!=m2.rows() || m1.columns()!=m2.columns())
244 throw size_mismatch("matrix==matrix");
246 for(unsigned i=0; i<m1.rows(); ++i)
247 for(unsigned j=0; j<m1.columns(); ++j)
248 if(m1(i, j)!=m2(i, j))
255 inline DynamicMatrix<T> &DynamicMatrix<T>::exchange_columns(unsigned i, unsigned j)
257 if(i>=columns_ || j>=columns_)
258 throw std::out_of_range("DynamicMatrix::exchange_columns");
261 for(unsigned k=0; k<columns_; ++k)
262 swap(element(k, i), element(k, j));
268 inline DynamicMatrix<T> &DynamicMatrix<T>::multiply_column(unsigned i, T s)
271 throw std::out_of_range("DynamicMatrix::multiply_column");
273 for(unsigned k=0; k<columns_; ++k)
280 inline DynamicMatrix<T> &DynamicMatrix<T>::add_column(unsigned i, unsigned j, T s)
282 if(i>=columns_ || j>=columns_)
283 throw std::out_of_range("DynamicMatrix::exchange_columns");
285 for(unsigned k=0; k<columns_; ++k)
286 element(k, j) += element(k, i)*s;
292 inline DynamicMatrix<T> transpose(const DynamicMatrix<T> &m)
294 DynamicMatrix<T> r(m.columns(), m.rows());
295 for(unsigned i=0; i<m.rows(); ++i)
296 for(unsigned j=0; j<m.columns(); ++j)
302 inline DynamicMatrix<T> &DynamicMatrix<T>::invert()
305 throw size_mismatch("matrix invert");
307 DynamicMatrix<T> r(rows_, columns_);
308 for(unsigned i=0; i<rows_; ++i)
311 gauss_jordan(*this, r);
316 inline DynamicMatrix<T> invert(const DynamicMatrix<T> &m)
318 DynamicMatrix<T> r = m;
323 inline std::ostream &operator<<(std::ostream &s, const DynamicMatrix<T> &m)
325 s << "DynamicMatrix" << m.rows() << 'x' << m.columns() << '(';
326 for(unsigned i=0; i<m.columns(); ++i)
331 for(unsigned j=0; j<m.rows(); ++j)