summaryrefslogtreecommitdiffstats
path: root/include/mm/experiments
diff options
context:
space:
mode:
Diffstat (limited to 'include/mm/experiments')
-rw-r--r--include/mm/experiments/mmdiag_matrix.hpp159
-rw-r--r--include/mm/experiments/mmsubdiag.h277
2 files changed, 436 insertions, 0 deletions
diff --git a/include/mm/experiments/mmdiag_matrix.hpp b/include/mm/experiments/mmdiag_matrix.hpp
new file mode 100644
index 0000000..124e4b3
--- /dev/null
+++ b/include/mm/experiments/mmdiag_matrix.hpp
@@ -0,0 +1,159 @@
+#pragma once
+
+namespace mm {
+
+ template<typename T>
+ class diag_component;
+
+ template<typename T, std::size_t N>
+ class multi_diag_matrix;
+}
+
+/*
+ * Optimized case of square matrix
+ * It's a matrix only composed by a diagonal
+ */
+
+template<class T>
+class mm::diag_component
+{
+public:
+ virtual int dimension() const = 0;
+};
+
+template<class T, std::size_t N>
+class mm::diag_vector
+{
+public:
+
+ // TODO, define constructor
+
+ virtual int dimension() const override
+ {
+ return N;
+ }
+
+private:
+ std::array<T, N - ((Diag < 0) ? -Diag : Diag)> vector;
+};
+
+template<typename T, std::size_t N>
+class mm::multi_diag_matrix {
+public:
+ using type = T;
+
+ template<typename U, std::size_t N>
+ friend class mm::multi_diag_matrix;
+
+ multi_diag_matrix() : shared_zero(0) {}
+ ~multi_diag_matrix();
+
+ // copyable and movable
+ multi_diag_matrix(const multi_diag_matrix<T, N>& other);
+ multi_diag_matrix(multi_diag_matrix<T, N>&& other);
+
+ // copy from another matrix
+ template<std::size_t N>
+ multi_diag_matrix(const multi_diag_matrix<T, N>& other);
+
+ // standard access data
+ T& at(std::size_t row, std::size_t col);
+ const T& at(std::size_t row, std::size_t col) const;
+
+ // allows to access a matrix M at row j col k with M[j][k]
+ auto operator[](std::size_t index);
+
+ // swap two diagonals
+ void swap_diags(std::size_t k, std::size_t l);
+
+ // diagonal construction or substitution
+ template<int Diag, int K = N - ((Diag < 0) ? -Diag : Diag)>
+ void put_diag(const mm::diag_vector<T, K>& diag)
+ {
+ //static_assert((Diag <= -N) || (Diag >= N),
+ static_assert(K < 1,
+ "Diagonal number must be bounded between ]-N,N[")
+
+ auto exist = diagonals.find(Diag);
+
+ if (exist != diagonals.end())
+ // copy
+ *exists = diag;
+ else
+ // create and copy
+ diagonals.insert(new mm::diag_vector<T, K>(diag));
+ }
+
+ // mathematical operations
+ virtual multi_diag_matrix<T, N> transposed() const;
+ inline multi_diag_matrix<T, N> td() const { return transposed(); }
+
+ // multiplication rhs and lhs
+ // TODO, need super class matrix abstraction and auto return type
+
+ // A * M, TODO abstraction virtual method
+ template <std::size_t Rows>
+ basic_matrix<Rows, N> rhs_mult(const mm::basic_matrix<T, Rows, N>& A) const;
+
+ // M * A, TODO abstraction virtual method
+ template <std::size_t Cols>
+ basic_matrix<N, Cols> lhs_mult(const mm::basic_matrix<T, N, Cols>& A) const;
+
+protected:
+ template<typename ConstIterator>
+ multi_diag_matrix(ConstIterator begin, ConstIterator end);
+
+private:
+ // return an arbitrary zero in non-const mode
+ T shared_zero;
+
+ // ordered set of diagonals
+ std::unordered_map<int, mm::diag_component<T>*> diagonals;
+};
+
+template<typename T, std::size_t N>
+T& mm::multi_diag_matrix<T, N>::at(std::size_t row, std::size_t col) {
+ assert(row < N); // "out of row bound"
+ assert(col < N); // "out of column bound"
+
+ const int k = row - col;
+ auto diag = diagonals.find(k);
+ const int line = (k > 0) ? col : row;
+
+ return (diag == diagonals.end()) ? (shared_zero = 0) : (*diag)[line];
+}
+
+template<typename T, std::size_t N>
+const T& mm::multi_diag_matrix<T, N>::at(std::size_t row, std::size_t col) const {
+ assert(row < N); // "out of row bound"
+ assert(col < N); // "out of column bound"
+
+ const int k = row - col;
+ auto diag = diagonals.find(k);
+ const int line = (k > 0) ? col : row;
+
+ return (diag == diagonals.end()) ? 0 : (*diag)[line];
+}
+
+template<typename T, std::size_t N>
+auto mm::multi_diag_matrix<T, N>::operator[](std::size_t index) {
+ assert(index < N)
+
+ // TODO, single row mapping
+}
+
+template <typename T, std::size_t N, std::size_t Rows>
+mm::basic_matrix<Rows, N> mm::multi_diag_matrix<T, N>::rhs_mult(const mm::basic_matrix<T, Rows, N>& A) const
+{
+ // TODO
+}
+
+template <typename T, std::size_t N, std::size_t Cols>
+mm::basic_matrix<N, Cols> mm::multi_diag_matrix<T, N>::lhs_mult(const mm::basic_matrix<T, N, Cols>& A) const
+{
+ mm::basic_matrix<N, Cols> out;
+
+
+}
+
+
diff --git a/include/mm/experiments/mmsubdiag.h b/include/mm/experiments/mmsubdiag.h
new file mode 100644
index 0000000..2f6c93c
--- /dev/null
+++ b/include/mm/experiments/mmsubdiag.h
@@ -0,0 +1,277 @@
+#pragma once
+
+#ifndef __TRIDIAG_H__
+#define __TRIDIAG_H__
+
+template<class T>
+class tridiag
+{
+ std::vector<T> diag, lower, upper;
+
+ struct row
+ {
+ std::size_t index;
+ tridiag& ref;
+
+ T shared_zero = 0;
+
+ row(std::size_t i, tridiag& _ref) : index(i), ref(_ref) {}
+
+ T& operator[](std::size_t j)
+ {
+ switch(static_cast<long int>(j) - index)
+ {
+ case -1:
+ return ref.lower[j];
+ case 0:
+ return ref.diag[j];
+ case 1:
+ return ref.upper[index];
+ default:
+ shared_zero = 0; // assure zero
+ return shared_zero;
+ }
+ }
+ };
+
+ struct const_row
+ {
+ std::size_t index;
+ const tridiag& ref;
+
+ const T shared_zero = 0;
+
+ const_row(std::size_t i, const tridiag& _ref) : index(i), ref(_ref) {}
+
+ const T& operator[](std::size_t j) const
+ {
+ switch(static_cast<long int>(j) - index)
+ {
+ case -1:
+ return ref.lower[j];
+ case 0:
+ return ref.diag[j];
+ case 1:
+ return ref.upper[index];
+ default:
+ return shared_zero;
+ }
+ }
+ };
+
+public:
+
+ tridiag(std::size_t N = 1) : diag(N), lower(N-1), upper(N-1)
+ {
+ }
+
+ void resize(std::size_t N)
+ {
+ diag.reserve(N);
+ lower.reserve(N-1);
+ upper.reserve(N-1);
+ }
+
+ void zeros()
+ {
+ for (std::size_t i = 0; i < size()-1; ++i)
+ diag[i] = lower[i] = upper[i] = 0;
+
+ diag[size()-1] = 0;
+ }
+
+ // set diagonal to zero
+ void zeros_diag()
+ {
+ for (std::size_t i = 0; i < size(); ++i)
+ diag[i] = 0;
+ }
+
+ row operator[](std::size_t i)
+ {
+ return row(i, *this);
+ }
+
+ const const_row operator[](std::size_t i) const
+ {
+ return const_row(i, *this);
+ }
+
+ std::size_t size() const
+ {
+ return diag.size();
+ }
+
+ template<template <typename> class V>
+ V<T> solve(const V<T>& rhs)
+ {
+ V<T> solution(diag.size());
+ V<T> new_diag(diag);
+ V<T> new_rhs(rhs);
+
+ for(std::size_t i(1); i < size(); ++i)
+ {
+ T pivot = lower[i-1]/new_diag[i-1];
+ new_diag[i] -= pivot * upper[i-1];
+ new_rhs[i] -= pivot * new_rhs[i-1];
+ }
+
+ solution[size()-1] = new_rhs[size()-1] / new_diag[size()-1];
+
+ for(int i(static_cast<int>(size()-2)); i>=0; --i)
+ solution[i] = (new_rhs[i] - upper[i]*solution[i+1]) / new_diag[i];
+
+ return solution;
+ }
+
+ static constexpr tridiag id(size_t N)
+ {
+ tridiag<T> diag;
+
+ for (size_t k = 0; k < N; ++k)
+ diag[k][k] = 1;
+
+ return diag;
+ }
+
+ const std::vector<T>& d() const
+ {
+ return diag;
+ }
+
+ const std::vector<T>& down() const
+ {
+ return lower;
+ }
+
+ const std::vector<T>& up() const
+ {
+ return upper;
+ }
+};
+
+// multiplication operator overloading
+/*template<class T, template <class> class V>
+V<T> operator*(const tridiag<T>& M, const V<T>& v)
+{
+ const std::size_t N = v.size();
+
+ V<T> u(N);
+
+ u[0] = M[0][0] * v[0] + M[0][1] * v[1];
+
+ for (std::size_t k = 1; k < N-1; ++k)
+ u[k] = M[k][k-1] * v[k-1] + M[k][k] * v[k] + M[k][k+1] * v[k+1];
+
+ u[N-1] = M[N-1][N-2] * v[N-2] + M[N-1][N-1] * v[N-1];
+
+ return u;
+}*/
+
+template<class T, template <class> class V>
+V<T> operator*(const tridiag<T>& M, const V<T>& v)
+{
+ const std::size_t N = v.size();
+
+ V<T> u(N);
+
+ u[0] = M.d()[0] * v[0] + M.up()[0] * v[1];
+
+ for (std::size_t k = 1; k < N-1; ++k)
+ u[k] = M.down()[k] * v[k-1] + M.d()[k] * v[k] + M.up()[k] * v[k+1];
+
+ u[N-1] = M.down()[N-1] * v[N-2] + M.d()[N-1] * v[N-1];
+
+ return u;
+}
+
+template<class T>
+tridiag<T> operator+(tridiag<T> M, const tridiag<T>& A)
+{
+ if (M.size() != A.size())
+ return M;
+
+ const std::size_t N = M.size();
+
+ M[0][0] += A[0][0];
+ M[0][1] += A[0][1];
+
+ for (std::size_t k = 1; k < N-1; ++k)
+ for (std::size_t j = k-1; j <= k+1; ++j)
+ M[k][j] += A[k][j];
+
+ M[N-1][N-2] += A[N-1][N-2];
+ M[N-1][N-1] += A[N-1][N-1];
+
+ return M;
+}
+
+template<class T>
+tridiag<T> operator*(tridiag<T> M, const T& value)
+{
+ const std::size_t N = M.size();
+
+ M[0][0] *= value;
+ M[0][1] *= value;
+
+ for (std::size_t k = 1; k < N-1; ++k)
+ for (std::size_t j = k-1; j <= k+1; ++j)
+ M[k][j] *= value;
+
+ M[N-1][N-2] *= value;
+ M[N-1][N-1] *= value;
+
+ return M;
+}
+
+// add value * Id
+template<class T>
+tridiag<T> operator+(tridiag<T> M, const T& value)
+{
+ const std::size_t N = M.size();
+
+ for (std::size_t k = 0; k < N; ++k)
+ M[k][k] += value;
+
+ return M;
+}
+
+template<class T>
+tridiag<T> operator*(const T& value, const tridiag<T>& M)
+{
+ return M * value;
+}
+
+// add value * Id
+template<class T>
+tridiag<T> operator+(const T& value, const tridiag<T>& M)
+{
+ return M + value;
+}
+
+#include <ostream>
+
+template<class T>
+std::ostream& operator<<(std::ostream& os, const tridiag<T>& M)
+{
+ const std::size_t N = M.size();
+
+
+ os << M[0][0] << " ";
+ os << M[0][1] << " ";
+ os << std::endl;
+
+ for (std::size_t k = 1; k < N-1; ++k) {
+ for (std::size_t j = k-1; j <= k+1; ++j)
+ os << M[k][j] << " ";
+ os << std::endl;
+ }
+
+
+ os << M[N-1][N-2] << " ";
+ os << M[N-1][N-1] << std::endl;
+
+ return os;
+}
+
+#endif
generated by cgit v1.2.1 (git 2.18.0) at 2025-01-08 20:58:24 +0100