--- /dev/null
+#ifndef MSP_LINAL_DYNAMICMATRIX_H_
+#define MSP_LINAL_DYNAMICMATRIX_H_
+
+#include <algorithm>
+#include "dynamicvector.h"
+#include "matrixops.h"
+
+namespace Msp {
+namespace LinAl {
+
+/**
+A general mathematical matrix. Size is specified at runtime. May be slower
+than a fixed-size matrix. There are no compile-time diagnostics for mismatched
+matrix sizes.
+*/
+template<typename T>
+class DynamicMatrix
+{
+public:
+ typedef T ElementType;
+
+private:
+ unsigned rows_;
+ unsigned columns_;
+ T *data;
+
+public:
+ DynamicMatrix(unsigned, unsigned);
+ DynamicMatrix(unsigned, unsigned, const T *);
+ DynamicMatrix(const DynamicMatrix &);
+ DynamicMatrix &operator=(const DynamicMatrix &);
+ ~DynamicMatrix();
+
+ unsigned rows() const { return rows_; }
+ unsigned columns() const { return columns_; }
+
+ T &element(unsigned, unsigned);
+ const T &element(unsigned, unsigned) const;
+ T &operator()(unsigned i, unsigned j) { return element(i, j); }
+ const T &operator()(unsigned i, unsigned j) const { return element(i, j); }
+
+ DynamicMatrix &operator*=(T);
+ DynamicMatrix &operator/=(T);
+ DynamicMatrix &operator+=(const DynamicMatrix &);
+ DynamicMatrix &operator-=(const DynamicMatrix &);
+
+ DynamicMatrix &exchange_rows(unsigned, unsigned);
+ DynamicMatrix &multiply_row(unsigned, T);
+ DynamicMatrix &add_row(unsigned, unsigned, T);
+
+ DynamicMatrix &invert();
+};
+
+
+template<typename T>
+inline DynamicMatrix<T>::DynamicMatrix(unsigned r, unsigned c):
+ rows_(r),
+ columns_(c),
+ data(new T[rows_*columns_])
+{
+ std::fill(data, data+rows_*columns_, T());
+}
+
+template<typename T>
+inline DynamicMatrix<T>::DynamicMatrix(unsigned r, unsigned c, const T *d):
+ rows_(r),
+ columns_(c),
+ data(new T[rows_*columns_])
+{
+ std::copy(d, d+rows_*columns_, data);
+}
+
+template<typename T>
+inline DynamicMatrix<T>::DynamicMatrix(const DynamicMatrix &other):
+ rows_(other.rows()),
+ columns_(other.columns()),
+ data(new T[rows_*columns_])
+{
+ std::copy(other.data, other.data+rows_*columns_, data);
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::operator=(const DynamicMatrix &other)
+{
+ if(rows_!=other.rows() || columns_!=other.columns())
+ {
+ delete[] data;
+ rows_ = other.rows();
+ columns_ = other.columns();
+ data = new T[rows_*columns_];
+ }
+
+ std::copy(other.data, other.data+rows_*columns_, data);
+
+ return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T>::~DynamicMatrix()
+{
+ delete[] data;
+}
+
+template<typename T>
+inline T &DynamicMatrix<T>::element(unsigned i, unsigned j)
+{
+ if(i>=rows_ || j>=columns_)
+ throw std::out_of_range("DynamicMatrix::element");
+
+ return data[i+rows_*j];
+}
+
+template<typename T>
+inline const T &DynamicMatrix<T>::element(unsigned i, unsigned j) const
+{
+ if(i>=rows_ || j>=columns_)
+ throw std::out_of_range("DynamicMatrix::element");
+
+ return data[i+rows_*j];
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::operator*=(T s)
+{
+ for(unsigned i=0; i<rows_*columns_; ++i)
+ data[i] *= s;
+ return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator*(const DynamicMatrix<T> &m, T s)
+{
+ DynamicMatrix<T> r(m);
+ return r *= s;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator*(T s, const DynamicMatrix<T> &m)
+{
+ return m*s;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator*(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
+{
+ if(m1.columns()!=m2.rows())
+ throw size_mismatch("matrix*matrix");
+
+ DynamicMatrix<T> r(m1.rows(), m2.columns());
+ for(unsigned i=0; i<m1.rows(); ++i)
+ for(unsigned j=0; j<m2.columns(); ++j)
+ for(unsigned k=0; k<m1.columns(); ++k)
+ r(i, j) += m1(i, k)*m2(k, j);
+
+ return r;
+}
+
+template<typename T>
+inline DynamicVector<T> operator*(const DynamicMatrix<T> &m, const DynamicVector<T> &v)
+{
+ if(m.columns()!=v.size())
+ throw size_mismatch("matrix*vector");
+
+ DynamicVector<T> r(m.rows());
+ for(unsigned i=0; i<m.rows(); ++i)
+ for(unsigned j=0; j<m.columns(); ++j)
+ r[i] += m(i, j)*v[j];
+
+ return r;
+}
+
+template<typename T>
+inline DynamicVector<T> operator*(const DynamicVector<T> &v, const DynamicMatrix<T> &m)
+{
+ if(v.size()!=m.rows())
+ throw size_mismatch("vector*matrix");
+
+ DynamicVector<T> r(m.columns());
+ for(unsigned j=0; j<m.columns(); ++j)
+ for(unsigned i=0; i<m.rows(); ++i)
+ r[i] += v[i]*m(i, j);
+
+ return r;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::operator/=(T s)
+{
+ for(unsigned i=0; i<rows_*columns_; ++i)
+ data[i] /= s;
+ return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator/(const DynamicMatrix<T> &m, T s)
+{
+ DynamicMatrix<T> r(m);
+ return r /= s;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::operator+=(const DynamicMatrix<T> &m)
+{
+ if(rows_!=m.rows_ || columns_!=m.columns_)
+ throw size_mismatch("matrix+matrix");
+
+ for(unsigned i=0; i<rows_*columns_; ++i)
+ data[i] += m.data[i];
+
+ return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator+(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
+{
+ DynamicMatrix<T> r(m1);
+ return r += m2;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::operator-=(const DynamicMatrix<T> &m)
+{
+ if(rows_!=m.rows_ || columns_!=m.columns_)
+ throw size_mismatch("matrix-matrix");
+
+ for(unsigned i=0; i<rows_*columns_; ++i)
+ data[i] += m.data[i];
+
+ return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator-(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
+{
+ DynamicMatrix<T> r(m1);
+ return r += m2;
+}
+
+template<typename T>
+inline bool operator==(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
+{
+ if(m1.rows()!=m2.rows() || m1.columns()!=m2.columns())
+ throw size_mismatch("matrix==matrix");
+
+ for(unsigned i=0; i<m1.rows(); ++i)
+ for(unsigned j=0; j<m1.columns(); ++j)
+ if(m1(i, j)!=m2(i, j))
+ return false;
+
+ return true;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::exchange_rows(unsigned i, unsigned j)
+{
+ if(i>=rows_ || j>=rows_)
+ throw std::out_of_range("DynamicMatrix::exchange_rows");
+
+ using std::swap;
+ for(unsigned k=0; k<columns_; ++k)
+ swap(element(i, k), element(j, k));
+
+ return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::multiply_row(unsigned i, T s)
+{
+ if(i>=rows_)
+ throw std::out_of_range("DynamicMatrix::multiply_row");
+
+ for(unsigned k=0; k<columns_; ++k)
+ element(i, k) *= s;
+
+ return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::add_row(unsigned i, unsigned j, T s)
+{
+ if(i>=rows_ || j>=rows_)
+ throw std::out_of_range("DynamicMatrix::exchange_rows");
+
+ for(unsigned k=0; k<columns_; ++k)
+ element(j, k) += element(i, k)*s;
+
+ return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> transpose(const DynamicMatrix<T> &m)
+{
+ DynamicMatrix<T> r(m.columns(), m.rows());
+ for(unsigned i=0; i<m.rows(); ++i)
+ for(unsigned j=0; j<m.columns(); ++j)
+ r(j, i) = m(i, j);
+ return r;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::invert()
+{
+ if(rows_!=columns_)
+ throw size_mismatch("matrix invert");
+
+ DynamicMatrix<T> r(rows_, columns_);
+ for(unsigned i=0; i<rows_; ++i)
+ r(i, i) = T(1);
+
+ return invert_matrix(*this, r);
+}
+
+template<typename T>
+inline DynamicMatrix<T> invert(const DynamicMatrix<T> &m)
+{
+ DynamicMatrix<T> r = m;
+ return r.invert();
+}
+
+} // namespace LinAl
+} // namespace Msp
+
+#endif
projects / libs / math.git / blobdiff