From a3d9772a9fd483278b6248a811dc6f4ce968892b Mon Sep 17 00:00:00 2001 From: Mikko Rasa Date: Tue, 14 May 2013 00:11:38 +0300 Subject: [PATCH] Another big batch of stuff This includes most of the missing things. --- source/linal/matrix.h | 139 ++++++++++++++++++++++++++++++------ source/linal/squarematrix.h | 86 +++++++++++++++++++++- source/linal/vector.h | 9 +++ source/linal/vector3.h | 8 +-- 4 files changed, 216 insertions(+), 26 deletions(-) diff --git a/source/linal/matrix.h b/source/linal/matrix.h index 4fc2bc6..94a6e09 100644 --- a/source/linal/matrix.h +++ b/source/linal/matrix.h @@ -1,6 +1,7 @@ #ifndef MSP_LINAL_MATRIX_H_ #define MSP_LINAL_MATRIX_H_ +#include #include "vector.h" namespace Msp { @@ -20,9 +21,12 @@ public: Matrix(const T *); template Matrix(const Matrix &); + static Matrix from_columns(const Vector *); static Matrix from_rows(const Vector *); + T &element(unsigned, unsigned); + const T &element(unsigned, unsigned) const; T &operator()(unsigned, unsigned); const T &operator()(unsigned, unsigned) const; @@ -31,21 +35,72 @@ public: Matrix &operator+=(const Matrix &); Matrix &operator-=(const Matrix &); - Matrix transpose() const; + Matrix &exchange_rows(unsigned, unsigned); + Matrix &multiply_row(unsigned, T); + Matrix &add_row(unsigned, unsigned, T); }; template -inline T &Matrix::operator()(unsigned i, unsigned j) +inline Matrix::Matrix() +{ + std::fill(data, data+M*N, T()); +} + +template +inline Matrix::Matrix(const T *d) +{ + std::copy(d, d+M*N, data); +} + +template +template +inline Matrix::Matrix(const Matrix &m) +{ + std::copy(m.data, m.data+M*N, data); +} + +template +inline Matrix Matrix::from_columns(const Vector *v) +{ + Matrix m; + for(unsigned i=0; i +inline Matrix Matrix::from_rows(const Vector *v) +{ + Matrix m; + for(unsigned i=0; i +inline T &Matrix::element(unsigned i, unsigned j) { return data[i+M*j]; } template -inline const T &Matrix::operator()(unsigned i, unsigned j) const +inline const T &Matrix::element(unsigned i, unsigned j) const { return data[i+M*j]; } +template +inline T &Matrix::operator()(unsigned i, unsigned j) +{ + return element(i, j); +} + +template +inline const T &Matrix::operator()(unsigned i, unsigned j) const +{ + return element(i, j); +} + template inline Matrix &Matrix::operator*=(T s) { @@ -67,6 +122,37 @@ inline Matrix operator*(T s, const Matrix &m) return m*s; } +template +inline Matrix operator*(const Matrix &m1, const Matrix &m2) +{ + Matrix r; + for(unsigned i=0; i +inline Vector operator*(const Matrix &m, const Vector &v) +{ + Vector r; + for(unsigned i=0; i +inline Vector operator*(const Vector &v, const Matrix &m) +{ + Vector r; + for(unsigned j=0; j inline Matrix &Matrix::operator/=(T s) { @@ -112,34 +198,47 @@ inline Matrix operator-(const Matrix &m1, const Matrix -Matrix operator*(const Matrix &m1, const Matrix &m2) +template +inline bool operator==(const Matrix &a, const Matrix &b) { - Matrix r; - for(unsigned i=0; i -Vector operator*(const Matrix &m, const Vector &v) +inline Matrix &Matrix::exchange_rows(unsigned i, unsigned j) { - Vector r; - for(unsigned i=0; i -Vector operator*(const Vector &v, const Matrix &m) +inline Matrix &Matrix::multiply_row(unsigned i, T s) { - Vector r; + for(unsigned k=0; k +inline Matrix &Matrix::add_row(unsigned i, unsigned j, T s) +{ + for(unsigned k=0; k +inline Matrix transpose(const Matrix &m) +{ + Matrix r; for(unsigned j=0; j #include "matrix.h" namespace Msp { namespace LinAl { +class not_invertible: public std::domain_error +{ +public: + not_invertible(): domain_error(std::string()) { } + virtual ~not_invertible() throw() { } +}; + template class SquareMatrix: public Matrix { public: - SquareMatrix(); + SquareMatrix() { } SquareMatrix(const T *); template SquareMatrix(const Matrix &); + static SquareMatrix identity(); SquareMatrix &operator*=(const SquareMatrix &); - void invert(); + SquareMatrix &invert(); }; +template +SquareMatrix::SquareMatrix(const T *d): + Matrix(d) +{ } + +template +template +SquareMatrix::SquareMatrix(const Matrix &m): + Matrix(m) +{ } + +template +inline SquareMatrix SquareMatrix::identity() +{ + SquareMatrix m; + for(unsigned i=0; i +SquareMatrix &SquareMatrix::operator*=(const SquareMatrix &m) +{ + Matrix::operator*=(m); + return *this; +} + +template +SquareMatrix &SquareMatrix::invert() +{ + SquareMatrix r = identity(); + for(unsigned i=0; ielement(i, i)==T(0)) + { + unsigned pivot = i; + for(unsigned j=i+1; jelement(j, i))>abs(this->element(pivot, i))) + pivot = j; + + if(pivot==i) + throw not_invertible(); + + 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; +} + +template +inline SquareMatrix invert(const SquareMatrix &m) +{ + SquareMatrix r = m; + return r.invert(); +} + } // namespace LinAl } // namespace Msp diff --git a/source/linal/vector.h b/source/linal/vector.h index 933611e..1a79bd8 100644 --- a/source/linal/vector.h +++ b/source/linal/vector.h @@ -141,6 +141,15 @@ inline Vector operator-(const Vector &v) return r; } +template +inline bool operator==(const Vector &v, const Vector &w) +{ + for(unsigned i=0; i inline T inner_product(const Vector &v1, const Vector &v2) { diff --git a/source/linal/vector3.h b/source/linal/vector3.h index 90a5ff5..0197dbc 100644 --- a/source/linal/vector3.h +++ b/source/linal/vector3.h @@ -40,18 +40,18 @@ inline Vector3::Vector3(const Vector &v): { } template -inline T dot(const Vector3 &v1, const Vector3 &v2) +inline T dot(const Vector &v1, const Vector &v2) { return inner_product(v1, v2); } template -inline Vector3 cross(const Vector3 &v1, const Vector3 &v2) +inline Vector cross(const Vector &v1, const Vector &v2) { - return Vector3(v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1]-v1[1]*v2[0]); + return Vector(v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1]-v1[1]*v2[0]); } } // namespace LinAl -} // namespace LinAl +} // namespace Msp #endif -- 2.45.2