libraries/libmatricxx
1
0
Fork 0
Code Issues Pull Requests Releases Activity

more operators more tests

Browse Source
This commit is contained in:
Marc Wäckerlin
2016-08-19 14:52:35 +00:00
parent 092d194b95
commit b9ccd5dee5
8 changed files with 363 additions and 59 deletions
Download Patch File Download Diff File Expand all files Collapse all files

180
src/matrix.hxx
View File

@@ -84,6 +84,41 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
///}
/// @name element access
///{
TYPE& at(size_t row, size_t column) {
assert(row<ROWS);
assert(column<COLUMNS);
return *((TYPE*)_c+row*COLUMNS+column);
}
const TYPE& at(size_t row, size_t column) const {
assert(row<ROWS);
assert(column<COLUMNS);
return *((TYPE*)_c+row*COLUMNS+column);
}
TYPE& operator()(size_t row, size_t column) {
return at(row, column);
}
const TYPE& operator()(size_t row, size_t column) const {
return at(row, column);
}
RowVector operator[](size_t row) {
assert(row<ROWS);
return RowVector(*this, (TYPE*)_c+row*COLUMNS);
}
const ConstRowVector operator[](size_t row) const {
assert(row<ROWS);
return ConstRowVector(*this, (TYPE*)_c+row*COLUMNS);
}
///}
/// @name operators
///{
@@ -138,33 +173,32 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
///}
/// @name element access
/// @name special operations
///{
TYPE& operator()(size_t row, size_t column) {
assert(row<ROWS);
assert(column<COLUMNS);
return *((TYPE*)_c+row*COLUMNS+column);
}
const TYPE& operator()(size_t row, size_t column) const {
assert(row<ROWS);
assert(column<COLUMNS);
return *((TYPE*)_c+row*COLUMNS+column);
}
RowVector operator[](size_t row) {
assert(row<ROWS);
return RowVector(*this, (TYPE*)_c+row*COLUMNS);
}
const ConstRowVector operator[](size_t row) const {
assert(row<ROWS);
return ConstRowVector(*this, (TYPE*)_c+row*COLUMNS);
}
///}
TYPE det() {
// rule of sarrus
TYPE res(0);
assert(ROWS==COLUMNS); // not really necessary
// 1) terms to add
for (size_t i(0); i<COLUMNS; ++i) {
TYPE term(1);
for (size_t j(0); j<ROWS; ++j)
term *= at((i+j)%COLUMNS, j);
res += term;
}
// 2) terms to subtract
for (size_t i(0); i<COLUMNS; ++i) {
TYPE term(1);
for (size_t j(0); j<ROWS; ++j)
term *= at((i+j)%COLUMNS, ROWS-j-1);
res -= term;
}
return res;
}
///}
//................................................................methods
protected:
@@ -196,18 +230,19 @@ template<typename TYPE, size_t TROWS=0, size_t TCOLUMNS=0> class Matrix:
/// @name construction
///{
Matrix(): Parent(TROWS, TCOLUMNS) {}
template<typename ...ARGS>
Matrix(ARGS...t): Parent(TROWS, TCOLUMNS, t...) {
static_assert(sizeof...(t)==TROWS*TCOLUMNS,
"wrong array size");
Matrix(): Parent(TROWS, TCOLUMNS) {
memset(Parent::_c, 0, Parent::MEM_SIZE);
}
Matrix(const Matrix& o): Matrix() {
memcpy(Parent::_c, o._c, Parent::MEM_SIZE);
}
template<typename ...ARGS>
Matrix(ARGS...t): Parent(TROWS, TCOLUMNS, t...) {
static_assert(sizeof...(t)==TROWS*TCOLUMNS, "wrong array size");
}
///}
/// @name operators
@@ -218,14 +253,6 @@ template<typename TYPE, size_t TROWS=0, size_t TCOLUMNS=0> class Matrix:
return *this;
}
Matrix<TYPE, TCOLUMNS, TROWS> T() const {
Matrix<TYPE, TCOLUMNS, TROWS> res;
for (size_t row(0); row<TROWS; ++row)
for (size_t column(0); column<TCOLUMNS; ++column)
res(column, row) = Parent::operator()(row, column);
return res;
}
Matrix operator-() const {
Matrix res;
for (TYPE *to((TYPE*)(res._c)+Parent::SIZE); to>(TYPE*)(Parent::_c); --to)
@@ -234,15 +261,30 @@ template<typename TYPE, size_t TROWS=0, size_t TCOLUMNS=0> class Matrix:
}
template<size_t NEWCOLUMNS>
Matrix operator*(const Matrix<TYPE, TCOLUMNS, NEWCOLUMNS>& o) {
Matrix<TYPE, TROWS, NEWCOLUMNS>
operator*(const Matrix<TYPE, TCOLUMNS, NEWCOLUMNS>& o) const {
Matrix<TYPE, TROWS, NEWCOLUMNS> res;
for (size_t i(0); i<TROWS; ++i)
for (size_t k(0); k<NEWCOLUMNS; ++k)
for (size_t j(0); j<TCOLUMNS; ++j)
res(i, k) += (*this)(i, j) * o(j, k);
res(i, k) += this->at(i, j) * o(j, k);
return res;
}
///}
/// @name special operations
///{
Matrix<TYPE, TCOLUMNS, TROWS> t() const {
Matrix<TYPE, TCOLUMNS, TROWS> res;
for (size_t row(0); row<TROWS; ++row)
for (size_t column(0); column<TCOLUMNS; ++column)
res(column, row) = this->at(row, column);
return res;
}
///}
};
//==============================================================================
@@ -259,12 +301,18 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
/// @name construction
///{
Matrix() = delete;
Matrix(size_t rows, size_t columns):
Parent(rows, columns) {
Parent::_c = new TYPE[rows*columns];
memset(Parent::_c, 0, Parent::MEM_SIZE);
}
Matrix(const Matrix& o): Matrix(o.ROWS, o.Parent::COLUMNS) {
memcpy(Parent::_c, o.Parent::_c, Parent::MEM_SIZE);
}
template<typename ...ARGS>
Matrix(size_t rows, size_t columns, ARGS...t):
@@ -273,10 +321,6 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
copy_args(Parent::_c, t...);
}
Matrix(const Matrix& o): Matrix(o.ROWS, o.Parent::COLUMNS) {
memcpy(Parent::_c, o.Parent::_c, Parent::MEM_SIZE);
}
///}
/// @name destruction
@@ -303,11 +347,11 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
return Parent::operator=(o);
}
Matrix T() const {
Matrix t() const {
Matrix res(Parent::COLUMNS, Parent::ROWS);
for (size_t row(0); row<Parent::ROWS; ++row)
for (size_t column(0); column<Parent::COLUMNS; ++column)
res(column, row) = Parent::operator()(row, column);
res(column, row) = this->at(row, column);
return res;
}
@@ -318,6 +362,16 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
return res;
}
Matrix operator*(const Matrix& o) const {
Matrix<TYPE> res(this->ROWS, o.COLUMNS);
assert(this->COLUMNS==o.ROWS);
for (size_t i(0); i<this->ROWS; ++i)
for (size_t k(0); k<o.COLUMNS; ++k)
for (size_t j(0); j<this->COLUMNS; ++j)
res(i, k) += this->at(i, j) * o(j, k);
return res;
}
///}
//................................................................methods
@@ -343,6 +397,38 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
//==============================================================================
template<typename TYPE, size_t ROWS, size_t COLUMNS>
Matrix<TYPE, ROWS, COLUMNS> operator+(const Matrix<TYPE, ROWS, COLUMNS>& a,
const Matrix<TYPE, ROWS, COLUMNS>& b) {
Matrix<TYPE, ROWS, COLUMNS> res(a);
res += b;
return res;
}
template<typename TYPE, size_t ROWS, size_t COLUMNS>
Matrix<TYPE, ROWS, COLUMNS> operator-(const Matrix<TYPE, ROWS, COLUMNS>& a,
const Matrix<TYPE, ROWS, COLUMNS>& b) {
Matrix<TYPE, ROWS, COLUMNS> res(a);
res -= b;
return res;
}
template<typename TYPE, size_t ROWS, size_t COLUMNS>
Matrix<TYPE, ROWS, COLUMNS> operator*(const TYPE& v,
const Matrix<TYPE, ROWS, COLUMNS>& m) {
Matrix<TYPE, ROWS, COLUMNS> res(m);
res *= v;
return res;
}
template<typename TYPE, size_t ROWS, size_t COLUMNS>
Matrix<TYPE, ROWS, COLUMNS> operator*(const Matrix<TYPE, ROWS, COLUMNS>& m,
const TYPE& v) {
Matrix<TYPE, ROWS, COLUMNS> res(m);
res *= v;
return res;
}
template<typename TYPE, size_t ROWS, size_t COLUMNS>
std::ostream& operator<<(std::ostream& s, const Matrix<TYPE, ROWS, COLUMNS>& m) {
for (size_t w = 0; w < m.ROWS; ++w) {