Optimising matrices access and operations

Creating the K-diagonal matrix

3 files changed, 832 insertions, 9 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
diff --git a/include/mm/mmmatrix.hpp b/include/mm/mmmatrix.hpp
index 5285429..c11b626 100644
--- a/include/mm/mmmatrix.hpp
+++ b/include/mm/mmmatrix.hpp
@@ -19,9 +19,14 @@
 #include <array>
 
 namespace mm {
+
     template<typename T, std::size_t Rows, std::size_t Cols>
     class basic_matrix;
 
+    // TODO, not sure it's a good idea
+    //template<typename T, std::size_t Rows, std::size_t Cols>
+    //class transposed_matrix;
+
     /* specialization of basic_matrx for Cols = 1 */
     template<typename T, std::size_t Rows>
     class row_vec;
@@ -38,10 +43,251 @@ namespace mm {
     template<typename T, std::size_t N>
     class square_matrix;
 
-    // template<typename T, std::size_t N>
-    // class diag_matrix;
+    template<typename T, std::size_t N>
+    class diagonal_matrix;
+
+    /* 
+     * Iterators 
+     */
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    class vector_iterator;
+
+    template<typename T, std::size_t N>
+    class diag_iterator;
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    class const_vector_iterator;
+
+    template<typename T, std::size_t N>
+    class const_diag_iterator;
 }
 
+/* Non-const Iterators */
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::vector_iterator
+{
+    std::size_t index; // variable index
+
+    mm::basic_matrix<T, Rows, Cols>& M;
+
+    const std::size_t position; // fixed index
+    const bool direction; // true = row, false = column
+
+public:
+    template<typename U, std::size_t ORows, std::size_t OCols>
+    friend class vector_iterator;
+
+    vector_iterator(mm::basic_matrix<T, Rows, Cols>& M, std::size_t position, bool direction);
+
+    mm::vector_iterator<T, Rows, Cols> operator++()
+    {
+        vector_iterator<T, Rows, Cols> it = *this;
+        ++index;
+        return it;
+    }
+
+    mm::vector_iterator<T, Rows, Cols> operator--()
+    {
+        vector_iterator<T, Rows, Cols> it = *this;
+        --index;
+        return it;
+    }
+
+    mm::vector_iterator<T, Rows, Cols>& operator++(int)
+    {
+        ++index;
+        return *this;
+    }
+
+    mm::vector_iterator<T, Rows, Cols>& operator--(int)
+    {
+        --index;
+        return *this;
+    }
+
+    bool operator==(const mm::vector_iterator<T, Rows, Cols>& other) const
+    {
+        return index == other.index;
+    }
+
+    bool operator=!(const mm::vector_iterator<T, Rows, Cols>& other) const
+    {
+        return index != other.index;
+    }
+
+    T& operator*() const;
+    T& operator[](std::size_t);
+};
+
+template<typename T, std::size_t N>
+class diag_iterator
+{
+    std::size_t index; // variable index
+
+    mm::square_matrix<T, N>& M;
+
+    const int position; // fixed diagonal index
+
+public:
+    template<typename U, std::size_t ON>
+    friend class diag_iterator;
+
+    diag_iterator(mm::square_matrix<T, N>& M, std::size_t position, bool direction);
+
+    mm::diag_iterator<T, N> operator++()
+    {
+        diag_iterator<T, N> it = *this;
+        ++index;
+        return it;
+    }
+
+    mm::diag_iterator<T, N> operator--()
+    {
+        diag_iterator<T, N> it = *this;
+        --index;
+        return it;
+    }
+
+    mm::diag_iterator<T, N>& operator++(int)
+    {
+        ++index;
+        return *this;
+    }
+
+    mm::diag_iterator<T, N>& operator--(int)
+    {
+        --index;
+        return *this;
+    }
+
+    bool operator==(const mm::diag_iterator<T, N>& other) const
+    {
+        return index == other.index;
+    }
+
+    bool operator=!(const mm::diag_iterator<T, N>& other) const
+    {
+        return index != other.index;
+    }
+
+    T& operator*() const;
+};
+
+/* Const Iterators */
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::const_vector_iterator
+{
+    std::size_t index; // variable index
+
+    const mm::basic_matrix<T, Rows, Cols>& M;
+
+    const std::size_t position; // fixed index
+    const bool direction; // true = row, false = column
+
+public:
+    const_vector_iterator(mm::basic_matrix<T, Rows, Cols>& M, std::size_t position, bool direction);
+
+    mm::const_vector_iterator<T, Rows, Cols> operator++()
+    {
+        vector_iterator<T, Rows, Cols> it = *this;
+        ++index;
+        return it;
+    }
+
+    mm::const_vector_iterator<T, Rows, Cols> operator--()
+    {
+        vector_iterator<T, Rows, Cols> it = *this;
+        --index;
+        return it;
+    }
+
+    mm::const_vector_iterator<T, Rows, Cols>& operator++(int)
+    {
+        ++index;
+        return *this;
+    }
+
+    mm::const_vector_iterator<T, Rows, Cols>& operator--(int)
+    {
+        --index;
+        return *this;
+    }
+
+    bool operator==(const mm::const_vector_iterator<T, Rows, Cols>& other) const
+    {
+        return index == other;
+    }
+
+    bool operator=!(const mm::const_vector_iterator<T, Rows, Cols>& other) const
+    {
+        return index != other;
+    }
+
+    const T& operator*() const;
+    const T& operator[](std::size_t) const;
+};
+
+template<typename T>
+class const_diag_iterator
+{
+    std::size_t index; // variable index
+
+    const mm::square_matrix<T, N>& M;
+
+    const int position; // fixed diagonal index
+
+public:
+    template<typename U, std::size_t ON>
+    friend class const_diag_iterator;
+
+    const_diag_iterator(const mm::square_matrix<T, N>& M, std::size_t position, bool direction);
+
+    mm::const_diag_iterator<T, N> operator++()
+    {
+        const_diag_iterator<T, N> it = *this;
+        ++index;
+        return it;
+    }
+
+    mm::const_diag_iterator<T, N> operator--()
+    {
+        const_diag_iterator<T, N> it = *this;
+        --index;
+        return it;
+    }
+
+    mm::const_diag_iterator<T, N>& operator++(int)
+    {
+        ++index;
+        return *this;
+    }
+
+    mm::const_diag_iterator<T, N>& operator--(int)
+    {
+        --index;
+        return *this;
+    }
+
+    bool operator==(const mm::const_diag_iterator<T, N>& other) const
+    {
+        return index == other.index;
+    }
+
+    bool operator=!(const mm::const_diag_iterator<T, N>& other) const
+    {
+        return index != other.index;
+    }
+
+    const T& operator*() const;
+};
+
+/*
+ * Matrix class
+ */
+
 template<typename T, std::size_t Rows, std::size_t Cols>
 class mm::basic_matrix {
 public:
@@ -50,6 +296,9 @@ public:
     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::vector_iterator;
+
     static constexpr std::size_t rows = Rows;
     static constexpr std::size_t cols = Cols;
 
@@ -67,21 +316,20 @@ public:
     basic_matrix(const basic_matrix<T, ORows, OCols>& other);
 
     // access data
-    T& at(std::size_t row, std::size_t col);
-    const T& at(std::size_t row, std::size_t col) const;
+    virtual T& at(std::size_t row, std::size_t col);
+    virtual 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);
+    virtual auto operator[](std::size_t index);
 
     void swap_rows(std::size_t x, std::size_t y);
     void swap_cols(std::size_t x, std::size_t y);
 
     // mathematical operations
+    // TODO, simply switch iteration mode
     virtual basic_matrix<T, Cols, Rows> transposed() const;
     inline basic_matrix<T, Cols, Rows> td() const { return transposed(); }
 
-    // bool is_invertible() const;
-    // basic_matrix<T, Rows, Cols> inverse() const;
-
 
     /// downcast to square matrix
     static inline constexpr bool is_square() { return (Rows == Cols); }
@@ -307,6 +555,9 @@ std::ostream& operator<<(std::ostream& os, const mm::basic_matrix<T, Rows, Cols>
 
 
 
+/* 
+ * derivated classes 
+ */
 
 /* row vector specialization */
 template<typename T, std::size_t Rows>
@@ -329,7 +580,35 @@ public:
     using mm::basic_matrix<T, Rows, Cols>::basic_matrix;
 };
 
-/* square matrix specializaiton */
+/* 
+ * transposed matrix format 
+ * TODO: write this class, or put a bool flag into the original one
+ */
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::transposed_matrix : public mm::basic_matrix<T, Rows, Cols>
+{
+public:
+    using mm::basic_matrix<T, Rows, Cols>::basic_matrix;
+
+    virtual T& at(std::size_t row, std::size_t col) override
+    {
+        return mm::basic_matrix<T, Rows, Cols>::at(col, row);
+    }
+
+    virtual const T& at(std::size_t row, std::size_t col) const override
+    {
+        return mm::basic_matrix<T, Rows, Cols>::at(col, row);
+    }
+
+    // allows to access a matrix M at row j col k with M[j][k]
+    virtual auto operator[](std::size_t index) override
+    {
+        // TODO, return other direction iterator
+    }
+}
+
+/* square matrix specialization */
 template<typename T, std::size_t N>
 class mm::square_matrix : public mm::basic_matrix<T, N, N> {
 public:
@@ -343,8 +622,20 @@ public:
     inline T tr() { return trace(); }
 
     /// in place inverse
+    // TODO, det != 0
+    // TODO, use gauss jordan for invertible ones
     void invert();
 
+
+    // TODO, downcast to K-diagonal, user defined cast
+    template<int K>
+    operator mm::diagonal_matrix<T, N, K>() const
+    {
+        // it's always possible to do it bidirectionally, 
+        // without loosing information
+        return dynamic_cast<mm::diagonal_matrix<T, N, K>>(*this);
+    }
+
     // get the identity of size N
     static inline constexpr square_matrix<T, N> identity() {
         square_matrix<T, N> i;
@@ -356,6 +647,20 @@ public:
     }
 };
 
+/* 
+ * K-diagonal square matrix format 
+ * K is bounded between ]-N, N[
+ */
+
+template<typename T, std::size_t N, int K>
+class mm::diagonal_matrix : public mm::square_matrix
+{
+public:
+    using mm::square_matrix<T, N>::square_matrix;
+
+    // TODO, redefine at, operator[]
+    // TODO, matrix multiplication
+};
 
 template<typename T, std::size_t N>
 void mm::square_matrix<T, N>::transpose() {
@@ -372,3 +677,85 @@ T mm::square_matrix<T, N>::trace() {
 
     return sum;
 }
+
+/* Iterators implementations */
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::vector_iterator<T, Rows, Cols>::vector_iterator(mm::basic_matrix<T, Rows, Cols>& _M, std::size_t pos, bool dir)
+    index(0), M(_M), position(pos), direction(dir)
+{
+    assert((dir && pos < Cols) || (!dir && pos < Rows))
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+T& mm::vector_iterator<T, Rows, Cols>::operator*() const
+{
+    return (direction) ?
+            M.data[position * Cols + index] :
+            M.data[index * Cols + position];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+T& mm::vector_iterator<T, Rows, Cols>::operator[](std::size_t i)
+{
+    return (direction) ?
+            M.data[position * Cols + i] :
+            M.data[i * Cols + position];
+}
+
+template<typename T, std::size_t N>
+mm::diag_iterator<T, N>::diag_iterator(mm::square_matrix<T, N>& _M, int pos)
+    index(0), M(_M), position(pos)
+{
+    assert(abs(pos) < N) // pos bounded between ]-N, N[ 
+}
+
+template<typename T, std::size_t N>
+T& mm::diag_iterator<T, N>::operator*() const
+{
+    return (k > 0) ?
+        M.data[(index + position) * Cols + index] :
+        M.data[index * Cols + (index - position)];
+}
+
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::const_vector_iterator<T, Rows, Cols>::const_vector_iterator(const mm::basic_matrix<T, Rows, Cols>& _M, std::size_t pos, bool dir)
+    index(0), M(_M), position(pos), direction(dir)
+{
+    assert((dir && pos < Cols) || (!dir && pos < Rows))
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+const T& mm::const_vector_iterator<T, Rows, Cols>::operator*() const
+{
+    return (direction) ?
+            M.data[position * Cols + index] :
+            M.data[index * Cols + position];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+const T& mm::const_vector_iterator<T, Rows, Cols>::operator[](std::size_t i) const
+{
+    return (direction) ?
+            M.data[position * Cols + i] :
+            M.data[i * Cols + position];
+}
+
+template<typename T, std::size_t N>
+mm::const_diag_iterator<T, N>::const_diag_iterator(const mm::square_matrix<T, N>& _M, int pos)
+    index(0), M(_M), position(pos)
+{
+    assert(abs(pos) < N) // pos bounded between ]-N, N[ 
+}
+
+template<typename T, std::size_t N>
+T& mm::const_diag_iterator<T, N>::operator*() const
+{
+    return (k > 0) ?
+        M.data[(index + position) * Cols + index] :
+        M.data[index * Cols + (index - position)];
+}
+
+
+