Optimized matrix section

    - Vector iterators: allow to iterate on rows, columns or diagonals
    - Transposition doesn't affect allocated space, O(1)

6 files changed, 992 insertions, 463 deletions

diff --git a/include/mm/experiments/iterators.hpp b/include/mm/experiments/iterators.hpp
index 4afd58f..fb5f568 100644
--- a/include/mm/experiments/iterators.hpp
+++ b/include/mm/experiments/iterators.hpp
@@ -619,4 +619,184 @@ T& mm::const_diag_iterator<T, N>::operator*() const
         M.data[index * Cols + (index - position)];
 }
 
+/*
+ * SECOND IMPLEMENTATION
+ *
+ */
+
+// TODO, short term solution
+#define MM_ROW_ITER 0
+#define MM_COL_ITER 1
+#define MM_DIAG_ITER 2
+
+template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
+class mm::vector_iterator
+{
+    std::size_t index; // variable index
+
+    Grid& M;
+
+    const int position; // fixed index, negative too for diagonal iterator
+
+public:
+
+    vector_iterator(Grid& M, int position, std::size_t index = 0);
+
+    operator T&()
+    {
+        return *(*this);
+    }
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid> operator++()
+    {
+        vector_iterator<T, Rows, Cols, IterType, Grid> it = *this;
+        ++index;
+        return it;
+    }
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid> operator--()
+    {
+        vector_iterator<T, Rows, Cols, IterType, Grid> it = *this;
+        --index;
+        return it;
+    }
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid>& operator++(int)
+    {
+        ++index;
+        return *this;
+    }
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid>& operator--(int)
+    {
+        --index;
+        return *this;
+    }
+
+    bool operator==(const mm::vector_iterator<T, Rows, Cols, IterType, Grid>& other) const
+    {
+        return index == other.index;
+    }
+
+    bool operator!=(const mm::vector_iterator<T, Rows, Cols, IterType, Grid>& other) const
+    {
+        return index != other.index;
+    }
+
+    bool ok() const
+    {
+        if constexpr(IterType == MM_ROW_ITER)
+            return index < Cols;
+        else 
+            return index < Rows;
+    }
+
+    T& operator*();
+    T& operator[](std::size_t);
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid> begin()
+    {
+        return mm::vector_iterator<T, Rows, Cols, IterType, Grid>(M, position, 0);
+    }
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid> end()
+    {
+        if constexpr(IterType == MM_ROW_ITER)
+            return mm::vector_iterator<T, Rows, Cols, IterType, Grid>(M, position, Cols);
+        else 
+            return mm::vector_iterator<T, Rows, Cols, IterType, Grid>(M, position, Rows);
+    }
+
+    /*
+     * Scalar product 
+     */
+
+    template<std::size_t P>
+    T operator*(const mm::vector_iterator<T, Rows, P, IterType, Grid>& v)
+    {
+        T out(0);
+
+        for (unsigned k(0); k < Rows; ++k)
+            out += (*this)[k] * v[k];
+
+        return out;
+    }
+
+    template<std::size_t P>
+    T operator*(const mm::vector_iterator<T, P, Cols, IterType, Grid>& v)
+    {
+        T out(0);
+
+        for (unsigned k(0); k < Cols; ++k)
+            out += (*this)[k] * v[k];
+
+        return out;
+    }
+};
+
+
+/* Row Iterators */
+
+namespace mm {
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    using row_iterator = vector_iterator<T, Rows, Cols, MM_ROW_ITER, mm::basic_matrix<T, Rows, Cols>>;
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    using col_iterator = vector_iterator<T, Rows, Cols, MM_COL_ITER, mm::basic_matrix<T, Rows, Cols>>;
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    using const_row_iterator = vector_iterator<typename std::add_const<T>::type, Rows, Cols, MM_ROW_ITER, typename std::add_const<mm::basic_matrix<T, Rows, Cols>>::type>;
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    using const_col_iterator = vector_iterator<typename std::add_const<T>::type, Rows, Cols, MM_COL_ITER, typename std::add_const<mm::basic_matrix<T, Rows, Cols>>::type>;
+
+    template<typename T, std::size_t N>
+    using diag_iterator = vector_iterator<T, N, N, MM_DIAG_ITER, mm::basic_matrix<T, N, N>>;
+
+    template<typename T, std::size_t N>
+    using const_diag_iterator = vector_iterator<typename std::add_const<T>::type, N, N, MM_DIAG_ITER, typename std::add_const<mm::basic_matrix<T, N, N>>::type>;
+}
+
+
+/* Iterators implementation */
+
+template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
+mm::vector_iterator<T, Rows, Cols, IterType, Grid>::vector_iterator(Grid& _M, int pos, std::size_t i)
+    : index(i), M(_M), position(pos)
+{
+    if constexpr (IterType == MM_ROW_ITER) {
+        assert(pos < Rows);
+    } else if constexpr (IterType == MM_COL_ITER) {
+        assert(pos < Cols);
+    } else if constexpr (IterType == MM_DIAG_ITER) {
+        assert(abs(pos) < Rows);
+    }
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
+T& mm::vector_iterator<T, Rows, Cols, IterType, Grid>::operator*()
+{
+    if constexpr (IterType == MM_ROW_ITER)
+        return M.data[position * Cols + index];
+    else if constexpr (IterType == MM_COL_ITER)
+        return M.data[index * Cols + position];
+    else if constexpr (IterType == MM_DIAG_ITER)
+        return (position > 0) ?
+            M.data[(index + position) * Cols + index] :
+            M.data[index * Cols + (index - position)];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
+T& mm::vector_iterator<T, Rows, Cols, IterType, Grid>::operator[](std::size_t i)
+{
+    if constexpr (IterType == MM_ROW_ITER)
+        return M.data[position * Cols + i];
+    else if constexpr (IterType == MM_COL_ITER)
+        return M.data[i * Cols + position];
+    else if constexpr (IterType == MM_DIAG_ITER)
+        return (position > 0) ?
+            M.data[(i + position) * Cols + i] :
+            M.data[i * Cols + (i - position)];
+}
 
diff --git a/include/mm/experiments/old.hpp b/include/mm/experiments/old.hpp
new file mode 100644
index 0000000..939f335
--- /dev/null
+++ b/include/mm/experiments/old.hpp
@@ -0,0 +1,58 @@
+#pragma once
+
+
+/* 
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::basic_matrix<T, Rows, Cols>::basic_matrix(
+    const mm::basic_matrix<T, Rows, Cols>& other
+) : data(other.data) {}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::basic_matrix<T, Rows, Cols>::basic_matrix(
+    mm::basic_matrix<T, Rows, Cols>&& other
+) : data(std::forward<decltype(other.data)>(other.data)) {}*/
+
+
+/* member functions */
+
+/*template<typename T, std::size_t Rows, std::size_t Cols>
+T& mm::basic_matrix<T, Rows, Cols>::at(std::size_t row, std::size_t col) {
+    assert(row < Rows); // "out of row bound"
+    assert(col < Cols); // "out of column bound"
+
+    return data[row * Cols + col];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+const T& mm::basic_matrix<T, Rows, Cols>::at(std::size_t row, std::size_t col) const {
+    assert(row < Rows); // "out of row bound"
+    assert(col < Cols); // "out of column bound"
+
+    return data[row * Cols + col];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+auto mm::basic_matrix<T, Rows, Cols>::operator[](std::size_t index) {
+    if constexpr (is_row_vec() || is_col_vec()) {
+        return data.at(index);
+    } else {
+        return mm::row_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+
+        //return row_vec<T, Rows>(
+         //   data.cbegin() + (index * Cols),
+          //  data.cbegin() + ((index + 1) * Cols) + 1
+        );
+    }
+}*/
+
+/*template<typename T, std::size_t M, std::size_t N>
+mm::basic_matrix<T, N, M> mm::basic_matrix<T, M, N>::transposed() const {
+    mm::basic_matrix<T, N, M> result;
+
+    for (unsigned row = 0; row < M; row++)
+        for (unsigned col = 0; col < N; col++)
+            result.at(col, row) = this->at(row, col);
+
+    return result;
+}*/
diff --git a/include/mm/experiments/tranpositions.hpp b/include/mm/experiments/tranpositions.hpp
new file mode 100644
index 0000000..f124a5c
--- /dev/null
+++ b/include/mm/experiments/tranpositions.hpp
@@ -0,0 +1,284 @@
+#pragma once
+
+template<typename T, std::size_t Rows, std::size_t Cols, bool Regular>
+class mm::access : virtual public mm::basic_matrix<T, Rows, Cols>
+{
+public:
+    template<typename U, std::size_t ORows, std::size_t OCols>
+    friend class mm::access;
+
+    using mm::basic_matrix<T, Rows, Cols>::basic_matrix;
+
+    access(mm::basic_matrix<T, Rows, Cols>&& m) 
+        : basic_matrix<T, Rows, Cols>(m) {}
+
+    virtual T& at(std::size_t row, std::size_t col) = 0;
+    virtual const T& at(std::size_t row, std::size_t col) const = 0;
+
+    static constexpr std::size_t rows = Regular ? Rows : Cols;
+    static constexpr std::size_t cols = Regular ? Cols : Rows;
+};
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::row_access : public mm::access<T, Rows, Cols, true>
+{
+public:
+    template<typename U, std::size_t ORows, std::size_t OCols>
+    friend class mm::row_access;
+
+    using mm::access<T, Rows, Cols, true>::access;
+
+    row_access(mm::basic_matrix<T, Rows, Cols>&& m) 
+        : access<T, Rows, Cols, true>(m) {}
+
+    virtual T& at(std::size_t row, std::size_t col) override;
+    virtual const T& at(std::size_t row, std::size_t col) const override;
+
+    row_iterator<T, Rows, Cols> operator[](std::size_t index);
+    const_row_iterator<T, Rows, Cols> operator[](std::size_t index) const;
+};
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::col_access : public mm::access<T, Rows, Cols, false>
+{
+public:
+    template<typename U, std::size_t ORows, std::size_t OCols>
+    friend class mm::access;
+
+    using mm::access<T, Rows, Cols, false>::access;
+
+    col_access(mm::basic_matrix<T, Rows, Cols>&& m) 
+        : access<T, Rows, Cols, false>(m) {}
+
+    virtual T& at(std::size_t row, std::size_t col) override;
+    virtual const T& at(std::size_t row, std::size_t col) const override;
+
+    col_iterator<T, Rows, Cols> operator[](std::size_t index);
+    const_col_iterator<T, Rows, Cols> operator[](std::size_t index) const;
+};
+
+/*
+ * Square interface
+ */
+
+template<typename T, std::size_t N>
+class mm::square_interface : virtual public mm::basic_matrix<T, N, N>
+{
+public:
+
+    template<typename U, std::size_t ON>
+    friend class mm::square_interface;
+
+    using mm::basic_matrix<T, N, N>::basic_matrix;
+
+    T trace();
+    inline T tr() { return trace(); }
+
+    mm::diag_iterator<T, N> diagonal_iter(int index = 0)
+    {
+        return mm::diag_iterator<T, N>(*this, index);
+    }
+
+    mm::diag_iterator<T, N> diagonal_iter(int index = 0) const
+    {
+        return mm::const_diag_iterator<T, N>(*this, index);
+    }
+
+    // TODO, determinant
+
+    /// in place inverse
+    // TODO, det != 0
+    // TODO, use gauss jordan for invertible ones
+    //void invert();, TODO, section algorithm
+};
+
+/* 
+ * derivated classes 
+ */
+
+// simple format matrix
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::matrix : public mm::row_access<T, Rows, Cols>
+{
+public:
+    using mm::row_access<T, Rows, Cols>::row_access;
+
+    matrix(mm::t_matrix<T, Rows, Cols>&& m)
+        : row_access<T, Rows, Cols>(std::move<mm::basic_matrix<T, Rows, Cols>>(m))
+    {
+                
+    }
+
+    operator mm::t_matrix<T, Rows, Cols>()
+    {
+        return mm::t_matrix<T, Rows, Cols>(*this); 
+    }
+
+    mm::t_matrix<T, Rows, Cols>& t()
+    {
+        return static_cast<mm::t_matrix<T, Rows, Cols>>(*this);
+    }
+
+    const mm::t_matrix<T, Rows, Cols>& t() const
+    {
+        return static_cast<mm::t_matrix<T, Rows, Cols>>(*this);
+    }
+};
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::t_matrix : public mm::col_access<T, Rows, Cols>
+{
+public:
+    using mm::col_access<T, Rows, Cols>::col_access;
+
+    t_matrix(mm::matrix<T, Rows, Cols>&& m)
+        : col_access<T, Rows, Cols>(std::move<mm::basic_matrix<T, Rows, Cols>>(m))
+    {
+                
+    }
+
+    operator mm::matrix<T, Rows, Cols>()
+    {
+        return mm::matrix<T, Rows, Cols>(*this); 
+    }
+
+    mm::matrix<T, Rows, Cols>& t()
+    {
+        return static_cast<mm::matrix<T, Rows, Cols>>(*this);
+    }
+
+    const mm::matrix<T, Rows, Cols>& t() const
+    {
+        return static_cast<mm::matrix<T, Rows, Cols>>(*this);
+    }
+};
+
+// transposed matrix
+/*template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::t_matrix : public mm::basic_matrix<T, Rows, Cols>, virtual public mm::access<T, Rows, Cols, false>
+{
+public:
+    using mm::basic_matrix<T, Rows, Cols>::basic_matrix;
+};*/
+
+/* row vector specialization */
+template<typename T, std::size_t Rows>
+class mm::row_vec : public mm::matrix<T, Rows, 1> {
+public:
+    using mm::matrix<T, Rows, 1>::matrix;
+
+    // TODO, begin, end
+};
+
+/* column vector specialization */
+template<typename T, std::size_t Cols>
+class mm::col_vec : public mm::matrix<T, 1, Cols> {
+public:
+    using mm::matrix<T, 1, Cols>::matrix;
+
+    // TODO, begin, end
+};
+
+/* square matrix specialization */
+template<typename T, std::size_t N>
+class mm::square_matrix : public mm::matrix<T, N, N>, public mm::square_interface<T, N> {
+public:
+    using mm::matrix<T, N, N>::matrix;
+
+    // get the identity of size N
+    static inline constexpr 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>
+class mm::t_square_matrix : public mm::t_matrix<T, N, N>, public mm::square_interface<T, N> {
+public:
+    using mm::t_matrix<T, N, N>::t_matrix;
+};
+
+/* 
+ * K-diagonal square matrix format 
+ * K is bounded between ]-N, N[
+ */
+
+/*template<typename T, std::size_t N, std::size_t K>
+class mm::diagonal_matrix : public mm::square_matrix<T, N>
+{
+public:
+    using mm::square_matrix<T, N>::square_matrix;
+
+    // TODO, redefine at, operator[]
+    // TODO, matrix multiplication
+};*/
+
+/*
+ * Accessors implementation
+ */
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+T& mm::row_access<T, Rows, Cols>::at(std::size_t row, std::size_t col)
+{
+    return mm::basic_matrix<T, Rows, Cols>::data[row * Cols + col];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+const T& mm::row_access<T, Rows, Cols>::at(std::size_t row, std::size_t col) const
+{
+    return mm::basic_matrix<T, Rows, Cols>::data[row * Cols + col];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+T& mm::col_access<T, Rows, Cols>::at(std::size_t row, std::size_t col)
+{
+    return mm::basic_matrix<T, Rows, Cols>::data[col * Cols + row]; // transpose
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+const T& mm::col_access<T, Rows, Cols>::at(std::size_t row, std::size_t col) const
+{
+    return mm::basic_matrix<T, Rows, Cols>::data[col * Cols + row]; // transpose
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::row_iterator<T, Rows, Cols> mm::row_access<T, Rows, Cols>::operator[](std::size_t index)
+{
+    return mm::row_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::const_row_iterator<T, Rows, Cols> mm::row_access<T, Rows, Cols>::operator[](std::size_t index) const
+{
+    return mm::const_row_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::col_iterator<T, Rows, Cols> mm::col_access<T, Rows, Cols>::operator[](std::size_t index)
+{
+    return mm::col_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::const_col_iterator<T, Rows, Cols> mm::col_access<T, Rows, Cols>::operator[](std::size_t index) const
+{
+    return mm::const_col_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+}
+
+/* Square interface implementation */
+
+template<typename T, std::size_t N>
+T mm::square_interface<T, N>::trace() 
+{
+    T sum = 0;
+    for (const auto& x : diagonal_iter())
+        sum += x;
+
+    return sum;
+}
+
+
diff --git a/include/mm/mmiterator.hpp b/include/mm/mmiterator.hpp
new file mode 100644
index 0000000..c67b92b
--- /dev/null
+++ b/include/mm/mmiterator.hpp
@@ -0,0 +1,199 @@
+#pragma once
+
+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) {}
+
+    operator T&()
+    {
+        return *(*this);
+    }
+
+    IterType operator++()
+    {
+        IterType it = *this;
+        ++index;
+        return it;
+    }
+
+    IterType operator--()
+    {
+        IterType it = *this;
+        --index;
+        return it;
+    }
+
+    IterType& operator++(int)
+    {
+        ++index;
+        return *this;
+    }
+
+    IterType& operator--(int)
+    {
+        --index;
+        return *this;
+    }
+
+    bool operator==(const IterType& other) const
+    {
+        return index == other.index;
+    }
+
+    bool operator!=(const IterType& other) const
+    {
+        return index != other.index;
+    }
+
+    virtual bool ok() const = 0;
+    
+    virtual T& operator*() = 0;
+    virtual T& operator[](std::size_t) = 0;
+
+    IterType begin()
+    {
+        return IterType(M, position, 0);
+    }
+
+    virtual IterType end() = 0;
+    
+    /*
+     * Scalar product 
+     */
+
+    template<std::size_t P>
+    T operator*(const mm::iter::vector_iterator<T, Rows, P, IterType, Grid>& v)
+    {
+        T out(0);
+
+        for (unsigned k(0); k < Rows; ++k)
+            out += (*this)[k] * v[k];
+
+        return out;
+    }
+
+    template<std::size_t P>
+    T operator*(const mm::iter::vector_iterator<T, P, Cols, IterType, Grid>& v)
+    {
+        T out(0);
+
+        for (unsigned k(0); k < Cols; ++k)
+            out += (*this)[k] * v[k];
+
+        return out;
+    }
+
+protected:
+
+    Grid& M; // grid mapping
+
+    const std::size_t position; // fixed index, negative too for diagonal iterator
+    std::size_t index; // variable index
+};
+
+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;
+
+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)
+    {
+        if (direction)
+           assert(pos < Rows);
+        else
+           assert(pos < Cols);
+    }
+
+    virtual bool ok() const override
+    {
+        return (direction) ? this->index < Cols : this->index < 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 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;
+
+public:
+
+    diag_iterator(Grid& A, signed long 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>(abs(pos)), _index), sign(pos >= 0)
+    {
+        assert(this->position < N);
+    }
+
+    virtual bool ok() const 
+    {
+        return this->index < N;
+    }
+
+    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 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 9910273..f88e90c 100644
--- a/include/mm/mmmatrix.hpp
+++ b/include/mm/mmmatrix.hpp
@@ -17,39 +17,20 @@
 #include <cassert>
 #include <initializer_list>
 #include <array>
+#include <memory>
+
+#include "mm/mmiterator.hpp"
 
 namespace mm {
 
     template<typename T, std::size_t Rows, std::size_t Cols>
     class basic_matrix;
 
-    /* specialization of basic_matrix for Rows == Cols */
-    template<typename T, std::size_t N>
-    class square_interface;
-
-    /*
-     * Most general type of matrix iterator
-     * IterType = Row, Column, Diagonal
-     * Grid = constness of mm::basic_matrix
-     */
-    template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
-    class vector_iterator;
-
-    /*
-     * Access methods interface
-     */
-
-    template<typename T, std::size_t Rows, std::size_t Cols, bool Regular>
-    class access;
-
     /* specialisations */
 
     template<typename T, std::size_t Rows, std::size_t Cols>
     class matrix; // simple matrix format
 
-    template<typename T, std::size_t Rows, std::size_t Cols>
-    class t_matrix; // transposed matrix format
-
     /* specialization of basic_matrx for Cols = 1 */
     template<typename T, std::size_t Rows>
     class row_vec;
@@ -61,155 +42,20 @@ namespace mm {
     template<typename T, std::size_t N>
     class square_matrix;
 
-    template<typename T, std::size_t N>
-    class t_square_matrix;
-
     /* specialisation of a square_matrix for a sub-diagonal composed matrix */
     template<typename T, std::size_t N, std::size_t K = 0>
     class diag_matrix;
-
-    template<typename T, std::size_t N, std::size_t K = 0>
-    class t_diag_matrix;
 }
 
-// TODO, short term solution
-#define MM_ROW_ITER 0
-#define MM_COL_ITER 1
-#define MM_DIAG_ITER 2
-
-template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
-class mm::vector_iterator
-{
-    std::size_t index; // variable index
-
-    Grid& M;
-
-    const int position; // fixed index, negative too for diagonal iterator
-
-public:
-    template<typename U, std::size_t ORows, std::size_t OCols, class OIterType, class OGrid>
-    friend class vector_iterator;
-
-    vector_iterator(Grid& M, int position, std::size_t index = 0);
-
-    mm::vector_iterator<T, Rows, Cols, IterType, Grid> operator++()
-    {
-        vector_iterator<T, Rows, Cols, IterType, Grid> it = *this;
-        ++index;
-        return it;
-    }
-
-    mm::vector_iterator<T, Rows, Cols, IterType, Grid> operator--()
-    {
-        vector_iterator<T, Rows, Cols, IterType, Grid> it = *this;
-        --index;
-        return it;
-    }
-
-    mm::vector_iterator<T, Rows, Cols, IterType, Grid>& operator++(int)
-    {
-        ++index;
-        return *this;
-    }
-
-    mm::vector_iterator<T, Rows, Cols, IterType, Grid>& operator--(int)
-    {
-        --index;
-        return *this;
-    }
-
-    bool operator==(const mm::vector_iterator<T, Rows, Cols, IterType, Grid>& other) const
-    {
-        return index == other.index;
-    }
-
-    bool operator!=(const mm::vector_iterator<T, Rows, Cols, IterType, Grid>& other) const
-    {
-        return index != other.index;
-    }
-
-    bool ok() const
-    {
-        if constexpr(IterType == MM_ROW_ITER)
-            return index < Cols;
-        else 
-            return index < Rows;
-    }
 
-    T& operator*() const;
-    T& operator[](std::size_t);
-};
-
-/* Row Iterators */
-
-namespace mm {
-
-    template<typename T, std::size_t Rows, std::size_t Cols>
-    using row_iterator = vector_iterator<T, Rows, Cols, MM_ROW_ITER, mm::basic_matrix<T, Rows, Cols>>;
-
-    template<typename T, std::size_t Rows, std::size_t Cols>
-    using col_iterator = vector_iterator<T, Rows, Cols, MM_COL_ITER, mm::basic_matrix<T, Rows, Cols>>;
-
-    template<typename T, std::size_t Rows, std::size_t Cols>
-    using const_row_iterator = vector_iterator<T, Rows, Cols, MM_ROW_ITER, std::add_const<mm::basic_matrix<T, Rows, Cols>>>;
-
-    template<typename T, std::size_t Rows, std::size_t Cols>
-    using const_col_iterator = vector_iterator<T, Rows, Cols, MM_COL_ITER, std::add_const<mm::basic_matrix<T, Rows, Cols>>>;
+/*namespace mm {
 
     template<typename T, std::size_t N>
     using diag_iterator = vector_iterator<T, N, N, MM_DIAG_ITER, mm::basic_matrix<T, N, N>>;
 
     template<typename T, std::size_t N>
-    using const_diag_iterator = vector_iterator<T, N, N, MM_DIAG_ITER, std::add_const<mm::basic_matrix<T, N, N>>>;
-}
-
-/*
- * Accessors
- */
-template<typename T, std::size_t Rows, std::size_t Cols, bool Regular>
-class mm::access
-{
-public:
-    
-    //access(mm::basic_matrix<T, Rows, Cols>& ref) : M(ref) {}
-
-    T& at(std::size_t row, std::size_t col);
-    const T& at(std::size_t row, std::size_t col) const;
-
-    auto operator[](std::size_t index);
-    auto operator[](std::size_t index) const;
-
-//private:
-//    mm::basic_matrix<T, Rows, Cols>& M;
-protected:
-    std::array<T, Rows * Cols> data;
-};
-
-/*
- * Square interface
- */
-
-template<typename T, std::size_t N>
-class mm::square_interface {
-public:
-
-    //square_interface(mm:basic_matrix<T, N, N>& _M) : M(_M) {}
-
-    T trace();
-    inline T tr() { return trace(); }
-
-    // TODO, determinant
-
-    /// in place inverse
-    // TODO, det != 0
-    // TODO, use gauss jordan for invertible ones
-    //void invert();, TODO, section algorithm
-
-//private:
-//    mm:basic_matrix<T, N, N>& M; // one information more than mm::matrix !
-protected:
-    std::array<T, N * N> data;
-};
+    using const_diag_iterator = vector_iterator<typename std::add_const<T>::type, N, N, MM_DIAG_ITER, typename std::add_const<mm::basic_matrix<T, N, N>>::type>;
+}*/
 
 
 /*
@@ -217,15 +63,25 @@ protected:
  */
 
 template<typename T, std::size_t Rows, std::size_t Cols>
-class mm::basic_matrix {
+class mm::basic_matrix
+{
 public:
     using type = T;
 
     template<typename U, std::size_t ORows, std::size_t OCols>
     friend class mm::basic_matrix;
 
-    static constexpr std::size_t rows = Rows;
-    static constexpr std::size_t cols = Cols;
+    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 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();
 
@@ -233,20 +89,13 @@ public:
     basic_matrix(std::initializer_list<std::initializer_list<T>> l);
 
     // copyable and movable
-    basic_matrix(const basic_matrix<T, Rows, Cols>& other);
-    basic_matrix(basic_matrix<T, Rows, Cols>&& other);
+    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);
 
-    // access data, basic definition
-    //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);
-
     void swap_rows(std::size_t x, std::size_t y);
     void swap_cols(std::size_t x, std::size_t y);
 
@@ -254,35 +103,15 @@ public:
     //virtual basic_matrix<T, Cols, Rows> transposed() const;
     //inline basic_matrix<T, Cols, Rows> td() const { return transposed(); }
 
-    /// 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 row_vector
-    static inline constexpr bool is_row_vec() { return (Cols == 1); }
-    inline constexpr row_vec<T, Rows> to_row_vec() const {
-        static_assert(is_row_vec());
-        return static_cast<row_vec<T, Rows>>(*this);
-    }
-
-    /// downcast to col_vector
-    static inline constexpr bool is_col_vec() { return (Rows == 1); }
-    inline constexpr col_vec<T, Cols> to_col_vec() const {
-        static_assert(is_col_vec());
-        return static_cast<col_vec<T, Cols>>(*this);
-    }
-
 protected:
     template<typename ConstIterator>
     basic_matrix(ConstIterator begin, ConstIterator end);
 
-//private:
+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);
@@ -301,16 +130,6 @@ mm::basic_matrix<T, Rows, Cols>::basic_matrix(
 }
 
 template<typename T, std::size_t Rows, std::size_t Cols>
-mm::basic_matrix<T, Rows, Cols>::basic_matrix(
-    const mm::basic_matrix<T, Rows, Cols>& other
-) : data(other.data) {}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::basic_matrix<T, Rows, Cols>::basic_matrix(
-    mm::basic_matrix<T, Rows, Cols>&& other
-) : data(std::forward<decltype(other.data)>(other.data)) {}
-
-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
@@ -339,40 +158,6 @@ mm::basic_matrix<T, Rows, Cols>::basic_matrix(
     std::copy(begin, end, data.begin());
 }
 
-
-/* member functions */
-
-/*template<typename T, std::size_t Rows, std::size_t Cols>
-T& mm::basic_matrix<T, Rows, Cols>::at(std::size_t row, std::size_t col) {
-    assert(row < Rows); // "out of row bound"
-    assert(col < Cols); // "out of column bound"
-
-    return data[row * Cols + col];
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-const T& mm::basic_matrix<T, Rows, Cols>::at(std::size_t row, std::size_t col) const {
-    assert(row < Rows); // "out of row bound"
-    assert(col < Cols); // "out of column bound"
-
-    return data[row * Cols + col];
-}
-
-template<typename T, std::size_t Rows, std::size_t Cols>
-auto mm::basic_matrix<T, Rows, Cols>::operator[](std::size_t index) {
-    if constexpr (is_row_vec() || is_col_vec()) {
-        return data.at(index);
-    } else {
-        return mm::row_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
-
-        //return row_vec<T, Rows>(
-         //   data.cbegin() + (index * Cols),
-          //  data.cbegin() + ((index + 1) * Cols) + 1
-        );
-    }
-}*/
-
-
 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)
@@ -387,274 +172,297 @@ 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++)
+    for (unsigned row = 0; row < Rows; row++)
         std::swap(this->at(row, x), this->at(row, y));
 }
 
-/*template<typename T, std::size_t M, std::size_t N>
-mm::basic_matrix<T, N, M> mm::basic_matrix<T, M, N>::transposed() const {
-    mm::basic_matrix<T, N, M> result;
+/*
+ * Matrix object
+ */
 
-    for (unsigned row = 0; row < M; row++)
-        for (unsigned col = 0; col < N; col++)
-            result.at(col, row) = this->at(row, col);
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::matrix
+{
+protected:
 
-    return result;
-}*/
+    // shallow construction
+    matrix(std::shared_ptr<mm::basic_matrix<T, Rows, Cols>> grid = nullptr, bool tr = false) : M(grid), transposed(tr) {}
 
+public:
 
-/* TODO, operator overloading */
-/*template<typename T, std::size_t Rows, std::size_t Cols>
-mm::basic_matrix<T, Rows, Cols> operator+(
-    const mm::basic_matrix<T, Rows, Cols>& a,
-    const mm::basic_matrix<T, Rows, Cols>& b
-) {
-    mm::basic_matrix<T, Rows, Cols> result;
+    //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>>;
 
-    for (unsigned row = 0; row < Rows; row++)
-        for (unsigned col = 0; col < Cols; col++)
-            result.at(row, col) = a.at(row, col) + b.at(row, col);
+    //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>;
     
-    return result;
-}
+    // 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) {}
 
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::basic_matrix<T, Rows, Cols> operator*(
-    const mm::basic_matrix<T, Rows, Cols>& m,
-    const T& scalar
-) {
-    mm::basic_matrix<T, Rows, Cols> result;
-    for (unsigned row = 0; row < Rows; row++)
-        for (unsigned col = 0; col < Cols; col++)
-            result.at(row, col) = m.at(row, col) * scalar;
+    // 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) {}
 
-    return result;
-}
+    matrix(basic_matrix<T, Rows, Cols>&& other) // move ptr
+        : M(other.M), transposed(other.transposed)
+    {
+        other.M = nullptr;
+    }
 
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::basic_matrix<T, Rows, Cols> operator*(
-    const T& scalar,
-    const mm::basic_matrix<T, Rows, Cols>& m
-) {
-    return m * scalar;
-}
+    // copy from another matrix
+    /*template<std::size_t ORows, std::size_t OCols>
+    matrix(const matrix<T, ORows, OCols>& other)
+        : M(std::make_shared<mm::basic_matrix<T, Rows, Cols>(*other.M)), transposed(other.transposed) {} */
+
+    matrix<T, Rows, Cols> operator=(const basic_matrix<T, Rows, Cols>& other) // deep copy
+    {         
+        *M = *other.M;
+        transposed = other.transposed;
+    }
 
-template<typename T, std::size_t M, std::size_t P1, std::size_t P2, std::size_t N>
-mm::basic_matrix<T, M, N> operator*(
-    const mm::basic_matrix<T, M, P1>& a,
-    const mm::basic_matrix<T, P2, N>& b
-) {
-    static_assert(P1 == P2, "invalid matrix multiplication");
-    mm::basic_matrix<T, M, N>  result;
+    /* 
+     * Transposition
+     */
 
-    // TODO: use a more efficient algorithm
-    for (unsigned row = 0; row < M; row++)
-        for (unsigned col = 0; col < N; col++)
-            for (unsigned k = 0; k < P1; k++)
-                result.at(row, col) = a.at(row, k) * b.at(k, col);
+    matrix<T, Rows, Cols>& transpose()
+    {
+        transposed = !transposed;
+        return *this;
+    }
 
-    return result;
-}
+    matrix<T, Rows, Cols> transpose() const
+    {
+        auto m = shallow_cpy();
+        m.transposed = !transposed;
 
+        return m;
+    }
 
-template<typename T, std::size_t Rows, std::size_t Cols>
-mm::basic_matrix<T, Rows, Cols> operator-(
-    const mm::basic_matrix<T, Rows, Cols>& a,
-    const mm::basic_matrix<T, Rows, Cols>& b
-) {
-    return a + (static_cast<T>(-1) * b);
-}
+    inline matrix<T, Rows, Cols>& t()
+    {
+        return transpose();
+    }
 
-template<typename T, std::size_t Rows, std::size_t Cols, unsigned NumW = 3>
-std::ostream& operator<<(std::ostream& os, const mm::basic_matrix<T, Rows, Cols>& m) {
-    for (unsigned row = 0; row < Rows; row++) {
-        os << "[ ";
-        for (unsigned col = 0; col < (Cols -1); col++) {
-            os << std::setw(NumW) << m.at(row, col) << ", ";
-        }
-        os << std::setw(NumW) << m.at(row, (Cols -1)) << " ]\n";
+    inline matrix<T, Rows, Cols> t() const
+    {
+        return transpose();
     }
 
-    return os;
-}*/
+    // strongly transpose
+    matrix<T, Cols, Rows> transpose_cpy() const
+    {
+        matrix<T, Cols, Rows> out(); // copy
+        // TODO
+    }
 
-/* 
- * derivated classes 
- */
+    /*
+     * Pointer status 
+     */
 
-// simple format matrix
-template<typename T, std::size_t Rows, std::size_t Cols>
-class mm::matrix : public mm::basic_matrix<T, Rows, Cols>, virtual public mm::access<T, Rows, Cols, true>
-{
-public:
-    using mm::basic_matrix<T, Rows, Cols>::basic_matrix;
-};
+    bool expired() const
+    {
+        return M == nullptr;
+    }
 
-// transposed matrix
-template<typename T, std::size_t Rows, std::size_t Cols>
-class mm::t_matrix : public mm::basic_matrix<T, Rows, Cols>, virtual public mm::access<T, Rows, Cols, false>
-{
-public:
-    using mm::basic_matrix<T, Rows, Cols>::basic_matrix;
-};
+    /*
+     * Downcasting conditions
+     */
 
-/* row vector specialization */
-template<typename T, std::size_t Rows>
-class mm::row_vec : public mm::matrix<T, Rows, 1> {
-public:
-    using mm::matrix<T, Rows, 1>::matrix;
+    /// 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);
+    }
 
-    // TODO, begin, end
-};
+    /// downcast to row_vector
+    static inline constexpr bool is_row_vec() { return (Cols == 1); }
+    inline constexpr row_vec<T, Rows> to_row_vec() const {
+        static_assert(is_row_vec());
+        return static_cast<row_vec<T, Rows>>(*this);
+    }
 
-/* column vector specialization */
-template<typename T, std::size_t Cols>
-class mm::col_vec : public mm::t_matrix<T, 1, Cols> {
-public:
-    using mm::t_matrix<T, 1, Cols>::t_matrix;
+    /// downcast to col_vector
+    static inline constexpr bool is_col_vec() { return (Rows == 1); }
+    inline constexpr col_vec<T, Cols> to_col_vec() const {
+        static_assert(is_col_vec());
+        return static_cast<col_vec<T, Cols>>(*this);
+    }
 
-    // TODO, begin, end
-};
+    /* Accessors */
 
-/* square matrix specialization */
-template<typename T, std::size_t N>
-class mm::square_matrix : public mm::matrix<T, N, N> , virtual public mm::square_interface<T, N> {
-public:
-    using mm::matrix<T, N, N>::matrix;
+    T& at(std::size_t row, std::size_t col)
+    {
+        return (transposed) ? M->data[col * Cols + row] : M->data[row * Cols + col];
+    }
 
-    // get the identity of size N
-    static inline constexpr 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;
+    const T& at(std::size_t row, std::size_t col) const
+    {
+        return (transposed) ? M->data[col * Cols + row] : M->data[row * Cols + col];
+    }
 
-        return i;
+    std::size_t rows() const {
+        return (transposed) ? Cols : Rows;
     }
-};
 
-template<typename T, std::size_t N>
-class mm::t_square_matrix : virtual public mm::t_matrix<T, N, N>, virtual public mm::square_interface<T, N> {
-public:
-    using mm::t_matrix<T, N, N>::t_matrix;
-};
+    std::size_t cols() const {
+        return (transposed) ? Rows : Cols;
+    }
 
-/* 
- * K-diagonal square matrix format 
- * K is bounded between ]-N, N[
- */
+    mm::matrix<T, Rows, Cols>::vec_iterator operator[](std::size_t index)
+    {
+        return mm::matrix<T, Rows, Cols>::vec_iterator(*M, index, !transposed);
+    }
 
-/*template<typename T, std::size_t N, std::size_t K>
-class mm::diagonal_matrix : public mm::square_matrix<T, N>
-{
-public:
-    using mm::square_matrix<T, N>::square_matrix;
+    mm::matrix<T, Rows, Cols>::const_vec_iterator operator[](std::size_t index) const
+    {
+        return mm::matrix<T, Rows, Cols>::const_vec_iterator(*M, index, !transposed);
+    }
+
+    /*
+     * Basic matematical operations (dimension indipendent)
+     */
 
-    // TODO, redefine at, operator[]
-    // TODO, matrix multiplication
-};*/
+    mm::matrix<T, Rows, Cols>& operator+=(const mm::matrix<T, Rows, Cols>& m) {
 
-/*template<typename T, std::size_t N>
-void mm::square_matrix<T, N>::transpose() {
-    for (unsigned row = 0; row < N; row++)
-        for (unsigned col = 0; col < row; col++)
-            std::swap(this->at(row, col), this->at(col, row));
-}*/
+        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;
+    }
 
-/* Iterators implementation */
+    mm::matrix<T, Rows, Cols>& operator-=(const mm::matrix<T, Rows, Cols>& m) {
 
-template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
-mm::vector_iterator<T, Rows, Cols, IterType, Grid>::vector_iterator(Grid& _M, int pos, std::size_t i)
-    : index(i), M(_M), position(pos)
-{
-    if constexpr (IterType == MM_ROW_ITER) {
-        assert(pos < Cols);
-    } else if constexpr (IterType == MM_COL_ITER) {
-        assert(pos < Rows);
-    } else if constexpr (IterType == MM_DIAG_ITER) {
-        assert(abs(pos) < Rows);
+        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;
     }
-}
 
-template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
-T& mm::vector_iterator<T, Rows, Cols, IterType, Grid>::operator*() const
-{
-    if constexpr (IterType == MM_ROW_ITER)
-        return M.data[position * Cols + index];
-    else if constexpr (IterType == MM_COL_ITER)
-        return M.data[index * Cols + position];
-    else if constexpr (IterType == MM_DIAG_ITER)
-        return (position > 0) ?
-            M.data[(index + position) * Cols + index] :
-            M.data[index * Cols + (index - position)];
-}
+    mm::matrix<T, Rows, Cols> operator*=(const T& k) {
 
-template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
-T& mm::vector_iterator<T, Rows, Cols, IterType, Grid>::operator[](std::size_t i)
-{
-    if constexpr (IterType == MM_ROW_ITER)
-        return M.data[position * Cols + i];
-    else if constexpr (IterType == MM_COL_ITER)
-        return M.data[i * Cols + position];
-    else if constexpr (IterType == MM_DIAG_ITER)
-        return (position > 0) ?
-            M.data[(i + position) * Cols + i] :
-            M.data[i * Cols + (i - position)];
-}
+        for (unsigned row = 0; row < rows(); ++row)
+            for (auto& x : (*this)[row])
+                x *= k;
+    
+        return *this;
+    }
 
-/*
- * Accessors implementation
- */
+protected:
 
-template<typename T, std::size_t Rows, std::size_t Cols, bool Regular>
-T& mm::access<T, Rows, Cols, Regular>::at(std::size_t row, std::size_t col)
-{
-    if constexpr (Regular)
-        return data[row * Cols + col];
-    else
-        return data[col * Cols + row]; // transpose
+    std::shared_ptr<mm::basic_matrix<T, Rows, Cols>> M;
+
+    matrix<T, Rows, Cols> shallow_cpy()
+    {
+        return matrix<T, Rows, Cols>(M, transposed);
+    }
+
+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, bool Regular>
-const T& mm::access<T, Rows, Cols, Regular>::at(std::size_t row, std::size_t col) const
-{
-    if constexpr (Regular)
-        return data[row * Cols + col];
-    else
-        return data[col * Cols + row]; // transpose
+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, bool Regular>
-auto mm::access<T, Rows, Cols, Regular>::operator[](std::size_t index)
-{
-    if constexpr (this->is_row_vec() || this->is_col_vec())
-        return data.at(index);
-    else if (Regular)
-        return mm::row_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
-    else
-        return mm::col_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+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, bool Regular>
-auto mm::access<T, Rows, Cols, Regular>::operator[](std::size_t index) const
-{
-    if constexpr (this->is_row_vec() || this->is_col_vec())
-        return data.at(index);
-    else if (Regular)
-        return mm::const_row_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
-    else
-        return mm::const_col_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+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;
 }
 
-/* Square interface implementation */
+// 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
+) {
+    static_assert(P1 == P2, "invalid matrix multiplication");
+    assert(a.cols() == b.rows());
 
-template<typename T, std::size_t N>
-T mm::square_interface<T, N>::trace() 
-{
-    T sum = 0;
-    for (mm::diag_iterator<T, N> it(*this, 0); it.ok(); ++it)
-        sum += *it;
+    mm::matrix<T, M, N> result;
+    mm::matrix<T, P2, N> bt = b.t(); // weak transposition
 
-    return sum;
+    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;
 }
 
+// transposed 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, P1, M>& a,
+    const mm::matrix<T, P2, N>& b
+) {
+    static_assert(P1 == P2, "invalid matrix multiplication");
+    assert(a.cols() == b.rows());
+
+    mm::matrix<T, M, N> result;
+    mm::matrix<T, P2, N> bt = b.t(); // weak transposition
+
+    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;
+}*/
+
+
+/*template<typename T, std::size_t N>
+void mm::square_matrix<T, N>::transpose() {
+    for (unsigned row = 0; row < N; row++)
+        for (unsigned col = 0; col < row; col++)
+            std::swap(this->at(row, col), this->at(col, row));
+}*/
+
+
+/*
+ * 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) << ", ";
+        }
+        os << std::setw(NumW) << m.at(index, m.cols()-1) << " ]\n";
+    }
+
+    return os;
+}
diff --git a/test/matrix_example.cpp b/test/matrix_example.cpp
index e7f6580..291feab 100644
--- a/test/matrix_example.cpp
+++ b/test/matrix_example.cpp
@@ -12,14 +12,15 @@ int main(int argc, char *argv[]) {
     std::cout << "a = \n" << a;
     std::cout << "b = \n" << b;
     std::cout << "c = \n" << c;
-    std::cout << std::endl;
-
+    std::cout << std::endl; 
+    
     // access elements
     std::cout << "Access elements" << std::endl;
     std::cout << "a.at(2,0) = " << a.at(2, 0) << std::endl;
     std::cout << "a[2][0]   = " << a[2][0] << std::endl;;
     std::cout << std::endl;
     
+    /*
     // basic operations
     std::cout << "Basic operations" << std::endl;
     std::cout << "a + b = \n" << a + b;
@@ -28,26 +29,25 @@ int main(int argc, char *argv[]) {
     std::cout << "a * 2 = \n" << a * 2;
     std::cout << "2 * a = \n" << 2 * a;
     std::cout << "a.td() = \n" << a.td(); // or a.trasposed();
-    std::cout << std::endl;
+    std::cout << std::endl;*/
 
     // special matrices
-    mm::square_matrix<std::complex<int>, 2> f {{{2, 3}, {1, 4}}, {{6, 1}, {-3, 4}}};
+    /*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;
+    std::cout << "tr(f) = " << f.tr();  //or f.trace()  << std::endl;
 
-    f.t();
-    std::cout << "after in place transpose f.t(), f = \n" << f; 
+    mm::t_square_matrix<std::complex<int>, 2>& 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;
+    std::cout << std::endl; */
 
 
 	return 0;