From: Mikko Rasa Date: Sun, 28 Aug 2016 23:57:09 +0000 (+0300) Subject: Add dynamically allocated versions of matrix and vector X-Git-Url: http://git.tdb.fi/?a=commitdiff_plain;h=5797fe0a5296952cb8f8643fdc6cabddee19a554;p=libs%2Fmath.git Add dynamically allocated versions of matrix and vector Sometimes it's not possible to know the dimensions at compile time, for example when calculating regressions from user-inputted data. --- diff --git a/source/linal/dummy.cpp b/source/linal/dummy.cpp index 05e6ff2..84d7395 100644 --- a/source/linal/dummy.cpp +++ b/source/linal/dummy.cpp @@ -1,6 +1,8 @@ #include "vector.h" #include "matrix.h" #include "squarematrix.h" +#include "dynamicvector.h" +#include "dynamicmatrix.h" /* This dummy file is needed to make Builder happy, and also serves as a syntax check by pulling in all the headers. */ diff --git a/source/linal/dynamicmatrix.h b/source/linal/dynamicmatrix.h new file mode 100644 index 0000000..3852424 --- /dev/null +++ b/source/linal/dynamicmatrix.h @@ -0,0 +1,323 @@ +#ifndef MSP_LINAL_DYNAMICMATRIX_H_ +#define MSP_LINAL_DYNAMICMATRIX_H_ + +#include +#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 +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 +inline DynamicMatrix::DynamicMatrix(unsigned r, unsigned c): + rows_(r), + columns_(c), + data(new T[rows_*columns_]) +{ + std::fill(data, data+rows_*columns_, T()); +} + +template +inline DynamicMatrix::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 +inline DynamicMatrix::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 +inline DynamicMatrix &DynamicMatrix::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 +inline DynamicMatrix::~DynamicMatrix() +{ + delete[] data; +} + +template +inline T &DynamicMatrix::element(unsigned i, unsigned j) +{ + if(i>=rows_ || j>=columns_) + throw std::out_of_range("DynamicMatrix::element"); + + return data[i+rows_*j]; +} + +template +inline const T &DynamicMatrix::element(unsigned i, unsigned j) const +{ + if(i>=rows_ || j>=columns_) + throw std::out_of_range("DynamicMatrix::element"); + + return data[i+rows_*j]; +} + +template +inline DynamicMatrix &DynamicMatrix::operator*=(T s) +{ + for(unsigned i=0; i +inline DynamicMatrix operator*(const DynamicMatrix &m, T s) +{ + DynamicMatrix r(m); + return r *= s; +} + +template +inline DynamicMatrix operator*(T s, const DynamicMatrix &m) +{ + return m*s; +} + +template +inline DynamicMatrix operator*(const DynamicMatrix &m1, const DynamicMatrix &m2) +{ + if(m1.columns()!=m2.rows()) + throw size_mismatch("matrix*matrix"); + + DynamicMatrix r(m1.rows(), m2.columns()); + for(unsigned i=0; i +inline DynamicVector operator*(const DynamicMatrix &m, const DynamicVector &v) +{ + if(m.columns()!=v.size()) + throw size_mismatch("matrix*vector"); + + DynamicVector r(m.rows()); + for(unsigned i=0; i +inline DynamicVector operator*(const DynamicVector &v, const DynamicMatrix &m) +{ + if(v.size()!=m.rows()) + throw size_mismatch("vector*matrix"); + + DynamicVector r(m.columns()); + for(unsigned j=0; j +inline DynamicMatrix &DynamicMatrix::operator/=(T s) +{ + for(unsigned i=0; i +inline DynamicMatrix operator/(const DynamicMatrix &m, T s) +{ + DynamicMatrix r(m); + return r /= s; +} + +template +inline DynamicMatrix &DynamicMatrix::operator+=(const DynamicMatrix &m) +{ + if(rows_!=m.rows_ || columns_!=m.columns_) + throw size_mismatch("matrix+matrix"); + + for(unsigned i=0; i +inline DynamicMatrix operator+(const DynamicMatrix &m1, const DynamicMatrix &m2) +{ + DynamicMatrix r(m1); + return r += m2; +} + +template +inline DynamicMatrix &DynamicMatrix::operator-=(const DynamicMatrix &m) +{ + if(rows_!=m.rows_ || columns_!=m.columns_) + throw size_mismatch("matrix-matrix"); + + for(unsigned i=0; i +inline DynamicMatrix operator-(const DynamicMatrix &m1, const DynamicMatrix &m2) +{ + DynamicMatrix r(m1); + return r += m2; +} + +template +inline bool operator==(const DynamicMatrix &m1, const DynamicMatrix &m2) +{ + if(m1.rows()!=m2.rows() || m1.columns()!=m2.columns()) + throw size_mismatch("matrix==matrix"); + + for(unsigned i=0; i +inline DynamicMatrix &DynamicMatrix::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 +inline DynamicMatrix &DynamicMatrix::multiply_row(unsigned i, T s) +{ + if(i>=rows_) + throw std::out_of_range("DynamicMatrix::multiply_row"); + + for(unsigned k=0; k +inline DynamicMatrix &DynamicMatrix::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 +inline DynamicMatrix transpose(const DynamicMatrix &m) +{ + DynamicMatrix r(m.columns(), m.rows()); + for(unsigned i=0; i +inline DynamicMatrix &DynamicMatrix::invert() +{ + if(rows_!=columns_) + throw size_mismatch("matrix invert"); + + DynamicMatrix r(rows_, columns_); + for(unsigned i=0; i +inline DynamicMatrix invert(const DynamicMatrix &m) +{ + DynamicMatrix r = m; + return r.invert(); +} + +} // namespace LinAl +} // namespace Msp + +#endif diff --git a/source/linal/dynamicvector.h b/source/linal/dynamicvector.h new file mode 100644 index 0000000..8839c33 --- /dev/null +++ b/source/linal/dynamicvector.h @@ -0,0 +1,249 @@ +#ifndef MSP_LINAL_DYNAMICVECTOR_H_ +#define MSP_LINAL_DYNAMICVECTOR_H_ + +#include +#include +#include + +namespace Msp { +namespace LinAl { + +class size_mismatch: public std::logic_error +{ +public: + size_mismatch(const std::string &w): std::logic_error(w) { } + virtual ~size_mismatch() throw() { } +}; + + +/** +A general mathematical vector. Size is specified at runtime. May be slower +than a fixed-size vector. There are no compile-time diagnostics for mismatched +vector sizes. +*/ +template +class DynamicVector +{ +public: + typedef T ElementType; + +private: + unsigned size_; + T *data; + +public: + DynamicVector(unsigned); + DynamicVector(unsigned, const T *); + DynamicVector(const DynamicVector &); + DynamicVector &operator=(const DynamicVector &); + ~DynamicVector(); + + unsigned size() const { return size_; } + + T &operator[](unsigned); + const T &operator[](unsigned) const; + + DynamicVector &operator*=(T); + DynamicVector &operator/=(T); + DynamicVector &operator+=(const DynamicVector &); + DynamicVector &operator-=(const DynamicVector &); + + T norm() const; + DynamicVector &normalize(); +}; + +template +inline DynamicVector::DynamicVector(unsigned s): + size_(s), + data(new T[size_]) +{ + std::fill(data, data+size_, T()); +} + +template +inline DynamicVector::DynamicVector(unsigned s, const T *d): + size_(s), + data(new T[size_]) +{ + std::copy(d, d+size_, data); +} + +template +inline DynamicVector::DynamicVector(const DynamicVector &other): + size_(other.size()), + data(new T[size_]) +{ + std::copy(other.data, other.data+size_, data); +} + +template +inline DynamicVector &DynamicVector::operator=(const DynamicVector &other) +{ + if(size_!=other.size()) + { + delete[] data; + size_ = other.size(); + data = new T[size_]; + } + + std::copy(other.data, other.data+size_, data); + + return *this; +} + +template +inline DynamicVector::~DynamicVector() +{ + delete[] data; +} + +template +inline T &DynamicVector::operator[](unsigned i) +{ + if(i>=size_) + throw std::out_of_range("DynamicVector::operator[]"); + + return data[i]; +} + +template +inline const T &DynamicVector::operator[](unsigned i) const +{ + if(i>=size_) + throw std::out_of_range("DynamicVector::operator[]"); + + return data[i]; +} + +template +inline DynamicVector &DynamicVector::operator*=(T s) +{ + for(unsigned i=0; i +inline DynamicVector operator*(const DynamicVector &v, T s) +{ + DynamicVector r(v); + return r *= s; +} + +template +inline DynamicVector operator*(T s, const DynamicVector &v) +{ + return v*s; +} + +template +inline DynamicVector &DynamicVector::operator/=(T s) +{ + for(unsigned i=0; i +inline DynamicVector operator/(const DynamicVector &v, T s) +{ + DynamicVector r(v); + return r /= s; +} + +template +inline DynamicVector &DynamicVector::operator+=(const DynamicVector &v) +{ + if(size_!=v.size()) + throw size_mismatch("vector+vector"); + + for(unsigned i=0; i +inline DynamicVector operator+(const DynamicVector &v1, const DynamicVector &v2) +{ + DynamicVector r(v1); + return r += v2; +} + +template +inline DynamicVector &DynamicVector::operator-=(const DynamicVector &v) +{ + if(size_!=v.size()) + throw size_mismatch("vector-vector"); + + for(unsigned i=0; i +inline DynamicVector operator-(const DynamicVector &v1, const DynamicVector &v2) +{ + DynamicVector r(v1); + return r -= v2; +} + +template +inline DynamicVector operator-(const DynamicVector &v) +{ + DynamicVector r(v); + for(unsigned i=0; i +inline bool operator==(const DynamicVector &v1, const DynamicVector &v2) +{ + if(v1.size()!=v2.size()) + throw size_mismatch("vector==vector"); + + for(unsigned i=0; i +inline T inner_product(const DynamicVector &v1, const DynamicVector &v2) +{ + if(v1.size()!=v2.size()) + throw size_mismatch("inner_product"); + + T r = T(); + for(unsigned i=0; i +inline T DynamicVector::norm() const +{ + using std::sqrt; + return sqrt(inner_product(*this, *this)); +} + +template +inline DynamicVector &DynamicVector::normalize() +{ + return *this /= norm(); +} + +template +inline DynamicVector normalize(const DynamicVector &v) +{ + DynamicVector r(v); + return r.normalize(); +} + +} // namespace LinAL +} // namespace Msp + +#endif diff --git a/source/linal/matrix.h b/source/linal/matrix.h index 929c482..a57ce08 100644 --- a/source/linal/matrix.h +++ b/source/linal/matrix.h @@ -13,6 +13,9 @@ A general mathematical matrix with M rows and N columns. template class Matrix { +public: + typedef T ElementType; + private: T data[M*N]; @@ -25,11 +28,17 @@ public: static Matrix from_columns(const Vector *); static Matrix from_rows(const Vector *); + unsigned rows() const { return M; } + unsigned columns() const { return N; } + T &element(unsigned i, unsigned j) { return data[i+M*j]; } const T &element(unsigned i, unsigned j) const { return data[i+M*j]; } T &operator()(unsigned i, unsigned j) { return element(i, j); } const T &operator()(unsigned i, unsigned j) const { return element(i, j); } + template + Matrix select(const Vector &, const Vector &) const; + template Matrix block(unsigned, unsigned) const; @@ -84,6 +93,17 @@ inline Matrix Matrix::from_rows(const Vector *v) return m; } +template +template +inline Matrix Matrix::select(const Vector &rows, const Vector &cols) const +{ + Matrix r; + for(unsigned j=0; j template inline Matrix Matrix::block(unsigned y, unsigned x) const diff --git a/source/linal/matrixops.h b/source/linal/matrixops.h new file mode 100644 index 0000000..790dfbf --- /dev/null +++ b/source/linal/matrixops.h @@ -0,0 +1,62 @@ +#ifndef MSP_LINAL_MATRIXOPS_H_ +#define MSP_LINAL_MATRIXOPS_H_ + +#include +#include + +namespace Msp { +namespace LinAl { + +class not_invertible: public std::domain_error +{ +public: + not_invertible(): domain_error(std::string()) { } + virtual ~not_invertible() throw() { } +}; + + +template +inline T &invert_matrix(T &m, T &r) +{ + typedef typename T::ElementType V; + using std::abs; + + for(unsigned i=0; iabs(m.element(pivot, i))) + pivot = j; + + if(m.element(pivot, i)==V(0)) + throw not_invertible(); + + if(pivot!=i) + { + m.exchange_rows(i, pivot); + r.exchange_rows(i, pivot); + } + + for(unsigned j=i+1; j0; ) + for(unsigned j=i; j-->0; ) + r.add_row(i, j, -m.element(j, i)); + + return m = r; +} + +} // namespace LinAl +} // namespace Msp + +#endif diff --git a/source/linal/squarematrix.h b/source/linal/squarematrix.h index 5117bb3..0a97e61 100644 --- a/source/linal/squarematrix.h +++ b/source/linal/squarematrix.h @@ -2,19 +2,12 @@ #define MSP_LINAL_SQUAREMATRIX_H_ #include -#include #include "matrix.h" +#include "matrixops.h" namespace Msp { namespace LinAl { -class not_invertible: public std::domain_error -{ -public: - not_invertible(): domain_error(std::string()) { } - virtual ~not_invertible() throw() { } -}; - /** A mathematical matrix with S rows and columns. Some operations are provided here that are only possible for square matrices. @@ -53,42 +46,8 @@ SquareMatrix &SquareMatrix::operator*=(const SquareMatrix &m) template SquareMatrix &SquareMatrix::invert() { - using std::abs; - SquareMatrix r = identity(); - for(unsigned i=0; ielement(j, i))>abs(this->element(pivot, i))) - pivot = j; - - if(this->element(pivot, i)==T(0)) - throw not_invertible(); - - if(pivot!=i) - { - this->exchange_rows(i, pivot); - r.exchange_rows(i, pivot); - } - - for(unsigned j=i+1; jelement(j, i)/this->element(i, i); - this->add_row(i, j, a); - r.add_row(i, j, a); - } - - T a = T(1)/this->element(i, i); - this->multiply_row(i, a); - r.multiply_row(i, a); - } - - for(unsigned i=S; i-->0; ) - for(unsigned j=i; j-->0; ) - r.add_row(i, j, -this->element(j, i)); - - return *this = r; + return invert_matrix(*this, r); } template diff --git a/source/linal/vector.h b/source/linal/vector.h index 6dfe6ca..152c26e 100644 --- a/source/linal/vector.h +++ b/source/linal/vector.h @@ -74,6 +74,8 @@ template class Vector: public VectorComponents { public: + typedef T ElementType; + Vector(); Vector(const T *); Vector(T, T); @@ -82,6 +84,8 @@ public: template Vector(const Vector &); + unsigned size() const { return N; } + template Vector slice(unsigned) const;