From cbcb6aeb1510ee1bd2d67d9d0b285df67f7088b5 Mon Sep 17 00:00:00 2001
From: Nao Pross <naopross@thearcway.org>
Date: Sat, 5 Oct 2019 17:59:32 +0200
Subject: New matrix data model (breaks everything)
---
include/mm/mmiterator.hpp | 248 -------------------
include/mm/mmmatrix.hpp | 595 +++++-----------------------------------------
include/mm/view.hpp | 66 +++++
test/matrix_example.cpp | 22 --
4 files changed, 132 insertions(+), 799 deletions(-)
delete mode 100644 include/mm/mmiterator.hpp
create mode 100644 include/mm/view.hpp
diff --git a/include/mm/mmiterator.hpp b/include/mm/mmiterator.hpp
deleted file mode 100644
index 4efe474..0000000
--- a/include/mm/mmiterator.hpp
+++ /dev/null
@@ -1,248 +0,0 @@
-#pragma once
-
-#include "debug.hpp"
-
-namespace mm::iter {
-
- template<typename T, std::size_t Rows, std::size_t Cols, class IterType, class Grid>
- class vector_iterator;
-
- template<typename T, std::size_t Rows, std::size_t Cols, class Grid>
- class basic_iterator;
-
- template<typename T, std::size_t N, class Grid>
- class diag_iterator;
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols, class IterType, class Grid>
-class mm::iter::vector_iterator
-{
-public:
-
- template<typename U, std::size_t R, std::size_t C, class G>
- friend class mm::iter::basic_iterator;
-
- template<typename U, std::size_t N, class G>
- friend class mm::iter::diag_iterator;
-
- vector_iterator(Grid& _M, std::size_t pos, std::size_t i = 0)
- : M(_M), position(pos), index(i) {}
-
-//#ifdef MM_IMPLICIT_CONVERSION_ITERATOR
- operator T&()
- {
- //npdebug("Calling +")
- return *(*this);
- }
-//#endif
-
- IterType operator++()
- {
- IterType it = cpy();
- ++index;
- return it;
- }
-
- IterType operator--()
- {
- IterType it = cpy();
- --index;
- return it;
- }
-
- IterType& operator++(int)
- {
- ++index;
- return ref();
- }
-
- IterType& operator--(int)
- {
- --index;
- return ref();
- }
-
- bool operator==(const IterType& other) const
- {
- return index == other.index;
- }
-
- bool operator!=(const IterType& other) const
- {
- return index != other.index;
- }
-
- bool ok() const
- {
- return index < size();
- }
-
- virtual std::size_t size() const = 0;
-
- virtual T& operator*() = 0;
- virtual T& operator[](std::size_t) = 0;
-
- virtual T& operator[](std::size_t) const = 0;
-
- IterType begin()
- {
- return IterType(M, position, 0);
- }
-
- virtual IterType end() = 0;
-
-protected:
-
- Grid& M; // grid mapping
-
- const std::size_t position; // fixed index, negative too for diagonal iterator
- std::size_t index; // variable index
-
- virtual IterType& ref() = 0;
- virtual IterType cpy() = 0;
-};
-
-
-/*
- * Scalar product
- */
-
-template<typename T,
- std::size_t R1, std::size_t C1,
- std::size_t R2, std::size_t C2,
- class IterType1, class IterType2,
- class Grid1, class Grid2>
-typename std::remove_const<T>::type operator*(const mm::iter::vector_iterator<T, R1, C1, IterType1, Grid1>& v,
- const mm::iter::vector_iterator<T, R2, C2, IterType2, Grid2>& w)
-{
- typename std::remove_const<T>::type out(0);
- const std::size_t N = std::min(v.size(), w.size());
-
- for(unsigned i = 0; i < N; ++i)
- out += v[i] * w[i];
-
- return out;
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols, class Grid>
-class mm::iter::basic_iterator : public mm::iter::vector_iterator<T, Rows, Cols, mm::iter::basic_iterator<T, Rows, Cols, Grid>, Grid>
-{
- bool direction;
-
- virtual mm::iter::basic_iterator<T, Rows, Cols, Grid>& ref() override
- {
- return *this;
- }
-
- virtual mm::iter::basic_iterator<T, Rows, Cols, Grid> cpy() override
- {
- return *this;
- }
-
-public:
-
- basic_iterator(Grid& A, std::size_t pos, std::size_t _index = 0, bool dir = true)
- : mm::iter::vector_iterator<T, Rows, Cols, mm::iter::basic_iterator<T, Rows, Cols, Grid>, Grid>
- (A, pos, _index), direction(dir)
- {
- //npdebug("Position: ", pos, ", Rows: ", Rows, " Cols: ", Cols, ", Direction: ", dir)
-
- if (direction)
- assert(pos < Rows);
- else
- assert(pos < Cols);
- }
-
- virtual std::size_t size() const
- {
- return (direction) ? Cols : Rows;
- }
-
-
- virtual T& operator*() override
- {
- return (direction) ?
- this->M.data[this->position * Cols + this->index] :
- this->M.data[this->index * Cols + this->position];
-
- }
-
- virtual T& operator[](std::size_t i) override
- {
- return (direction) ?
- this->M.data[this->position * Cols + i] :
- this->M.data[i * Cols + this->position];
- }
-
- virtual T& operator[](std::size_t i) const override
- {
- return (direction) ?
- this->M.data[this->position * Cols + i] :
- this->M.data[i * Cols + this->position];
- }
-
- virtual mm::iter::basic_iterator<T, Rows, Cols, Grid> end()
- {
- return mm::iter::basic_iterator<T, Rows, Cols, Grid>(this->M, this->position,
- (direction) ? Cols : Rows);
- }
-};
-
-template<typename T, std::size_t N, class Grid>
-class mm::iter::diag_iterator : public mm::iter::vector_iterator<T, N, N, mm::iter::diag_iterator<T, N, Grid>, Grid>
-{
- bool sign;
-
- virtual mm::iter::diag_iterator<T, N, Grid>& ref() override
- {
- return *this;
- }
-
- virtual mm::iter::diag_iterator<T, N, Grid> cpy() override
- {
- return *this;
- }
-
-public:
-
- diag_iterator(Grid& A, signed long int pos, std::size_t _index = 0)
- : mm::iter::vector_iterator<T, N, N, mm::iter::diag_iterator<T, N, Grid>, Grid>
- (A, static_cast<std::size_t>(labs(pos)), _index), sign(pos >= 0)
- {
- assert(this->position < N);
- }
-
- virtual std::size_t size() const
- {
- return N - this->position;
- }
-
- virtual T& operator*() override
- {
- return (sign) ?
- this->M.data[(this->index - this->position) * N + this->index] :
- this->M.data[this->index * N + (this->index + this->position)];
- }
-
- virtual T& operator[](std::size_t i) override
- {
- return (sign) ?
- this->M.data[(i - this->position) * N + i] :
- this->M.data[i * N + (i + this->position)];
- }
-
- virtual T& operator[](std::size_t i) const override
- {
- return (sign) ?
- this->M.data[(i - this->position) * N + i] :
- this->M.data[i * N + (i + this->position)];
- }
-
-
- virtual mm::iter::diag_iterator<T, N, Grid> end()
- {
- return mm::iter::diag_iterator<T, N, Grid>(this->M, this->position, N);
- }
-};
-
-
diff --git a/include/mm/mmmatrix.hpp b/include/mm/mmmatrix.hpp
index facbb48..d653d72 100644
--- a/include/mm/mmmatrix.hpp
+++ b/include/mm/mmmatrix.hpp
@@ -12,580 +12,117 @@
#pragma once
#include <iostream>
-#include <iomanip>
-#include <cstring>
#include <cassert>
#include <initializer_list>
#include <array>
#include <memory>
-#include "mm/mmiterator.hpp"
namespace mm {
-
- /* basic grid structure */
+ using index = std::size_t;
template<typename T, std::size_t Rows, std::size_t Cols>
class basic_matrix;
- /* basic wrapper */
-
- template<typename T, std::size_t Rows, std::size_t Cols>
- class matrix; // simple matrix format
-
/* specialisations */
-
- /* specialization of a matrix */
- template<typename T, std::size_t N>
- class vector; // by default, set a column vector
+ template<typename T, std::size_t Rows, std::size_t Cols>
+ class matrix;
template<typename T, std::size_t N>
class square_matrix;
- /* specialisation of a square_matrix for a sub-diagonal composed matrix */
- template<typename T, std::size_t N, signed long ... Diags>
- class multi_diag_matrix;
+ template<typename T, std::size_t N>
+ class diagonal_matrix;
}
/*
* Matrix class, no access methods
*/
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-class mm::basic_matrix
-{
-public:
- using type = T;
-
- template<typename U, std::size_t ORows, std::size_t OCols>
- friend class mm::basic_matrix;
-
- template<typename U, std::size_t ORows, std::size_t OCols>
- friend class mm::matrix;
-
- template<typename U, std::size_t ORows, std::size_t OCols, class Grid>
- friend class mm::iter::basic_iterator;
-
- template<typename U, std::size_t ON, class Grid>
- friend class mm::iter::diag_iterator;
-
- //template<typename U, std::size_t ORows, std::size_t OCols, class Grid>
- //friend class mm::iter::basic_iterator<T, Rows, Cols, mm::basic_matrix<T, Rows, Cols>>;
-
- //template<typename U, std::size_t ORows, std::size_t OCols, class Grid>
- //friend class mm::iter::basic_iterator<typename std::add_const<T>::type, Rows, Cols, typename std::add_const<mm::basic_matrix<T, Rows, Cols>>::type>;
-
- basic_matrix();
-
- // from initializer_list
- basic_matrix(std::initializer_list<std::initializer_list<T>> l);
-
- // copyable and movable
- basic_matrix(const basic_matrix<T, Rows, Cols>& other) = default;
- basic_matrix(basic_matrix<T, Rows, Cols>&& other) = default;
-
- // copy from another matrix
- template<std::size_t ORows, std::size_t OCols>
- basic_matrix(const basic_matrix<T, ORows, OCols>& other);
-
- void swap_rows(std::size_t x, std::size_t y);
- void swap_cols(std::size_t x, std::size_t y);
-
- // mathematical operations
- //virtual basic_matrix<T, Cols, Rows> transposed() const;
- //inline basic_matrix<T, Cols, Rows> td() const { return transposed(); }
-
-protected:
- template<typename ConstIterator>
- basic_matrix(ConstIterator begin, ConstIterator end);
-
-private:
- std::array<T, Rows * Cols> data;
-};
-
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::basic_matrix<T, Rows, Cols>::basic_matrix() {
- std::fill(data.begin(), data.end(), 0);
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::basic_matrix<T, Rows, Cols>::basic_matrix(
- std::initializer_list<std::initializer_list<T>> l
-) {
- assert(l.size() == Rows);
- auto data_it = data.begin();
-
- for (auto&& row : l) {
- data_it = std::copy(row.begin(), row.end(), data_it);
- }
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-template<std::size_t ORows, std::size_t OCols>
-mm::basic_matrix<T, Rows, Cols>::basic_matrix(
- const mm::basic_matrix<T, ORows, OCols>& other
-) {
- static_assert((ORows <= Rows),
- "cannot copy a taller matrix into a smaller one"
- );
-
- static_assert((OCols <= Cols),
- "cannot copy a larger matrix into a smaller one"
- );
-
- std::fill(data.begin(), data.end(), 0);
- for (unsigned row = 0; row < Rows; row++)
- for (unsigned col = 0; col < Cols; col++)
- this->at(row, col) = other.at(row, col);
-}
-
-/* protected construtor */
-template<typename T, std::size_t Rows, std::size_t Cols>
-template<typename ConstIterator>
-mm::basic_matrix<T, Rows, Cols>::basic_matrix(
- ConstIterator begin, ConstIterator end
-) {
- assert(static_cast<unsigned>(std::distance(begin, end)) >= ((Rows * Cols)));
- std::copy(begin, end, data.begin());
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-void mm::basic_matrix<T, Rows, Cols>::swap_rows(std::size_t x, std::size_t y) {
- if (x == y)
- return;
-
- for (unsigned col = 0; col < Cols; col++)
- std::swap(this->at(x, col), this->at(y, col));
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-void mm::basic_matrix<T, Rows, Cols>::swap_cols(std::size_t x, std::size_t y) {
- if (x == y)
- return;
-
- for (unsigned row = 0; row < Rows; row++)
- std::swap(this->at(row, x), this->at(row, y));
-}
-
-/*
- * Matrix object
- */
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-class mm::matrix
-{
-public:
-
- //template<typename U, std::size_t ORows, std::size_t OCols>
- using vec_iterator = mm::iter::basic_iterator<T, Rows, Cols, mm::basic_matrix<T, Rows, Cols>>;
-
- //template<typename U, std::size_t ORows, std::size_t OCols>
- using const_vec_iterator = mm::iter::basic_iterator<typename std::add_const<T>::type, Rows, Cols, typename std::add_const<mm::basic_matrix<T, Rows, Cols>>::type>;
-
- // default zeros constructor
- matrix() : M(std::make_shared<mm::basic_matrix<T, Rows, Cols>>()), transposed(false) {}
-
- // from initializer_list
- matrix(std::initializer_list<std::initializer_list<T>> l)
- : M(std::make_shared<mm::basic_matrix<T, Rows, Cols>>(l)), transposed(false) {}
-
- // copyable and movable
- matrix(const matrix<T, Rows, Cols>& other) // deep copy
- : M(std::make_shared<mm::basic_matrix<T, Rows, Cols>>(*other.M)), transposed(other.transposed) {}
-
- matrix(basic_matrix<T, Rows, Cols>&& other) // move ptr
- : M(other.M), transposed(other.transposed)
- {
- other.M = nullptr;
- }
-
- matrix<T, Rows, Cols> operator=(const basic_matrix<T, Rows, Cols>& other) // deep copy
- {
- *M = *other.M;
- transposed = other.transposed;
- }
-
- /*
- * Transposition
- */
-
- matrix<T, Rows, Cols>& transpose_d()
- {
- transposed = !transposed;
- return *this;
- }
-
- const matrix<T, Rows, Cols> transpose() const
- {
- return matrix<T, Rows, Cols>(M, !transposed);
- }
-
- inline matrix<T, Rows, Cols>& td()
- {
- return transpose();
- }
-
- inline matrix<T, Rows, Cols> t() const
- {
- return transpose();
- }
-
- // strongly transpose
- matrix<T, Cols, Rows> transpose_cpy() const
- {
- matrix<T, Cols, Rows> out(); // copy
- // TODO
- }
-
- /*
- * Pointer status
- */
-
- bool expired() const
+namespace mm {
+ template<typename T, std::size_t Rows, std::size_t Cols>
+ class basic_matrix
{
- return M == nullptr;
- }
-
- /*
- * Downcasting conditions
- */
-
- /// downcast to square matrix
- static inline constexpr bool is_square() { return (Rows == Cols); }
- inline constexpr square_matrix<T, Rows> to_square() const {
- static_assert(is_square());
- return static_cast<square_matrix<T, Rows>>(*this);
- }
-
- /// downcast to col_vector
- static inline constexpr bool is_vector() { return (Rows == 1 || Cols == 1); }
- inline vector<T, Cols> to_vector() const {
- if constexpr(Cols == 1)
- return static_cast<vector<T, Rows>>(*this);
- else if (Rows == 1)
- return vector<T, Cols>(*this); // copy into column vector
+ public:
+ using type = T;
- }
+ template<typename U, std::size_t ORows, std::size_t OCols>
+ friend class mm::matrix;
- /* Accessors */
+ // copy from another matrix
+ template<std::size_t ORows, std::size_t OCols>
+ matrix(const basic_matrix<T, ORows, OCols>& other);
- virtual T& at(std::size_t row, std::size_t col)
- {
- return (transposed) ? M->data[col * Cols + row] : M->data[row * Cols + col];
- }
+ virtual T& at(index row, index col) = 0;
+ virtual const T& at(index row, index col) const = 0;
+ };
- virtual const T& at(std::size_t row, std::size_t col) const
- {
- return (transposed) ? M->data[col * Cols + row] : M->data[row * Cols + col];
- }
- std::size_t rows() const {
- return (transposed) ? Cols : Rows;
- }
-
- std::size_t cols() const {
- return (transposed) ? Rows : Cols;
- }
-
- virtual mm::matrix<T, Rows, Cols>::vec_iterator operator[](std::size_t index)
- {
- return mm::matrix<T, Rows, Cols>::vec_iterator(*M, index, 0, !transposed);
- }
+ /* Specializations */
- virtual mm::matrix<T, Rows, Cols>::const_vec_iterator operator[](std::size_t index) const
+ template<typename T, std::size_t Rows, std::size_t Cols>
+ struct matrix : public basic_matrix<T, N>
{
- return mm::matrix<T, Rows, Cols>::const_vec_iterator(*M, index, 0, !transposed);
- }
-
- /*
- * Basic matematical operations (dimension indipendent)
- */
-
- mm::matrix<T, Rows, Cols>& operator+=(const mm::matrix<T, Rows, Cols>& m) {
-
- for (unsigned row = 0; row < std::min(rows(), m.rows()); ++row)
- for (unsigned col = 0; col < std::min(cols(), m.cols()); ++col)
- at(row, col) += m.at(row, col);
-
- return *this;
- }
-
- mm::matrix<T, Rows, Cols>& operator-=(const mm::matrix<T, Rows, Cols>& m) {
-
- for (unsigned row = 0; row < std::min(rows(), m.rows()); ++row)
- for (unsigned col = 0; col < std::min(cols(), m.cols()); ++col)
- at(row, col) -= m.at(row, col);
-
- return *this;
- }
-
- mm::matrix<T, Rows, Cols> operator*=(const T& k) {
-
- for (unsigned row = 0; row < rows(); ++row)
- for (auto& x : (*this)[row])
- x *= k;
-
- return *this;
- }
-
-protected:
-
- std::shared_ptr<mm::basic_matrix<T, Rows, Cols>> M;
-
- // shallow construction
- matrix(std::shared_ptr<mm::basic_matrix<T, Rows, Cols>> grid, bool tr = false) : M(grid), transposed(tr) {}
-
-private:
-
- bool transposed;
-};
-
-/* Basic operator overloading (dimension indipendent) */
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::matrix<T, Rows, Cols> operator+(
- mm::matrix<T, Rows, Cols> a,
- const mm::matrix<T, Rows, Cols>& b
-) {
- return a += b;
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::matrix<T, Rows, Cols> operator-(
- mm::matrix<T, Rows, Cols> a,
- const mm::matrix<T, Rows, Cols>& b
-) {
- return a -= b;
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::matrix<T, Rows, Cols> operator*(
- mm::matrix<T, Rows, Cols> a,
- const T& k
-) {
- return a *= k;
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::matrix<T, Rows, Cols> operator*(
- const T& k,
- mm::matrix<T, Rows, Cols> a
-) {
- return a *= k;
-}
-
-// simple multiplication
-template<typename T, std::size_t M, std::size_t P1, std::size_t P2, std::size_t N>
-mm::matrix<T, M, N> operator*(
- const mm::matrix<T, M, P1>& a,
- const mm::matrix<T, P2, N>& b
-) {
- // TODO, adjust asserts for transposed cases
- static_assert(P1 == P2, "invalid matrix multiplication");
- assert(a.cols() == b.rows());
-
- mm::matrix<T, M, N> result;
- const mm::matrix<T, P2, N> bt = b.t(); // weak transposition
-
- //npdebug("Calling *")
-
- for (unsigned row = 0; row < M; row++)
- for (unsigned col = 0; col < N; col++)
- result.at(row, col) = a[row] * bt[col]; // scalar product
-
- return result;
-}
-
-/*
- * Matrix operator <<
- */
-
-template<typename T, std::size_t Rows, std::size_t Cols, unsigned NumW = 3>
-std::ostream& operator<<(std::ostream& os, const mm::matrix<T, Rows, Cols>& m) {
-
- for (unsigned index = 0; index < m.rows(); index++) {
- os << "[ ";
- for (unsigned col = 0; col < m.cols()-1; ++col) {
- os << std::setw(NumW) << m.at(index, col) << ", ";
+ public:
+ virtual T& at(index row, index col) override {
+ return m_data[row * Cols + col];
}
- os << std::setw(NumW) << m.at(index, m.cols()-1) << " ]\n";
- }
-
- return os;
-}
-
-/*
- * Vector, TODO better manage column and row
- */
-
-template<typename T, std::size_t N>
-class mm::vector : public mm::matrix<T, N, 1>
-{
-public:
-
- using mm::matrix<T, N, 1>::matrix;
-
- vector(std::initializer_list<T> l)
- : mm::matrix<T, N, 1>(l) {}
-};
-template<typename T>
-mm::vector<T, 3> operator^(const mm::vector<T, 3>& v, const mm::vector<T, 3>& w)
-{
- mm::vector<T, 3> out;
-
- out[0] = v[1] * w[2] - v[2] * w[2];
- out[1] = v[2] * w[0] - v[0] * w[2];
- out[2] = v[0] * w[1] - v[1] * w[0];
-
- return out;
-}
-
-/*
- * Square matrix
- */
-
-template<typename T, std::size_t N>
-class mm::square_matrix : public mm::matrix<T, N, N>
-{
-public:
-
- using mm::matrix<T, N, N>::matrix;
-
- using diag_iterator = mm::iter::diag_iterator<T, N, mm::basic_matrix<T, N, N>>;
-
- using const_diag_iterator = mm::iter::diag_iterator<typename std::add_const<T>::type, N, typename std::add_const<mm::basic_matrix<T, N, N>>::type>;
-
- virtual T trace();
- inline T tr() { return trace(); }
-
- virtual mm::square_matrix<T, N>::diag_iterator diag_beg(int row = 0)
- {
- return diag_iterator(*(this->M), row, 0);
- }
-
- virtual mm::square_matrix<T, N>::const_diag_iterator diag_end(int row = 0) const
- {
- return const_diag_iterator(*(this->M), row, N);
- }
-
- // TODO, determinant
-
- /// in place inverse
- // TODO, det != 0
- // TODO, use gauss jordan for invertible ones
- //void invert();, TODO, section algorithm
-
- /*
- * Generate the identity
- */
-
- static inline constexpr mm::square_matrix<T, N> identity() {
- mm::square_matrix<T, N> i;
- for (unsigned row = 0; row < N; row++)
- for (unsigned col = 0; col < N; col++)
- i.at(row, col) = (row == col) ? 1 : 0;
-
- return i;
- }
-};
-
-template<typename T, std::size_t N>
-T mm::square_matrix<T, N>::trace()
-{
- T sum = 0;
- for (const auto& x : diag_beg())
- sum += x;
+ virtual const T& at(index row, index col) const override {
+ return at(row, col);
+ }
- return sum;
-}
+ private:
+ std::array<T, Rows * Cols> m_data;
+ };
-// TODO, static assert, for all: Diags > -N, Diags < N
-// TODO, force Diags to be ordered
-template<typename T, std::size_t N, signed long ... Diags>
-class mm::multi_diag_matrix : public mm::square_matrix<T, N>
-{
- T& shared_zero = 0;
-public:
- using mm::square_matrix<T, N>::square_matrix;
+ template<typename T, std::size_t N>
+ struct vector : public matrix<T, 1, N> {};
- // TODO, ordered case: dichotomy search O(log(M))
- // M = parameter pack size
- static inline bool constexpr is_in(std::size_t i, std::size_t j)
+ template<typename T, std::size_t N>
+ struct square_matrix : public basic_matrix<T, N>
{
- auto t = std::make_tuple(Diags...);
+ public:
+ virtual T& at(index row, index col) override {
+ return m_data[row * N + col];
+ }
- for(unsigned k(0); k < sizeof...(Diags); ++k)
- if ((i - j) == std::get<k>(t))
- return true;
+ virtual const T& at(index row, index col) const override {
+ return at(row, col);
+ }
- return false;
- }
+ private:
+ std::array<T, N*N> m_data;
+ };
- virtual T& at(std::size_t row, std::size_t col) override
+ template<typename T, std::size_t N>
+ struct identity_matrix : public basic_matrix<T, N, N>
{
- if (is_in(row, col))
- return mm::square_matrix<T, N>::at(row, col);
+ public:
+ const T& at(index row, index col) const override {
+ return (row != col) ? static_cast<T>(1) : static_cast<T>(0);
+ }
- shared_zero = 0;
- return shared_zero;
+ private:
+ T m_useless;
+ T& at(index row, index col) { return m_useless; }
}
- virtual const T& at(std::size_t row, std::size_t col) const override
+ template<typename T, std::size_t N>
+ struct diagonal_matrix : public basic_matrix<T, N, N>
{
- if (is_in(row, col))
- return mm::square_matrix<T, N>::at(row, col);
-
- return 0;
- }
-
-
-
- // TODO, implement limited iterators
-};
+ public:
+ T& at(index row, index col) override {
+ n_null_element = static_cast<T>(0);
+ return (row != col) ? m_data[row] : n_null_element;
+ }
-/*template<typename T, std::size_t N, signed long ... Diags, std::size_t P, std::size_t M>
-void constexpr diag_mult(const mm::multi_diag_matrix<T, N, Diags...>& a,
- const mm::matrix<T, P, M>& b, mm::matrix<T, M, N>& result)
-{
- static_assert(N == P && N == M, "invalid diagonal multiplication");
+ const T& at(index row, index col) const override {
+ return (row != col) ? m_data[row] : static_cast<T>(0);
+ }
- auto d = a.diagonal(Diags);
- if constexpr (Diags < 0) {
- for (unsigned k = 0; k < M; ++k)
- for (unsigned i = -Diags; i < N; ++i)
- result.at(i + Diags, k) += d[i + Diags] * b.at(i, k);
- } else {
- for (unsigned k = 0; k < M; ++k)
- for (unsigned i = Diags; i < N; ++i)
- result.at(i, k) += d[i - Diags] * b.at(i - Diags, k);
+ private:
+ T m_null_element;
+ std::array<T, N> m_data;
}
}
-
-template<typename T, std::size_t N, signed long ... Diags, std::size_t P, std::size_t M>
-mm::matrix<T, M, N> operator*(
- const mm::multi_diag_matrix<T, N, Diags...>& a,
- const mm::matrix<T, P, M>& b
-) {
- static_assert(N == P && N == M, "invalid matrix multiplication");
- assert(a.cols() == b.rows());
-
- mm::matrix<T, M, N> result;
- ((
- auto d = a.diagonal(Diags);
- if constexpr (Diags < 0) {
- for (unsigned k = 0; k < M; ++k)
- for (unsigned i = -Diags; i < N; ++i)
- result.at(i + Diags, k) += d[i + Diags] * b.at(i, k);
- } else {
- for (unsigned k = 0; k < M; ++k)
- for (unsigned i = Diags; i < N; ++i)
- result.at(i, k) += d[i - Diags] * b.at(i - Diags, k);
- }
- ) ...);
-
- return result;
-}*/
-
diff --git a/include/mm/view.hpp b/include/mm/view.hpp
new file mode 100644
index 0000000..910c16a
--- /dev/null
+++ b/include/mm/view.hpp
@@ -0,0 +1,66 @@
+#pragma once
+
+#include <mmmatrix.hpp>
+
+
+namespace mm::alg {
+
+ template <
+ template<typename, std::size_t, std::size_t> typename Matrix,
+ typename T, std::size_t Rows, std::size_t Cols
+ >
+ struct visitor
+ {
+ using type = T;
+
+ // copy constructible
+ visitor(const visitor& other) = default;
+
+ T& operator()(const Matrix<T, Rows, Cols>& m, index row, index col) {
+ return m.at(row, col);
+ }
+
+ const T& operator()(const Matrix<T, Rows, Cols>& m, index row, index col) {
+ return operator()(m, row, col);
+ }
+ };
+
+ template <
+ template<typename, std::size_t, std::size_t> typename Matrix,
+ typename T, std::size_t Rows, std::size_t Cols
+ >
+ struct transpose : public visitor<Matrix, T, Rows, Cols>
+ {
+ T& operator()(const Matrix<T, Rows, Cols> m, index row, index col) {
+ // assert(col < Rows)
+ // assert(row < Cols)
+ return m.at(col, row);
+ }
+ };
+}
+
+namespace mm {
+ template <
+ template<typename, std::size_t, std::size_t> typename Matrix,
+ typename T, std::size_t Rows, std::size_t Cols
+ >
+ struct view
+ {
+ Matrix<T, Rows, Cols>& m;
+ // std::stack<std::unique_ptr<alg::visitor>> visitors;
+ std::unique_ptr<alg::visitor> visitor;
+
+ T& at(index row, index col) {
+ return visitor(m, row, col);
+ }
+
+ view& operator|=(const alg::visitor& other) {
+ // visitors.push(std::move(std::make_unique<alg::visitor>(other)));
+ visitor = std::make_unique<alg::visistor>(other);
+ }
+ };
+
+ view operator|(const view& left, const alg::visitor& right) {
+ return left |= right;
+ }
+}
diff --git a/test/matrix_example.cpp b/test/matrix_example.cpp
index a3f0eac..96ba67d 100644
--- a/test/matrix_example.cpp
+++ b/test/matrix_example.cpp
@@ -32,27 +32,5 @@ int main(int argc, char *argv[]) {
std::cout << "a.td() = \n" << a.t(); // or a.trasposed();
std::cout << std::endl;
- // square matrix
- mm::square_matrix<std::complex<int>, 2> f {{{2, 3}, {1, 4}}, {{6, 1}, {-3, 4}}};
-
- std::cout << "Square matrix" << std::endl;
- std::cout << "f = \n" << f;
-
- std::cout << "tr(f) = " << f.tr(); //or f.trace() << std::endl;
-
- auto ft = f.t();
- std::cout << "after in place transpose f.t(), f = \n" << ft;
- std::cout << std::endl;
-
- auto identity = mm::square_matrix<int, 3>::identity();
-
- std::cout << "Identity matrix" << std::endl;
- std::cout << "I = \n" << identity;
- std::cout << std::endl;
-
- // vector
-
- //
-
return 0;
}
--
cgit v1.2.1