diff --git a/CMakeLists.txt b/CMakeLists.txt
index d88bafa..a2a5ad2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -29,7 +29,7 @@ endif()
 ## Sources & Public Headers lists
 ##======================================================================================================================
 set ( SOURCES
-      src/matrix_impl64.cpp
+      src/matrix.cpp
     )
 
 ##======================================================================================================================
diff --git a/include/rotgen/detail/generators.hpp b/include/rotgen/detail/generators.hpp
new file mode 100644
index 0000000..049a249
--- /dev/null
+++ b/include/rotgen/detail/generators.hpp
@@ -0,0 +1,12 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#pragma once
+
+#define ROTGEN_CAT(a, ...)             ROTGEN_PRIMITIVE_CAT(a, __VA_ARGS__)
+#define ROTGEN_PRIMITIVE_CAT(a, ...)   a ## __VA_ARGS__
+#define ROTGEN_MATRIX_NAME(BASE,ID,SUFFIX) ROTGEN_CAT(BASE,ROTGEN_CAT(ID,SUFFIX))
diff --git a/include/rotgen/detail/static_info.hpp b/include/rotgen/detail/static_info.hpp
new file mode 100644
index 0000000..39fc391
--- /dev/null
+++ b/include/rotgen/detail/static_info.hpp
@@ -0,0 +1,27 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#pragma once
+
+namespace rotgen
+{
+  struct matrix_static_info
+  {
+    int rows      = -1;
+    int cols      = -1;
+    int options   =  0;
+    int max_rows  = -1;
+    int max_cols  = -1;
+  };
+
+  inline constexpr int Dynamic    = -1;
+  inline constexpr int Infinity   = -1;
+  inline constexpr int AutoAlign  =  0;
+  inline constexpr int ColMajor   =  0;
+  inline constexpr int RowMajor   =  1;
+  inline constexpr int DotnAlgin  =  2;
+}
\ No newline at end of file
diff --git a/include/rotgen/impl/matrix.hpp b/include/rotgen/impl/matrix.hpp
new file mode 100644
index 0000000..98379b2
--- /dev/null
+++ b/include/rotgen/impl/matrix.hpp
@@ -0,0 +1,55 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#pragma once
+
+#include <rotgen/detail/generators.hpp>
+#include <rotgen/detail/static_info.hpp>
+#include <initializer_list>
+#include <cstddef>
+#include <memory>
+
+namespace rotgen
+{
+  #define SIZE 64
+  #define TYPE double
+
+  #define MATRIX ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col)
+  #include <rotgen/impl/matrix_model.hpp>
+  #undef MATRIX
+
+  #define MATRIX ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row)
+  #include <rotgen/impl/matrix_model.hpp>
+  #undef MATRIX
+
+  #undef SIZE
+  #undef TYPE
+
+  #define SIZE 32
+  #define TYPE float
+
+  #define MATRIX ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col)
+  #include <rotgen/impl/matrix_model.hpp>
+  #undef MATRIX
+
+  #define MATRIX ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row)
+  #include <rotgen/impl/matrix_model.hpp>
+  #undef MATRIX
+
+  #undef SIZE
+  #undef TYPE
+
+  template<typename Scalar,int Options> struct find_matrix_impl;
+
+  template<> struct find_matrix_impl<float , ColMajor> { using type = matrix_impl32_col; };
+  template<> struct find_matrix_impl<float , RowMajor> { using type = matrix_impl32_row; };
+  template<> struct find_matrix_impl<double, ColMajor> { using type = matrix_impl64_col; };
+  template<> struct find_matrix_impl<double, RowMajor> { using type = matrix_impl64_row; };
+
+  template<typename Scalar,int Options>
+  using find_matrix = typename find_matrix_impl<Scalar,(Options & 1)>::type;
+}
diff --git a/include/rotgen/impl/matrix_impl64.hpp b/include/rotgen/impl/matrix_impl64.hpp
deleted file mode 100644
index d731805..0000000
--- a/include/rotgen/impl/matrix_impl64.hpp
+++ /dev/null
@@ -1,86 +0,0 @@
-//==================================================================================================
-/*
-  ROTGEN - Runtime Overlay for Eigen
-  Copyright : CODE RECKONS
-  SPDX-License-Identifier: BSL-1.0
-*/
-//==================================================================================================
-#pragma once
-
-#include <memory>
-#include <cstddef>
-#include <initializer_list>
-
-namespace rotgen
-{
-  class matrix_impl64
-  {
-    public:
-    matrix_impl64(std::size_t rows  = 0, std::size_t cols = 0);
-    matrix_impl64(std::size_t rows, std::size_t cols,std::initializer_list<double> init);
-
-    matrix_impl64(std::initializer_list<std::initializer_list<double>> init);
-
-    matrix_impl64(matrix_impl64 const& other);
-    matrix_impl64(matrix_impl64&&) noexcept;
-
-    matrix_impl64& operator=(matrix_impl64 const& other);
-    matrix_impl64& operator=(matrix_impl64&&) noexcept;
-
-    ~matrix_impl64();
-
-    std::size_t rows() const;
-    std::size_t cols() const;
-
-    std::size_t size() const;
-    void resize(std::size_t new_rows, std::size_t new_cols);
-    void conservativeResize(std::size_t new_rows, std::size_t new_cols);
-
-    matrix_impl64 transpose() const;
-    matrix_impl64 conjugate() const;
-    matrix_impl64 adjoint() const;
-
-    void transposeInPlace();
-    void adjointInPlace();
-
-    double sum() const;
-    double prod() const;
-    double mean() const;
-    double maxCoeff() const;
-    double maxCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const;
-    double minCoeff() const;
-    double minCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const;
-    double trace() const;
-
-    double squaredNorm() const;
-    double norm() const;
-    double lpNorm(int p) const;
-
-    double&       operator()(std::size_t i, std::size_t j);
-    double const& operator()(std::size_t i, std::size_t j) const;
-
-    double&       operator()(std::size_t index);
-    double const& operator()(std::size_t index) const;
-
-    matrix_impl64& operator+=(matrix_impl64 const& rhs);
-    matrix_impl64& operator-=(matrix_impl64 const& rhs);
-    matrix_impl64 operator-() const;
-    matrix_impl64& operator*=(matrix_impl64 const& rhs);
-    matrix_impl64& operator*=(double  d);
-    matrix_impl64& operator/=(double  d);
-
-    friend std::ostream& operator<<(std::ostream&,matrix_impl64 const&);
-    friend bool operator==(matrix_impl64 const& lhs, matrix_impl64 const& rhs);
-
-    const double* data() const;
-
-    static matrix_impl64 Zero(std::size_t rows, std::size_t cols);
-    static matrix_impl64 Constant(std::size_t rows, std::size_t cols, double value);
-    static matrix_impl64 Random(std::size_t rows, std::size_t cols);
-    static matrix_impl64 Identity(std::size_t rows, std::size_t cols);
-
-    private:
-    struct payload;
-    std::unique_ptr<payload> storage_;
-  };
-}
diff --git a/include/rotgen/impl/matrix_model.hpp b/include/rotgen/impl/matrix_model.hpp
new file mode 100644
index 0000000..43441ce
--- /dev/null
+++ b/include/rotgen/impl/matrix_model.hpp
@@ -0,0 +1,83 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+
+//==================================================================================================
+/*
+  This file is a X-File to generate various matrix_impl_* definitions variant
+*/
+//==================================================================================================
+class MATRIX
+{
+  public:
+  MATRIX(std::size_t rows  = 0, std::size_t cols = 0);
+  MATRIX(std::size_t rows, std::size_t cols,std::initializer_list<TYPE> init);
+
+  MATRIX(std::initializer_list<std::initializer_list<TYPE>> init);
+
+  MATRIX(MATRIX const& other);
+  MATRIX(MATRIX&&) noexcept;
+
+  MATRIX& operator=(MATRIX const& other);
+  MATRIX& operator=(MATRIX&&) noexcept;
+
+  ~MATRIX();
+
+  std::size_t rows() const;
+  std::size_t cols() const;
+
+  std::size_t size() const;
+  void resize(std::size_t new_rows, std::size_t new_cols);
+  void conservativeResize(std::size_t new_rows, std::size_t new_cols);
+
+  MATRIX transpose() const;
+  MATRIX conjugate() const;
+  MATRIX adjoint() const;
+
+  void transposeInPlace();
+  void adjointInPlace();
+
+  TYPE sum()      const;
+  TYPE prod()     const;
+  TYPE mean()     const;
+  TYPE trace()    const;
+  TYPE maxCoeff() const;
+  TYPE minCoeff() const;
+  TYPE maxCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const;
+  TYPE minCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const;
+
+  TYPE squaredNorm()  const;
+  TYPE norm()         const;
+  TYPE lpNorm(int p)  const;
+
+  TYPE&       operator()(std::size_t i, std::size_t j);
+  TYPE const& operator()(std::size_t i, std::size_t j) const;
+
+  TYPE&       operator()(std::size_t index);
+  TYPE const& operator()(std::size_t index) const;
+
+  MATRIX& operator+=(MATRIX const& rhs);
+  MATRIX& operator-=(MATRIX const& rhs);
+  MATRIX operator-() const;
+  MATRIX& operator*=(MATRIX const& rhs);
+  MATRIX& operator*=(TYPE  d);
+  MATRIX& operator/=(TYPE  d);
+
+  friend std::ostream& operator<<(std::ostream&,MATRIX const&);
+  friend bool operator==(MATRIX const& lhs, MATRIX const& rhs);
+
+  const TYPE* data() const;
+
+  static MATRIX Zero(std::size_t rows, std::size_t cols);
+  static MATRIX Constant(std::size_t rows, std::size_t cols, TYPE value);
+  static MATRIX Random(std::size_t rows, std::size_t cols);
+  static MATRIX Identity(std::size_t rows, std::size_t cols);
+
+  private:
+  struct payload;
+  std::unique_ptr<payload> storage_;
+};
\ No newline at end of file
diff --git a/include/rotgen/matrix.hpp b/include/rotgen/matrix.hpp
index 58638b9..efa35e8 100644
--- a/include/rotgen/matrix.hpp
+++ b/include/rotgen/matrix.hpp
@@ -7,28 +7,30 @@
 //==================================================================================================
 #pragma once
 
-#include <rotgen/impl/matrix_impl64.hpp>
+#include <rotgen/impl/matrix.hpp>
 #include <initializer_list>
 #include <cassert>
 
 namespace rotgen
 {
-  inline constexpr int Infinity = -1;
-
-  template< typename Scalar, int Rows    = -1  , int Cols    = -1
-          , int Options = 0, int MaxRows = Rows, int MaxCols = Cols
+  template< typename Scalar
+          , int Rows    = Dynamic , int Cols    = Dynamic
+          , int Options = ColMajor
+          , int MaxRows = Rows    , int MaxCols = Cols
           >
-  class matrix : public matrix_impl64
+  class matrix : public find_matrix<Scalar,Options>
   {
-    using parent = matrix_impl64;
+    using parent = find_matrix<Scalar,Options>;
 
     public:
+    using value_type = Scalar;
+
     matrix()                              : parent(Rows==-1?0:Rows,Cols==-1?0:Cols) {}
 
     matrix(int r, int c) : parent(r, c)
     {
-       if constexpr(Rows != -1) assert(r == Rows &&  "Mismatched between dynamic and static row size");
-       if constexpr(Cols != -1) assert(c == Cols &&  "Mismatched between dynamic and static column size");
+      if constexpr(Rows != -1) assert(r == Rows &&  "Mismatched between dynamic and static row size");
+      if constexpr(Cols != -1) assert(c == Cols &&  "Mismatched between dynamic and static column size");
     }
 
     matrix(parent const& base) : parent(base) {}
@@ -50,14 +52,12 @@ namespace rotgen
           : parent(Rows, Cols, {s0,static_cast<Scalar>(init)...})
     {}
 
-    void resize(int new_rows, int new_cols)
-    requires(Rows == -1 && Cols == -1)
+    void resize(int new_rows, int new_cols) requires(Rows == -1 && Cols == -1)
     {
       parent::resize(new_rows, new_cols);
     }
 
-    void conservativeResize(int new_rows, int new_cols)
-    requires(Rows == -1 && Cols == -1)
+    void conservativeResize(int new_rows, int new_cols) requires(Rows == -1 && Cols == -1)
     {
       parent::conservativeResize(new_rows, new_cols);
     }
@@ -78,8 +78,7 @@ namespace rotgen
     }
 
     void transposeInPlace() { parent::transposeInPlace(); }
-
-    void adjointInPlace() { parent::adjointInPlace(); }
+    void adjointInPlace()   { parent::adjointInPlace();   }
 
     friend bool operator==(matrix const& lhs, matrix const& rhs)
     {
@@ -109,20 +108,19 @@ namespace rotgen
       return *this;
     }
 
-    matrix& operator*=(double rhs)
+    matrix& operator*=(Scalar rhs)
     {
       static_cast<parent&>(*this) *= rhs;
       return *this;
     }
 
-    matrix& operator/=(double rhs)
+    matrix& operator/=(Scalar rhs)
     {
       static_cast<parent&>(*this) /= rhs;
       return *this;
     }
 
-    static matrix Zero()
-    requires (Rows != -1 && Cols != -1)
+    static matrix Zero() requires (Rows != -1 && Cols != -1)
     {
       return parent::Zero(Rows, Cols);
     }
@@ -134,21 +132,19 @@ namespace rotgen
       return parent::Zero(rows, cols);
     }
 
-    static matrix Constant(Scalar value)
-    requires (Rows != -1 && Cols != -1)
+    static matrix Constant(Scalar value) requires (Rows != -1 && Cols != -1)
     {
-      return parent::Constant(Rows, Cols, static_cast<double>(value));
+      return parent::Constant(Rows, Cols, static_cast<Scalar>(value));
     }
 
     static matrix Constant(int rows, int cols, Scalar value)
     {
       if constexpr(Rows != -1) assert(rows == Rows &&  "Mismatched between dynamic and static row size");
       if constexpr(Cols != -1) assert(cols == Cols &&  "Mismatched between dynamic and static column size");
-      return parent::Constant(rows, cols, static_cast<double>(value));
+      return parent::Constant(rows, cols, static_cast<Scalar>(value));
     }
 
-    static matrix Random()
-    requires (Rows != -1 && Cols != -1)
+    static matrix Random() requires (Rows != -1 && Cols != -1)
     {
       return parent::Random(Rows, Cols);
     }
@@ -160,8 +156,7 @@ namespace rotgen
       return parent::Random(rows, cols);
     }
 
-    static matrix Identity()
-    requires (Rows != -1 && Cols != -1)
+    static matrix Identity() requires (Rows != -1 && Cols != -1)
     {
       return parent::Identity(Rows, Cols);
     }
@@ -174,7 +169,8 @@ namespace rotgen
     }
 
     template<int P>
-    double lpNorm() const {
+    Scalar lpNorm() const
+    {
       assert(P == 1 || P == 2 || P == Infinity);
       return parent::lpNorm(P);
     }
@@ -203,20 +199,20 @@ namespace rotgen
   }
 
   template<typename S, int R, int C, int O, int MR, int MC>
-  matrix<S,R,C,O,MR,MC> operator*(matrix<S,R,C,O,MR,MC> const& lhs, double rhs)
+  matrix<S,R,C,O,MR,MC> operator*(matrix<S,R,C,O,MR,MC> const& lhs, S rhs)
   {
     matrix<S,R,C,O,MR,MC> that(lhs);
     return that *= rhs;
   }
 
   template<typename S, int R, int C, int O, int MR, int MC>
-  matrix<S,R,C,O,MR,MC> operator*(double lhs, matrix<S,R,C,O,MR,MC> const&  rhs)
+  matrix<S,R,C,O,MR,MC> operator*(S lhs, matrix<S,R,C,O,MR,MC> const&  rhs)
   {
     return rhs * lhs;
   }
 
   template<typename S, int R, int C, int O, int MR, int MC>
-  matrix<S,R,C,O,MR,MC> operator/(matrix<S,R,C,O,MR,MC> const& lhs, double rhs)
+  matrix<S,R,C,O,MR,MC> operator/(matrix<S,R,C,O,MR,MC> const& lhs, S rhs)
   {
     matrix<S,R,C,O,MR,MC> that(lhs);
     return that /= rhs;
diff --git a/src/matrix.cpp b/src/matrix.cpp
new file mode 100644
index 0000000..66ca20d
--- /dev/null
+++ b/src/matrix.cpp
@@ -0,0 +1,49 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#include <rotgen/detail/generators.hpp>
+#include <rotgen/impl/matrix.hpp>
+#include <Eigen/Dense>
+
+namespace rotgen
+{
+  #define SIZE 64
+  #define TYPE double
+  #define STORAGE_ORDER Eigen::ColMajor
+
+  #define MATRIX ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col)
+  #include "matrix_model.cpp"
+  #undef MATRIX
+  #undef STORAGE_ORDER
+
+  #define STORAGE_ORDER Eigen::RowMajor
+  #define MATRIX ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row)
+  #include "matrix_model.cpp"
+  #undef MATRIX
+  #undef STORAGE_ORDER
+
+  #undef SIZE
+  #undef TYPE
+
+  #define SIZE 32
+  #define TYPE float
+  #define STORAGE_ORDER Eigen::ColMajor
+
+  #define MATRIX ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_col)
+  #include "matrix_model.cpp"
+  #undef MATRIX
+  #undef STORAGE_ORDER
+
+  #define STORAGE_ORDER Eigen::RowMajor
+  #define MATRIX ROTGEN_MATRIX_NAME(matrix_impl,SIZE,_row)
+  #include "matrix_model.cpp"
+  #undef MATRIX
+  #undef STORAGE_ORDER
+
+  #undef SIZE
+  #undef TYPE
+}
diff --git a/src/matrix_impl64.cpp b/src/matrix_impl64.cpp
deleted file mode 100644
index 6568e76..0000000
--- a/src/matrix_impl64.cpp
+++ /dev/null
@@ -1,219 +0,0 @@
-//==================================================================================================
-/*
-  ROTGEN - Runtime Overlay for Eigen
-  Copyright : CODE RECKONS
-  SPDX-License-Identifier: BSL-1.0
-*/
-//==================================================================================================
-#include <rotgen/matrix.hpp>
-#include <Eigen/Dense>
-
-namespace rotgen
-{
-  //================================================================================================
-  // Internal paylod
-  //================================================================================================
-  struct matrix_impl64::payload
-  {
-    using data_type = Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>;
-
-    data_type data;
-    payload(std::size_t r=0, std::size_t c=0) : data(r, c) {}
-    payload(std::initializer_list<std::initializer_list<double>> init) : data(init) {}
-    payload(data_type&& matrix) : data(std::move(matrix)) {}
-  };
-
-  //==================================================================================================
-  // Constructors & Special Members
-  //==================================================================================================
-  matrix_impl64::matrix_impl64(std::size_t r, std::size_t c) : storage_(std::make_unique<payload>(r,c)) {}
-
-  matrix_impl64::matrix_impl64(std::initializer_list<std::initializer_list<double>> init)
-                : storage_(std::make_unique<payload>(init))
-  {}
-
-  matrix_impl64::matrix_impl64(std::size_t r, std::size_t c,std::initializer_list<double> init)
-  : matrix_impl64(r,c)
-  {
-    auto first = init.begin();
-    for(std::size_t i=0; i < init.size(); i++)
-      (*this)(i) = first[i];
-  }
-
-  matrix_impl64::matrix_impl64(matrix_impl64 const& o)
-                : matrix_impl64(o.rows(),o.cols())
-  {
-    storage_->data = o.storage_->data;
-  }
-
-  matrix_impl64::matrix_impl64(matrix_impl64&&) noexcept            = default;
-
-
-  matrix_impl64&  matrix_impl64::operator=(matrix_impl64 const& o)
-  {
-    if (this != &o) storage_->data = o.storage_->data;
-    return *this;
-  }
-
-
-  matrix_impl64& matrix_impl64::operator=(matrix_impl64&&) noexcept = default;
-
-  matrix_impl64::~matrix_impl64()                                   = default;
-
-  //==================================================================================================
-  // Matrix API
-  //==================================================================================================
-  std::size_t matrix_impl64::rows() const { return static_cast<std::size_t>(storage_->data.rows()); }
-  std::size_t matrix_impl64::cols() const { return static_cast<std::size_t>(storage_->data.cols()); }
-
-  std::size_t matrix_impl64::size() const { return static_cast<std::size_t>(storage_->data.size()); }
-  void matrix_impl64::resize(std::size_t new_rows, std::size_t new_cols) { storage_->data.resize(new_rows, new_cols); }
-  void matrix_impl64::conservativeResize(std::size_t new_rows, std::size_t new_cols) { storage_->data.conservativeResize(new_rows, new_cols); }
-
-  double&       matrix_impl64::operator()(std::size_t i, std::size_t j)       { return storage_->data(i,j); }
-  double const& matrix_impl64::operator()(std::size_t i, std::size_t j) const { return storage_->data(i,j); }
-
-  double&       matrix_impl64::operator()(std::size_t index)       { return storage_->data(index); }
-  double const& matrix_impl64::operator()(std::size_t index) const { return storage_->data(index); }
-
-  const double* matrix_impl64::data() const { return storage_->data.data(); }
-
-  matrix_impl64 matrix_impl64::transpose() const
-  {
-    matrix_impl64 result(*this);
-    result.storage_->data.transposeInPlace();
-    return result;
-
-  }
-
-  matrix_impl64 matrix_impl64::conjugate() const
-  {
-    matrix_impl64 result(*this);
-    result.storage_->data = storage_->data.conjugate();
-    return result;
-  }
-
-  matrix_impl64 matrix_impl64::adjoint() const
-  {
-    matrix_impl64 result(*this);
-    result.storage_->data.adjointInPlace();
-    return result;
-  }
-
-  void matrix_impl64::transposeInPlace()
-  {
-    storage_->data.transposeInPlace();
-  }
-
-  void matrix_impl64::adjointInPlace()
-  {
-    storage_->data.adjointInPlace();
-  }
-
-  double matrix_impl64::sum() const { return storage_->data.sum(); }
-  double matrix_impl64::prod() const { return storage_->data.prod(); }
-  double matrix_impl64::mean() const { return storage_->data.mean(); }
-
-  double matrix_impl64::minCoeff() const { return storage_->data.minCoeff(); }
-  double matrix_impl64::minCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const { return storage_->data.minCoeff(row, col); }
-
-  double matrix_impl64::maxCoeff() const { return storage_->data.maxCoeff(); }
-  double matrix_impl64::maxCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const { return storage_->data.maxCoeff(row, col); }
-
-  double matrix_impl64::trace() const { return storage_->data.trace(); }
-
-  double matrix_impl64::squaredNorm() const { return storage_->data.squaredNorm(); }
-  double matrix_impl64::norm() const { return storage_->data.norm(); }
-
-  double matrix_impl64::lpNorm(int p) const {
-    if (p == 1)
-      return storage_->data.lpNorm<1>();
-    else if (p == 2)
-      return storage_->data.lpNorm<2>();
-    else
-      return storage_->data.lpNorm<Eigen::Infinity>();
-  }
-
-  //==================================================================================================
-  // Operators
-  //==================================================================================================
-  std::ostream& operator<<(std::ostream& os,matrix_impl64 const& m)
-  {
-    return os << m.storage_->data;
-  }
-
-  bool operator==(matrix_impl64 const& lhs, matrix_impl64 const& rhs)
-  {
-    return lhs.storage_->data == rhs.storage_->data;
-  }
-
-  matrix_impl64& matrix_impl64::operator+=(matrix_impl64 const& rhs)
-  {
-    storage_->data += rhs.storage_->data;
-    return *this;
-  }
-
-  matrix_impl64& matrix_impl64::operator-=(matrix_impl64 const& rhs)
-  {
-    storage_->data -= rhs.storage_->data;
-    return *this;
-  }
-
-  matrix_impl64 matrix_impl64::operator-() const
-  {
-    matrix_impl64 result(*this);
-    result.storage_->data = -result.storage_->data;
-    return result;
-  }
-
-  matrix_impl64& matrix_impl64::operator*=(matrix_impl64 const& rhs)
-  {
-    storage_->data *= rhs.storage_->data;
-    return *this;
-  }
-
-  matrix_impl64& matrix_impl64::operator*=(double s)
-  {
-    storage_->data *= s;
-    return *this;
-  }
-
-  matrix_impl64& matrix_impl64::operator/=(double s)
-  {
-    storage_->data /= s;
-    return *this;
-  }
-
-  //==================================================================================================
-  //==================================================================================================
-  // Static functions
-  //==================================================================================================
-
-  matrix_impl64 matrix_impl64::Zero(std::size_t rows, std::size_t cols) {
-    matrix_impl64 m;
-    m.storage_ = std::make_unique<payload>(payload::data_type::Zero(rows, cols));
-    return m;
-  }
-
-  matrix_impl64 matrix_impl64::Constant(std::size_t rows, std::size_t cols, double value)
-  {
-    matrix_impl64 m;
-    m.storage_ = std::make_unique<payload>(payload::data_type::Constant(rows, cols, value));
-    return m;
-  }
-
-  matrix_impl64 matrix_impl64::Random(std::size_t rows, std::size_t cols)
-  {
-    matrix_impl64 m;
-    m.storage_ = std::make_unique<payload>(payload::data_type::Random(rows, cols));
-    return m;
-  }
-
-  matrix_impl64 matrix_impl64::Identity(std::size_t rows, std::size_t cols)
-  {
-    matrix_impl64 m;
-    m.storage_ = std::make_unique<payload>(payload::data_type::Identity(rows, cols));
-    return m;
-  }
-
-}
\ No newline at end of file
diff --git a/src/matrix_model.cpp b/src/matrix_model.cpp
new file mode 100644
index 0000000..c8dbac3
--- /dev/null
+++ b/src/matrix_model.cpp
@@ -0,0 +1,222 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+
+//==================================================================================================
+/*
+  This file is a X-File to generate various matrix_impl_* definitions variant
+*/
+//==================================================================================================
+
+//==================================================================================================
+// Internal payload
+//==================================================================================================
+struct MATRIX::payload
+{
+  using data_type = Eigen::Matrix<TYPE,Eigen::Dynamic,Eigen::Dynamic,STORAGE_ORDER>;
+
+  data_type data;
+  payload(std::size_t r=0, std::size_t c=0) : data(r, c) {}
+  payload(std::initializer_list<std::initializer_list<TYPE>> init) : data(init) {}
+  payload(data_type&& matrix) : data(std::move(matrix)) {}
+};
+
+//==================================================================================================
+// Constructors & Special Members
+//==================================================================================================
+MATRIX::MATRIX(std::size_t r, std::size_t c) : storage_(std::make_unique<payload>(r,c)) {}
+
+MATRIX::MATRIX(std::initializer_list<std::initializer_list<TYPE>> init)
+              : storage_(std::make_unique<payload>(init))
+{}
+
+MATRIX::MATRIX(std::size_t r, std::size_t c,std::initializer_list<TYPE> init)
+: MATRIX(r,c)
+{
+  auto first = init.begin();
+  for(std::size_t i=0; i < init.size(); i++)
+    (*this)(i) = first[i];
+}
+
+MATRIX::MATRIX(MATRIX const& o)
+              : MATRIX(o.rows(),o.cols())
+{
+  storage_->data = o.storage_->data;
+}
+
+MATRIX::MATRIX(MATRIX&&) noexcept            = default;
+
+
+MATRIX&  MATRIX::operator=(MATRIX const& o)
+{
+  if (this != &o) storage_->data = o.storage_->data;
+  return *this;
+}
+
+MATRIX& MATRIX::operator=(MATRIX&&) noexcept = default;
+
+MATRIX::~MATRIX()                                   = default;
+
+//==================================================================================================
+// Matrix API
+//==================================================================================================
+std::size_t MATRIX::rows() const { return static_cast<std::size_t>(storage_->data.rows()); }
+std::size_t MATRIX::cols() const { return static_cast<std::size_t>(storage_->data.cols()); }
+std::size_t MATRIX::size() const { return static_cast<std::size_t>(storage_->data.size()); }
+
+void MATRIX::resize(std::size_t new_rows, std::size_t new_cols)
+{
+  storage_->data.resize(new_rows, new_cols);
+}
+
+void MATRIX::conservativeResize(std::size_t new_rows, std::size_t new_cols)
+{
+  storage_->data.conservativeResize(new_rows, new_cols);
+}
+
+TYPE&       MATRIX::operator()(std::size_t i, std::size_t j)       { return storage_->data(i,j); }
+TYPE const& MATRIX::operator()(std::size_t i, std::size_t j) const { return storage_->data(i,j); }
+
+TYPE&       MATRIX::operator()(std::size_t index)       { return storage_->data(index); }
+TYPE const& MATRIX::operator()(std::size_t index) const { return storage_->data(index); }
+
+const TYPE* MATRIX::data() const { return storage_->data.data(); }
+
+MATRIX MATRIX::transpose() const
+{
+  MATRIX result(*this);
+  result.storage_->data.transposeInPlace();
+  return result;
+
+}
+
+MATRIX MATRIX::conjugate() const
+{
+  MATRIX result(*this);
+  result.storage_->data = storage_->data.conjugate();
+  return result;
+}
+
+MATRIX MATRIX::adjoint() const
+{
+  MATRIX result(*this);
+  result.storage_->data.adjointInPlace();
+  return result;
+}
+
+void MATRIX::transposeInPlace()
+{
+  storage_->data.transposeInPlace();
+}
+
+void MATRIX::adjointInPlace()
+{
+  storage_->data.adjointInPlace();
+}
+
+TYPE MATRIX::sum()    const { return storage_->data.sum(); }
+TYPE MATRIX::prod()   const { return storage_->data.prod(); }
+TYPE MATRIX::mean()   const { return storage_->data.mean(); }
+TYPE MATRIX::trace()  const { return storage_->data.trace(); }
+
+TYPE MATRIX::minCoeff() const { return storage_->data.minCoeff(); }
+TYPE MATRIX::maxCoeff() const { return storage_->data.maxCoeff(); }
+
+TYPE MATRIX::minCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const { return storage_->data.minCoeff(row, col); }
+TYPE MATRIX::maxCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const { return storage_->data.maxCoeff(row, col); }
+
+TYPE MATRIX::squaredNorm()  const { return storage_->data.squaredNorm(); }
+TYPE MATRIX::norm()         const { return storage_->data.norm(); }
+
+TYPE MATRIX::lpNorm(int p) const
+{
+  if (p == 1)       return storage_->data.lpNorm<1>();
+  else if (p == 2)  return storage_->data.lpNorm<2>();
+  else              return storage_->data.lpNorm<Eigen::Infinity>();
+}
+
+//==================================================================================================
+// Operators
+//==================================================================================================
+std::ostream& operator<<(std::ostream& os,MATRIX const& m)
+{
+  return os << m.storage_->data;
+}
+
+bool operator==(MATRIX const& lhs, MATRIX const& rhs)
+{
+  return lhs.storage_->data == rhs.storage_->data;
+}
+
+MATRIX& MATRIX::operator+=(MATRIX const& rhs)
+{
+  storage_->data += rhs.storage_->data;
+  return *this;
+}
+
+MATRIX& MATRIX::operator-=(MATRIX const& rhs)
+{
+  storage_->data -= rhs.storage_->data;
+  return *this;
+}
+
+MATRIX MATRIX::operator-() const
+{
+  MATRIX result(*this);
+  result.storage_->data = -result.storage_->data;
+  return result;
+}
+
+MATRIX& MATRIX::operator*=(MATRIX const& rhs)
+{
+  storage_->data *= rhs.storage_->data;
+  return *this;
+}
+
+MATRIX& MATRIX::operator*=(TYPE s)
+{
+  storage_->data *= s;
+  return *this;
+}
+
+MATRIX& MATRIX::operator/=(TYPE s)
+{
+  storage_->data /= s;
+  return *this;
+}
+
+//==================================================================================================
+// Static functions
+//==================================================================================================
+
+MATRIX MATRIX::Zero(std::size_t rows, std::size_t cols)
+{
+  MATRIX m;
+  m.storage_ = std::make_unique<payload>(payload::data_type::Zero(rows, cols));
+  return m;
+}
+
+MATRIX MATRIX::Constant(std::size_t rows, std::size_t cols, TYPE value)
+{
+  MATRIX m;
+  m.storage_ = std::make_unique<payload>(payload::data_type::Constant(rows, cols, value));
+  return m;
+}
+
+MATRIX MATRIX::Random(std::size_t rows, std::size_t cols)
+{
+  MATRIX m;
+  m.storage_ = std::make_unique<payload>(payload::data_type::Random(rows, cols));
+  return m;
+}
+
+MATRIX MATRIX::Identity(std::size_t rows, std::size_t cols)
+{
+  MATRIX m;
+  m.storage_ = std::make_unique<payload>(payload::data_type::Identity(rows, cols));
+  return m;
+}
\ No newline at end of file
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index b9dddab..87d5a48 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -17,4 +17,4 @@ target_link_libraries(rotgen_test INTERFACE rotgen Eigen3::Eigen)
 ##======================================================================================================================
 ## Unit Tests registration
 ##======================================================================================================================
-rotgen_glob_unit(PATTERN "basic/*.cpp"  INTERFACE rotgen_test)
\ No newline at end of file
+rotgen_glob_unit(PATTERN "unit/basic/*.cpp"  INTERFACE rotgen_test)
\ No newline at end of file
diff --git a/test/basic/arithmetic_functions.cpp b/test/basic/arithmetic_functions.cpp
deleted file mode 100644
index c61f505..0000000
--- a/test/basic/arithmetic_functions.cpp
+++ /dev/null
@@ -1,121 +0,0 @@
-//==================================================================================================
-/*
-  ROTGEN - Runtime Overlay for Eigen
-  Copyright : CODE RECKONS
-  SPDX-License-Identifier: BSL-1.0
-*/
-//==================================================================================================
-#define TTS_MAIN
-#include <rotgen/matrix.hpp>
-#include "tts.hpp"
-#include <Eigen/Dense>
-
-template <typename MatrixType>
-struct MatrixDescriptor
-{
-  std::size_t rows, cols;
-  std::function<void(MatrixType &, std::size_t, std::size_t)> init_fn;
-};
-
-template <typename MatrixType>
-void test_matrix_unary_ops(std::size_t rows, std::size_t cols,
-                           const std::function<void(MatrixType &, std::size_t, std::size_t)> &init_fn)
-{
-
-  MatrixType original(rows, cols);
-  MatrixType transposed_matrix(cols, rows);
-
-  for (std::size_t r = 0; r < rows; ++r)
-    for (std::size_t c = 0; c < cols; ++c)
-      init_fn(original, r, c);
-
-  for (std::size_t r = 0; r < rows; ++r)
-    for (std::size_t c = 0; c < cols; ++c)
-      transposed_matrix(c, r) = original(r, c);
-
-  TTS_EQUAL(original.transpose(), transposed_matrix);
-  TTS_EQUAL(original.conjugate(), original);
-  TTS_EQUAL(original.adjoint(), transposed_matrix);
-
-  original.transposeInPlace();
-  TTS_EQUAL(original, transposed_matrix);
-
-  original.transposeInPlace();
-  original.adjointInPlace();
-  TTS_EQUAL(original, transposed_matrix);
-}
-
-template <typename MatrixType>
-void test_matrix_scalar_reductions(std::size_t rows, std::size_t cols,
-                                   const std::function<void(MatrixType&, std::size_t, std::size_t)>& init_fn)
-{
-  using EigenMatrix = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
-
-  MatrixType original(rows, cols);
-  EigenMatrix ref(rows, cols);
-
-  for (std::size_t r = 0; r < rows; ++r)
-    for (std::size_t c = 0; c < cols; ++c)
-    {
-      init_fn(original, r, c);
-      ref(r, c) = original(r, c);
-    }
-
-  TTS_EQUAL(original.sum(), ref.sum());
-  TTS_EQUAL(original.prod(), ref.prod());
-  TTS_EQUAL(original.mean(), ref.mean());
-  TTS_EQUAL(original.maxCoeff(), ref.maxCoeff());
-  TTS_EQUAL(original.minCoeff(), ref.minCoeff());
-  TTS_EQUAL(original.trace(), ref.trace());
-
-  std::ptrdiff_t row, col, ref_row, ref_col;
-
-  double cmin = original.minCoeff(&row, &col);
-  double rmin = ref.minCoeff(&ref_row, &ref_col);
-
-  TTS_EQUAL(cmin, rmin);
-  TTS_EQUAL(row, ref_row);
-  TTS_EQUAL(col, ref_col);
-
-  double cmax = original.maxCoeff(&row, &col);
-  double rmax = ref.maxCoeff(&ref_row, &ref_col);
-
-  TTS_EQUAL(cmax, rmax);
-  TTS_EQUAL(row, ref_row);
-  TTS_EQUAL(col, ref_col);
-}
-
-TTS_CASE("Matrix unary operations: transpose, adjoint, conjugate")
-{
-  std::vector<MatrixDescriptor<rotgen::matrix<double>>> test_matrices = {
-      {3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r + 3 * c - 2.5; }},
-      {4, 9, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r*r + 3.12 * c + 6.87; }},
-      {2, 7, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 1.1 * (r - c); }},
-      {1, 5, [](auto &m, std::size_t r, std::size_t c) { m(r, c) = 9.99; }},
-      {4, 2, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 0.0; }},
-      {3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r == c) ? 1.0 : 0.0; }},
-      {2, 2, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r + c) * 1e-10; }},
-      {2, 2, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r + 1) * 1e+10; }},
-  };
-
-  for (auto const &desc : test_matrices)
-    test_matrix_unary_ops<rotgen::matrix<double>>(desc.rows, desc.cols, desc.init_fn);
-};
-
-TTS_CASE("Basic arithmetic reduction operations")
-{
-  std::vector<MatrixDescriptor<rotgen::matrix<double>>> test_matrices = {
-    {3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r + c; }},
-    {3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 0.0; }},
-    {2, 4, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = -r -c*c - 1234; }},
-    {4, 4, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 7.0; }},
-    {1, 1, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 42.0; }},
-    {4, 2, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = std::sin(r + c); }},
-    {1, 5, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = -1.5 * r + 2.56 * c; }},
-    {5, 7, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r == c ? 1.0 : 0.0); }},
-  };
-
-  for (const auto& [rows, cols, init_fn] : test_matrices)
-    test_matrix_scalar_reductions<rotgen::matrix<double>>(rows, cols, init_fn);
-};
-
diff --git a/test/basic/constructors.cpp b/test/basic/constructors.cpp
deleted file mode 100644
index ee01db8..0000000
--- a/test/basic/constructors.cpp
+++ /dev/null
@@ -1,175 +0,0 @@
-//==================================================================================================
-/*
-  ROTGEN - Runtime Overlay for Eigen
-  Copyright : CODE RECKONS
-  SPDX-License-Identifier: BSL-1.0
-*/
-//==================================================================================================
-#define TTS_MAIN
-#include <rotgen/matrix.hpp>
-#include "tts.hpp"
-
-TTS_CASE("Default matrix dynamic constructor")
-{
-  rotgen::matrix<double> matrix;
-
-  TTS_EQUAL(matrix.rows(), 0ULL);
-  TTS_EQUAL(matrix.cols(), 0ULL);
-};
-
-TTS_CASE("Default matrix static constructor")
-{
-  rotgen::matrix<double, 4, 9> matrix;
-
-  TTS_EQUAL(matrix.rows(), 4ULL);
-  TTS_EQUAL(matrix.cols(), 9ULL);
-};
-
-TTS_CASE("Dynamic matrix constructor with row and columns")
-{
-  rotgen::matrix<float> matrix(10, 5);
-
-  TTS_EQUAL(matrix.rows(), 10ULL);
-  TTS_EQUAL(matrix.cols(), 5ULL);
-};
-
-TTS_CASE("Static matrix constructor with row and columns")
-{
-  rotgen::matrix<double, 6, 11> matrix(6, 11);
-
-  TTS_EQUAL(matrix.rows(), 6ULL);
-  TTS_EQUAL(matrix.cols(), 11ULL);
-};
-
-TTS_CASE("Copy constructor produces identical but independent matrix")
-{
-  rotgen::matrix<double> a(3, 3);
-
-  for(std::size_t r=0;r<a.rows();r++)
-    for(std::size_t c=0;c<a.cols();c++)
-      a(r,c)   = r + c;
-
-  rotgen::matrix<double> b = a;
-
-  TTS_EQUAL(b.rows(), a.rows());
-  TTS_EQUAL(b.cols(), a.cols());
-
-  for(std::size_t r=0;r<a.rows();r++)
-    for(std::size_t c=0;c<a.cols();c++)
-      TTS_EQUAL(b(r,c), a(r,c));
-
-  TTS_EQUAL(b(0, 0), 0.0);
-  TTS_EQUAL(b(1, 0), 1.0);
-  TTS_EQUAL(b(2, 1), 3.0);
-
-  a(0,0) = 42.0;
-  TTS_NOT_EQUAL(b(0,0), a(0,0));
-};
-
-TTS_CASE("Copy constructor on default matrix")
-{
-  rotgen::matrix<double> a;
-  rotgen::matrix<double> b = a;
-  TTS_EQUAL(b.rows(), 0ULL);
-  TTS_EQUAL(b.cols(), 0ULL);
-};
-
-TTS_CASE("Copy constructor from const matrix")
-{
-  const rotgen::matrix<double> a(2, 2);
-  auto b = a;
-  TTS_EQUAL(b.rows(), 2ULL);
-  TTS_EQUAL(b.cols(), 2ULL);
-};
-
-TTS_CASE("Copy constructor on static matrix")
-{
-  rotgen::matrix<double, 2, 5> a;
-  rotgen::matrix<double, 2, 5> b = a;
-  TTS_EQUAL(b.rows(), 2ULL);
-  TTS_EQUAL(b.cols(), 5ULL);
-};
-
-/*
-TTS_CASE("Copy constructor on static/dynamic matrix")
-{
-  rotgen::matrix<double, 11, 4> a;
-  rotgen::matrix<double> b = a;
-  TTS_EQUAL(b.rows(), 11);
-  TTS_EQUAL(b.cols(), 4);
-};
-
-TTS_CASE("Copy constructor on dynamic/static matrix")
-{
-  rotgen::matrix<double> a(5, 7);
-  rotgen::matrix<double, 5, 7> b = a;
-  TTS_EQUAL(b.rows(), 5);
-  TTS_EQUAL(b.cols(), 7);
-};*/
-
-TTS_CASE("Move constructor transfers contents")
-{
-  rotgen::matrix<double> a(3, 3);
-  a(1,1) = 7;
-  auto ptr = a.data();
-
- rotgen::matrix<double> b = std::move(a);
-
-  TTS_EQUAL(b(1,1), 7);
-  TTS_EQUAL(b.rows(), 3ULL);
-  TTS_EQUAL(b.cols(), 3ULL);
-  TTS_EXPECT(b.data() == ptr);
-//  TTS_EXPECT(a.data() == nullptr);
-};
-
-TTS_CASE("Move constructor from Rvalue")
-{
-  rotgen::matrix<double> b = rotgen::matrix<double>(2, 2);
-  TTS_EQUAL(b.rows(), 2ULL);
-  TTS_EQUAL(b.cols(), 2ULL);
-};
-
-TTS_CASE("Constructor from Initializer list")
-{
-  rotgen::matrix<double,1,1> b1 = {3.5};
-  TTS_EQUAL(b1.rows(), 1ULL);
-  TTS_EQUAL(b1.cols(), 1ULL);
-  TTS_EQUAL(b1(0)    , 3.5);
-
-  rotgen::matrix<double,5,1> b9 = {0.01, 0.1, 1, 10, 100};
-  TTS_EQUAL(b9.rows(), 5ULL);
-  TTS_EQUAL(b9.cols(), 1ULL);
-
-  double i = 0.01;
-  for(std::size_t r=0;r<b9.rows();++r) {
-    TTS_EQUAL(b9(r,0), i);
-    i *= 10;
-  }
-
-  rotgen::matrix<double,1,3> b13 = {1.2,2.3,3.4};
-  TTS_EQUAL(b13.rows(), 1ULL);
-  TTS_EQUAL(b13.cols(), 3ULL);
-  TTS_EQUAL(b13(0)    , 1.2);
-  TTS_EQUAL(b13(1)    , 2.3);
-  TTS_EQUAL(b13(2)    , 3.4);
-};
-
-TTS_CASE("Constructor from Initializer list of rows")
-{
-  rotgen::matrix<double> b1 = {{3.5}};
-  TTS_EQUAL(b1.rows(), 1ULL);
-  TTS_EQUAL(b1.cols(), 1ULL);
-  TTS_EQUAL(b1(0,0)  , 3.5);
-
-  rotgen::matrix<double> b23  = { {1.2,2.3,3.4}
-    , {10,200,3000}
-  };
-  TTS_EQUAL(b23.rows(), 2ULL);
-  TTS_EQUAL(b23.cols(), 3ULL);
-  TTS_EQUAL(b23(0,0)    , 1.2);
-  TTS_EQUAL(b23(0,1)    , 2.3);
-  TTS_EQUAL(b23(0,2)    , 3.4);
-  TTS_EQUAL(b23(1,0)    , 10);
-  TTS_EQUAL(b23(1,1)    , 200);
-  TTS_EQUAL(b23(1,2)    , 3000);
-};
\ No newline at end of file
diff --git a/test/basic/norms.cpp b/test/basic/norms.cpp
deleted file mode 100644
index 09ce7cd..0000000
--- a/test/basic/norms.cpp
+++ /dev/null
@@ -1,61 +0,0 @@
-//==================================================================================================
-/*
-  ROTGEN - Runtime Overlay for Eigen
-  Copyright : CODE RECKONS
-  SPDX-License-Identifier: BSL-1.0
-*/
-//==================================================================================================
-#define TTS_MAIN
-#include <rotgen/matrix.hpp>
-#include "tts.hpp"
-#include <Eigen/Dense>
-
-template <typename MatrixType>
-struct MatrixDescriptor
-{
-  std::size_t rows, cols;
-  std::function<void(MatrixType &, std::size_t, std::size_t)> init_fn;
-};
-
-template <typename MatrixType>
-void test_matrix_norms(std::size_t rows, std::size_t cols,
-                       const std::function<void(MatrixType &, std::size_t, std::size_t)> &init_fn)
-{
-  MatrixType matrix(rows, cols);
-  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> ref(rows, cols);
-
-  for (std::size_t r = 0; r < rows; ++r)
-    for (std::size_t c = 0; c < cols; ++c)
-    {
-      init_fn(matrix, r, c);
-      ref(r, c) = matrix(r, c);
-    }
-
-  constexpr double epsilon = 1e-10;
-
-  TTS_RELATIVE_EQUAL(matrix.norm(), ref.norm(), epsilon);
-  TTS_RELATIVE_EQUAL(matrix.squaredNorm(), ref.squaredNorm(), epsilon);
-
-  TTS_RELATIVE_EQUAL(matrix.template lpNorm<1>(), ref.lpNorm<1>(), epsilon);
-  TTS_RELATIVE_EQUAL(matrix.template lpNorm<2>(), ref.lpNorm<2>(), epsilon);
-  TTS_RELATIVE_EQUAL(matrix.template lpNorm<-1>(), ref.lpNorm<Eigen::Infinity>(), epsilon);
-}
-
-TTS_CASE("Matrix norm-related operations")
-{
-  std::vector<MatrixDescriptor<rotgen::matrix<double>>> test_matrices = {
-      {3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r + c; }},
-      {2, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 0.0; }},
-      {2, 7, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = -r*r*r - c*c - 1.23; }},
-      {17, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r*c + 0.98; }},
-      {1, 1, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 42.0; }},
-      {10, 11, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = std::tan(r + c); }},
-      {1, 5, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = -1.5*r + 2.56*c + 3.33; }},
-      {7, 1, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r == c ? 1.0 : 0.0); }},
-      {0, 0, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r == c ? 1.0 : 0.0); }},
-  };
-
-  for (const auto& [rows, cols, init_fn] : test_matrices)
-    test_matrix_norms<rotgen::matrix<double>>(rows, cols, init_fn);
-
-};
\ No newline at end of file
diff --git a/test/basic/operators.cpp b/test/basic/operators.cpp
deleted file mode 100644
index a3051ad..0000000
--- a/test/basic/operators.cpp
+++ /dev/null
@@ -1,265 +0,0 @@
-//==================================================================================================
-/*
-  ROTGEN - Runtime Overlay for Eigen
-  Copyright : CODE RECKONS
-  SPDX-License-Identifier: BSL-1.0
-*/
-//==================================================================================================
-#define TTS_MAIN
-#include <rotgen/matrix.hpp>
-#include "tts.hpp"
-
-template <typename MatrixType>
-void test_matrix_scalar_multiplication(std::size_t rows, std::size_t cols, double scalar, auto init_fn)
-{
-  MatrixType a(rows, cols);
-  MatrixType ref(rows, cols);
-
-  for (std::size_t r = 0; r < rows; ++r)
-  {
-    for (std::size_t c = 0; c < cols; ++c)
-    {
-      init_fn(a, r, c);
-      ref(r, c) = a(r, c) * scalar;
-    }
-  }
-
-  TTS_EQUAL(a * scalar, ref);
-  TTS_EQUAL(scalar * a, ref);
-
-  a *= scalar;
-  TTS_EQUAL(a, ref);
-}
-
-template <typename MatrixType>
-void test_matrix_scalar_division(std::size_t rows, std::size_t cols, double scalar, auto init_fn)
-{
-  MatrixType a(rows, cols);
-  MatrixType ref(rows, cols);
-
-  for (std::size_t r = 0; r < rows; ++r)
-  {
-    for (std::size_t c = 0; c < cols; ++c)
-    {
-      init_fn(a, r, c);
-      ref(r, c) = a(r, c) / scalar;
-    }
-  }
-
-  TTS_EQUAL(a / scalar, ref);
-  a /= scalar;
-  TTS_EQUAL(a, ref);
-}
-
-template <typename MatrixType>
-void test_matrix_multiplication(std::size_t n, std::size_t m, std::size_t p, auto a_init_fn, auto b_init_fn)
-{
-  MatrixType a(n, m);
-  MatrixType b(m, p);
-  MatrixType ref(n, p);
-
-  for (std::size_t r = 0; r < n; ++r)
-    for (std::size_t c = 0; c < m; ++c)
-      a_init_fn(a, r, c);
-
-  for (std::size_t r = 0; r < m; ++r)
-    for (std::size_t c = 0; c < p; ++c)
-      b_init_fn(b, r, c);
-
-  for (std::size_t i = 0; i < n; ++i)
-    for (std::size_t j = 0; j < p; ++j)
-    {
-      ref(i, j) = 0;
-      for (std::size_t k = 0; k < m; ++k)
-        ref(i, j) += a(i, k) * b(k, j);
-    }
-
-  TTS_EQUAL(a * b, ref);
-  a *= b;
-  TTS_EQUAL(a, ref);
-}
-
-template <typename MatrixType>
-void test_matrix_addition(std::size_t rows, std::size_t cols, auto a_init_fn, auto b_init_fn)
-{
-  MatrixType a(rows, cols);
-  MatrixType b(rows, cols);
-  MatrixType ref(rows, cols);
-
-  for (std::size_t r = 0; r < rows; ++r)
-  {
-    for (std::size_t c = 0; c < cols; ++c)
-    {
-      a_init_fn(a, r, c);
-      b_init_fn(b, r, c);
-      ref(r, c) = a(r,c) + b(r,c);
-    }
-  }
-
-  TTS_EQUAL(a + b, ref);
-  TTS_EQUAL(b + a, ref);
-  a += b;
-  TTS_EQUAL(a, ref);
-}
-
-template <typename MatrixType>
-void test_matrix_substraction(std::size_t rows, std::size_t cols, auto a_init_fn, auto b_init_fn)
-{
-  MatrixType a(rows, cols);
-  MatrixType b(rows, cols);
-  MatrixType ref(rows, cols);
-  MatrixType a_minus_ref(rows, cols);
-  MatrixType b_minus_ref(rows, cols);
-
-  for (std::size_t r = 0; r < rows; ++r)
-  {
-    for (std::size_t c = 0; c < cols; ++c)
-    {
-      a_init_fn(a, r, c);
-      b_init_fn(b, r, c);
-      ref(r, c) = a(r,c) - b(r,c);
-      a_minus_ref(r, c) = -a(r,c);
-      b_minus_ref(r, c) = -b(r,c);
-    }
-  }
-
-  MatrixType a_unary = -a;
-  MatrixType b_unary = -b;
-
-  TTS_EQUAL(a - b, ref);
-
-  TTS_EQUAL(a_unary, a_minus_ref);
-  TTS_EQUAL(-a, a_minus_ref);
-  TTS_EQUAL(-b, b_minus_ref);
-  TTS_EQUAL(-(-a), a);
-  TTS_EQUAL(-(-b), b);
-
-  a -= b;
-  TTS_EQUAL(a, ref);
-};
-
-TTS_CASE("Check matrix * scalar and scalar * matrix with default values")
-{
-  test_matrix_scalar_multiplication<rotgen::matrix<double>>(2, 2, 10.5,
-      [](auto& a, std::size_t r, std::size_t c)
-      { a(r, c) = (1 + c) + 10 * (1 + r); });
-};
-
-TTS_CASE("Check matrix * scalar with zero scalar multiplication")
-{
-  test_matrix_scalar_multiplication<rotgen::matrix<double>>(3, 2, 0.0,
-      [](auto& a, std::size_t r, std::size_t c)
-      { a(r, c) = 5 * c - r; });
-};
-
-TTS_CASE("Check matrix * scalar with one scalar multiplication")
-{
-  test_matrix_scalar_multiplication<rotgen::matrix<double>>(3, 2, 1,
-      [](auto& a, std::size_t r, std::size_t c)
-      { a(r, c) = 3.3*r - 6; });
-};
-
-TTS_CASE("Check matrix - scalar with negative scalar multiplication")
-{
-  test_matrix_scalar_multiplication<rotgen::matrix<double>>(3, 2, -36.2,
-      [](auto& a, std::size_t r, std::size_t c)
-      { a(r, c) = r * r - c + 3.9; });
-};
-
-TTS_CASE("Check static matrix - scalar with float scalar multiplication")
-{
-  test_matrix_scalar_multiplication<rotgen::matrix<double, 3, 4>>(3, 4, 5.6,
-  [](auto& a, std::size_t r, std::size_t c)
-  { a(r, c) = 1.2*r+4.5*c-6; });
-};
-
-TTS_CASE("Check matrix / scalar with default values")
-{
-  test_matrix_scalar_division<rotgen::matrix<double>>(6, 7, 10.5,
-      [](auto& a, std::size_t r, std::size_t c)
-      { a(r, c) = (1 + c) + 10 * (1 + r); });
-};
-
-TTS_CASE("Check matrix * scalar with one scalar multiplication")
-{
-  test_matrix_scalar_division<rotgen::matrix<double>>(3, 2, 1,
-      [](auto& a, std::size_t r, std::size_t c)
-      { a(r, c) = 3.3*r - 4.5*c + 1.1; });
-};
-
-TTS_CASE("Check matrix - scalar with negative scalar multiplication")
-{
-  test_matrix_scalar_division<rotgen::matrix<double>>(2, 7, -36.2,
-      [](auto& a, std::size_t r, std::size_t c)
-      { a(r, c) = 3.4 * r * r - c + 3.9; });
-};
-
-TTS_CASE("Check static matrix - scalar with float scalar multiplication")
-{
-  test_matrix_scalar_division<rotgen::matrix<double, 3, 4>>(3, 4, 5.6,
-      [](auto& a, std::size_t r, std::size_t c)
-      { a(r, c) = 1.2*r+4.5*c-6; });
-};
-
-TTS_CASE("Matrix multiplication")
-{
-  test_matrix_multiplication<rotgen::matrix<int>>(2, 3, 4,
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = r + c; },
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = r * c; });
-};
-
-TTS_CASE("Matrix multiplication with floats")
-{
-  test_matrix_multiplication<rotgen::matrix<double>>(3, 6, 3,
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 3.4 * r + 5.4 * c*c - 5.2; },
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = r * c + 5 * r - 3.33; });
-};
-
-TTS_CASE("Matrix multiplication with zero matrix")
-{
-  test_matrix_multiplication<rotgen::matrix<double>>(3, 4, 5,
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = (-4) * r + 3*c; },
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 0; });
-};
-
-TTS_CASE("Matrix multiplication with itself")
-{
-  test_matrix_multiplication<rotgen::matrix<double>>(3, 3, 3,
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 3*r - 7*r*c + 14.4; },
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 3*r - 7*r*c + 14.4; });
-};
-
-TTS_CASE("Matrix multiplication with identity matrix")
-{
-  test_matrix_multiplication<rotgen::matrix<double>>(4, 4, 4,
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = r*r*r - c + 12; },
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = (r==c); });
-};
-
-TTS_CASE("Matrix addition")
-{
-  test_matrix_addition<rotgen::matrix<double>>(3, 4,
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 9.9*r*r*r - 6*c -12; },
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 3.1*r + 4.2*c - 12.3; });
-};
-
-TTS_CASE("Matrix addition with zero matrix")
-{
-  test_matrix_addition<rotgen::matrix<double>>(3, 4,
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 2*r*r + c + 3.4; },
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 0; });
-};
-
-TTS_CASE("Matrix subtraction")
-{
-  test_matrix_substraction<rotgen::matrix<double>>(3, 4,
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 6.78*r - 5.2*c - 0.01; },
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 3.1*r + 33.456*c*c; });
-};
-
-TTS_CASE("Matrix subtraction with zero matrix")
-{
-  test_matrix_substraction<rotgen::matrix<double>>(3, 4,
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = r + c*c*c*c - 56.6; },
-      [](auto& mat, std::size_t r, std::size_t c) { mat(r, c) = 0; });
-};
diff --git a/test/basic/static_functions.cpp b/test/basic/static_functions.cpp
deleted file mode 100644
index 4330b63..0000000
--- a/test/basic/static_functions.cpp
+++ /dev/null
@@ -1,96 +0,0 @@
-//==================================================================================================
-/*
-  ROTGEN - Runtime Overlay for Eigen
-  Copyright : CODE RECKONS
-  SPDX-License-Identifier: BSL-1.0
-*/
-//==================================================================================================
-#define TTS_MAIN
-#include <rotgen/matrix.hpp>
-#include "tts.hpp"
-
-template <typename MatrixType>
-void test_zero(const MatrixType& matrix, std::size_t rows, std::size_t cols)
-{
-  for(std::size_t r=0;r<rows;++r)
-    for(std::size_t c=0;c<cols;++c)
-      TTS_EQUAL(matrix(r, c), 0);
-}
-
-template <typename MatrixType>
-void test_constant(const MatrixType& matrix, std::size_t rows, std::size_t cols, double constant)
-{
-  for(std::size_t r=0;r<rows;++r)
-    for(std::size_t c=0;c<cols;++c)
-      TTS_EQUAL(matrix(r, c), constant);
-
-  constexpr double epsilon = 1e-10;
-  TTS_RELATIVE_EQUAL(rows * cols * constant, matrix.sum(), epsilon);
-}
-
-template <typename MatrixType>
-void test_random(const MatrixType& matrix, std::size_t rows, std::size_t cols)
-{
-  for(std::size_t r=0;r<rows;++r)
-    for(std::size_t c=0;c<cols;++c)
-    {
-      TTS_GREATER_EQUAL(matrix(r, c), -1.0);
-      TTS_LESS_EQUAL(matrix(r, c), 1.0);
-    }
-}
-
-template <typename MatrixType>
-void test_identity(const MatrixType& matrix, std::size_t rows, std::size_t cols)
-{
-  for(std::size_t r=0;r<rows;++r)
-    for(std::size_t c=0;c<cols;++c)
-      TTS_EQUAL(matrix(r, c), r==c ? 1 : 0);
-
-
-  std::size_t minimum = std::min(rows, cols);
-  TTS_EQUAL(matrix.trace(), minimum);
-}
-
-TTS_CASE("Test zero")
-{
-  test_zero(rotgen::matrix<double, 3, 4>::Zero(), 3, 4);
-  test_zero(rotgen::matrix<double, 1, 1>::Zero(), 1, 1);
-  test_zero(rotgen::matrix<double, 10, 10>::Zero(), 10, 10);
-  test_zero(rotgen::matrix<double>::Zero(3, 4), 3, 4);
-  test_zero(rotgen::matrix<double, 7, 5>::Zero(7, 5), 7, 5);
-  test_zero(rotgen::matrix<double, 9>::Zero(9, 3), 9, 3);
-  test_zero(rotgen::matrix<double, -1, 3>::Zero(2, 3), 2, 3);
-};
-
-TTS_CASE("Test constant")
-{
-  test_constant(rotgen::matrix<double, 3, 8>::Constant(5.12), 3, 8, 5.12);
-  test_constant(rotgen::matrix<double, 1, 1>::Constant(2.2), 1, 1, 2.2);
-  test_constant(rotgen::matrix<double, 11, 12>::Constant(13), 11, 12, 13);
-  test_constant(rotgen::matrix<double>::Constant(2, 7, 5.6), 2, 7, 5.6);
-  test_constant(rotgen::matrix<double, 2, 2>::Constant(2, 2, 2.0), 2, 2, 2.0);
-  test_constant(rotgen::matrix<double, 9>::Constant(9, 3, 1.1), 9, 3, 1.1);
-  test_constant(rotgen::matrix<double, -1, 9>::Constant(5, 9, 42), 5, 9, 42);
-};
-
-TTS_CASE("Test random")
-{
-  test_random(rotgen::matrix<double, 2, 3>::Random(), 2, 3);
-  test_random(rotgen::matrix<double, 1, 1>::Random(), 1, 1);
-  test_random(rotgen::matrix<double, 11, 17>::Random(), 11, 17);
-  test_random(rotgen::matrix<double>::Random(7, 3), 7, 3);
-  test_random(rotgen::matrix<double, 2, 2>::Random(2, 2), 2, 2);
-  test_random(rotgen::matrix<double, 4, -1>::Random(4, 3), 4, 3);
-  test_random(rotgen::matrix<double, -1, 5>::Random(5, 5), 5, 5);
-};
-
-TTS_CASE("Test identity")
-{
-  test_identity(rotgen::matrix<double, 4, 5>::Identity(), 4, 5);
-  test_identity(rotgen::matrix<double, 1, 1>::Identity(), 1, 1);
-  test_identity(rotgen::matrix<double, 21, 3>::Identity(), 21, 3);
-  test_identity(rotgen::matrix<double>::Identity(2, 7), 2, 7);
-  test_identity(rotgen::matrix<double, 2, 2>::Identity(2, 2), 2, 2);
-  test_identity(rotgen::matrix<double, 3>::Identity(3, 3), 3, 3);
-  test_identity(rotgen::matrix<double, -1, 11>::Identity(5, 11), 5, 11);
-};
diff --git a/test/unit/basic/arithmetic_functions.cpp b/test/unit/basic/arithmetic_functions.cpp
new file mode 100644
index 0000000..5b50daf
--- /dev/null
+++ b/test/unit/basic/arithmetic_functions.cpp
@@ -0,0 +1,114 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#include "unit/tests.hpp"
+#include <rotgen/matrix.hpp>
+#include <Eigen/Dense>
+
+template<typename MatrixType>
+void test_matrix_unary_ops(std::size_t rows, std::size_t cols, auto const &init_fn)
+{
+  MatrixType original(rows, cols);
+  MatrixType transposed_matrix(cols, rows);
+
+  for (std::size_t r = 0; r < rows; ++r)
+    for (std::size_t c = 0; c < cols; ++c)
+      original(r,c) = init_fn(r, c);
+
+  for (std::size_t r = 0; r < rows; ++r)
+    for (std::size_t c = 0; c < cols; ++c)
+      transposed_matrix(c, r) = original(r, c);
+
+  TTS_EQUAL(original.transpose(), transposed_matrix);
+  TTS_EQUAL(original.conjugate(), original);
+  TTS_EQUAL(original.adjoint(), transposed_matrix);
+
+  original.transposeInPlace();
+  TTS_EQUAL(original, transposed_matrix);
+
+  original.transposeInPlace();
+  original.adjointInPlace();
+  TTS_EQUAL(original, transposed_matrix);
+}
+
+TTS_CASE_TPL("Test transpotion related operations", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  std::vector<rotgen::tests::matrix_descriptor> test_matrices =
+  {
+    {3, 3, [](auto r, auto c) { return r + 3 * c - 2.5; }},
+    {4, 9, [](auto r, auto c) { return r*r + 3.12 * c + 6.87; }},
+    {2, 7, [](auto r, auto c) { return 1.1 * (r - c); }},
+    {1, 5, [](auto  , auto  ) { return 9.99; }},
+    {4, 2, [](auto  , auto  ) { return 0.0; }},
+    {3, 3, [](auto r, auto c) { return (r == c) ? 1.0 : 0.0; }},
+    {2, 2, [](auto r, auto c) { return (r + c) * 1e-10; }},
+    {2, 2, [](auto r, auto  ) { return (r + 1) * 1e+10; }}
+  };
+
+  for (const auto& [rows, cols, fn] : test_matrices)
+  {
+    test_matrix_unary_ops<rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>>(rows, cols, fn);
+  }
+};
+
+template <typename MatrixType>
+void test_matrix_reductions(std::size_t rows, std::size_t cols, auto const& init_fn)
+{
+  using EigenMatrix = Eigen::Matrix<typename MatrixType::value_type, Eigen::Dynamic, Eigen::Dynamic>;
+
+  MatrixType original(rows, cols);
+  EigenMatrix ref(rows, cols);
+
+  for (std::size_t r = 0; r < rows; ++r)
+    for (std::size_t c = 0; c < cols; ++c)
+      ref(r, c) = original(r,c) = init_fn(r, c);
+
+  TTS_EQUAL(original.sum(), ref.sum());
+  TTS_EQUAL(original.prod(), ref.prod());
+  TTS_EQUAL(original.mean(), ref.mean());
+  TTS_EQUAL(original.maxCoeff(), ref.maxCoeff());
+  TTS_EQUAL(original.minCoeff(), ref.minCoeff());
+  TTS_EQUAL(original.trace(), ref.trace());
+
+  std::ptrdiff_t row, col, ref_row, ref_col;
+
+  double cmin = original.minCoeff(&row, &col);
+  double rmin = ref.minCoeff(&ref_row, &ref_col);
+
+  TTS_EQUAL(cmin, rmin);
+  TTS_EQUAL(row, ref_row);
+  TTS_EQUAL(col, ref_col);
+
+  double cmax = original.maxCoeff(&row, &col);
+  double rmax = ref.maxCoeff(&ref_row, &ref_col);
+
+  TTS_EQUAL(cmax, rmax);
+  TTS_EQUAL(row, ref_row);
+  TTS_EQUAL(col, ref_col);
+}
+
+TTS_CASE_TPL("Test reductions", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  std::vector<rotgen::tests::matrix_descriptor> test_matrices =
+  {
+    {3, 3, [](auto r, auto c) {return r + c; }},
+    {3, 3, [](auto  , auto  ) {return 0.0; }},
+    {2, 4, [](auto r, auto c) {return -r -c*c - 1234; }},
+    {4, 4, [](auto  , auto  ) {return 7.0; }},
+    {1, 1, [](auto  , auto  ) {return 42.0; }},
+    {4, 2, [](auto r, auto c) {return std::sin(r + c); }},
+    {1, 5, [](auto r, auto c) {return -1.5 * r + 2.56 * c; }},
+    {5, 7, [](auto r, auto c) {return (r == c ? 1.0 : 0.0); }}
+  };
+
+  for (const auto& [rows, cols, fn] : test_matrices)
+  {
+    test_matrix_reductions<rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>>(rows, cols, fn);
+  }
+};
diff --git a/test/unit/basic/constructors.cpp b/test/unit/basic/constructors.cpp
new file mode 100644
index 0000000..5506c34
--- /dev/null
+++ b/test/unit/basic/constructors.cpp
@@ -0,0 +1,187 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#include "unit/tests.hpp"
+#include <rotgen/matrix.hpp>
+
+TTS_CASE_TPL("Default matrix dynamic constructor", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> matrix;
+
+  TTS_EQUAL(matrix.rows(), 0ULL);
+  TTS_EQUAL(matrix.cols(), 0ULL);
+};
+
+TTS_CASE_TPL("Default matrix static constructor", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,4,9,O::value> matrix;
+
+  TTS_EQUAL(matrix.rows(), 4ULL);
+  TTS_EQUAL(matrix.cols(), 9ULL);
+};
+
+TTS_CASE_TPL("Dynamic matrix constructor with row and columns", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> matrix(10, 5);
+
+  TTS_EQUAL(matrix.rows(), 10ULL);
+  TTS_EQUAL(matrix.cols(), 5ULL);
+};
+
+TTS_CASE_TPL("Static matrix constructor with row and columns", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,6,11,O::value> matrix(6, 11);
+
+  TTS_EQUAL(matrix.rows(), 6ULL);
+  TTS_EQUAL(matrix.cols(), 11ULL);
+};
+
+TTS_CASE_TPL("Copy constructor produces identical but independent matrix", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> a(3, 3);
+
+  for(std::size_t r=0;r<a.rows();r++)
+    for(std::size_t c=0;c<a.cols();c++)
+      a(r,c)   = r + c;
+
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> b = a;
+
+  TTS_EQUAL(b.rows(), a.rows());
+  TTS_EQUAL(b.cols(), a.cols());
+
+  for(std::size_t r=0;r<a.rows();r++)
+    for(std::size_t c=0;c<a.cols();c++)
+      TTS_EQUAL(b(r,c), a(r,c));
+
+  TTS_EQUAL(b(0, 0), 0.0);
+  TTS_EQUAL(b(1, 0), 1.0);
+  TTS_EQUAL(b(2, 1), 3.0);
+
+  a(0,0) = 42.0;
+  TTS_NOT_EQUAL(b(0,0), a(0,0));
+};
+
+TTS_CASE_TPL("Copy constructor on default matrix", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> a;
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> b = a;
+  TTS_EQUAL(b.rows(), 0ULL);
+  TTS_EQUAL(b.cols(), 0ULL);
+};
+
+TTS_CASE_TPL("Copy constructor from const matrix", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  const rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> a(2, 2);
+  auto b = a;
+  TTS_EQUAL(b.rows(), 2ULL);
+  TTS_EQUAL(b.cols(), 2ULL);
+};
+
+TTS_CASE_TPL("Copy constructor on static matrix", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+   rotgen::matrix<T, 2, 5> a;
+   rotgen::matrix<T, 2, 5> b = a;
+  TTS_EQUAL(b.rows(), 2ULL);
+  TTS_EQUAL(b.cols(), 5ULL);
+};
+
+/*
+TTS_CASE_TPL("Copy constructor on static/dynamic matrix", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T, 11, 4,O::value> a;
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> b = a;
+  TTS_EQUAL(b.rows(), 11);
+  TTS_EQUAL(b.cols(), 4);
+};
+
+TTS_CASE_TPL("Copy constructor on dynamic/static matrix", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> a(5, 7);
+  rotgen::matrix<T, 5, 7,O::value> b = a;
+  TTS_EQUAL(b.rows(), 5);
+  TTS_EQUAL(b.cols(), 7);
+};*/
+
+TTS_CASE_TPL("Move constructor transfers contents", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> a(3, 3);
+  a(1,1) = 7;
+  auto ptr = a.data();
+
+ rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> b = std::move(a);
+
+  TTS_EQUAL(b(1,1), 7);
+  TTS_EQUAL(b.rows(), 3ULL);
+  TTS_EQUAL(b.cols(), 3ULL);
+  TTS_EXPECT(b.data() == ptr);
+};
+
+TTS_CASE_TPL("Move constructor from Rvalue", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> b = rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>(2, 2);
+  TTS_EQUAL(b.rows(), 2ULL);
+  TTS_EQUAL(b.cols(), 2ULL);
+};
+
+TTS_CASE_TPL("Constructor from Initializer list", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,1,1,O::value> b1 = {3.5};
+  TTS_EQUAL(b1.rows(), 1ULL);
+  TTS_EQUAL(b1.cols(), 1ULL);
+  TTS_EQUAL(b1(0)    , 3.5);
+
+  rotgen::matrix<T,5,1,O::value> b9 = {0.25, 0.5, 1, 2, 4};
+  TTS_EQUAL(b9.rows(), 5ULL);
+  TTS_EQUAL(b9.cols(), 1ULL);
+
+  T i = 0.25;
+  for(std::size_t r=0;r<b9.rows();++r) {
+    TTS_EQUAL(b9(r,0), i);
+    i *= 2;
+  }
+
+  rotgen::matrix<T,1,3,O::value> b13 = {1.2,2.3,3.4};
+  TTS_EQUAL(b13.rows(), 1ULL);
+  TTS_EQUAL(b13.cols(), 3ULL);
+  TTS_EQUAL(b13(0)    , T(1.2));
+  TTS_EQUAL(b13(1)    , T(2.3));
+  TTS_EQUAL(b13(2)    , T(3.4));
+};
+
+TTS_CASE_TPL("Constructor from Initializer list of rows", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> b1 = {{3.5}};
+  TTS_EQUAL(b1.rows(), 1ULL);
+  TTS_EQUAL(b1.cols(), 1ULL);
+  TTS_EQUAL(b1(0,0)  , T(3.5));
+
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> b23  = { {1.2,2.3,3.4}
+    , {10,200,3000}
+  };
+  TTS_EQUAL(b23.rows(), 2ULL);
+  TTS_EQUAL(b23.cols(), 3ULL);
+  TTS_EQUAL(b23(0,0)    , T(1.2));
+  TTS_EQUAL(b23(0,1)    , T(2.3));
+  TTS_EQUAL(b23(0,2)    , T(3.4));
+  TTS_EQUAL(b23(1,0)    , T(10));
+  TTS_EQUAL(b23(1,1)    , T(200));
+  TTS_EQUAL(b23(1,2)    , T(3000));
+};
\ No newline at end of file
diff --git a/test/basic/functions.cpp b/test/unit/basic/functions.cpp
similarity index 58%
rename from test/basic/functions.cpp
rename to test/unit/basic/functions.cpp
index 7d5c3a4..5de37fa 100644
--- a/test/basic/functions.cpp
+++ b/test/unit/basic/functions.cpp
@@ -5,16 +5,16 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#define TTS_MAIN
+#include "unit/tests.hpp"
 #include <rotgen/matrix.hpp>
-#include "tts.hpp"
 
-TTS_CASE("Function size")
+TTS_CASE_TPL("Function size", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
 {
-  rotgen::matrix<double> empty_matrix;
-  rotgen::matrix<double> matrix(3,4);
-  rotgen::matrix<double> row_vector(9,1);
-  rotgen::matrix<double> column_vector(1,5);
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> empty_matrix;
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> matrix(3,4);
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> row_vector(9,1);
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> column_vector(1,5);
 
   TTS_EQUAL(empty_matrix.size(), 0ULL);
   TTS_EQUAL(matrix.size(), 12ULL);
@@ -22,9 +22,10 @@ TTS_CASE("Function size")
   TTS_EQUAL(column_vector.size(), 5ULL);
 };
 
-TTS_CASE("Resizing dynamic matrix")
+TTS_CASE_TPL("Resizing dynamic matrix", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
 {
-  rotgen::matrix<double> a(2, 3);
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> a(2, 3);
 
   for(std::size_t r=0;r<a.rows();++r)
     for(std::size_t c=0;c<a.cols();++c)
@@ -41,21 +42,19 @@ TTS_CASE("Resizing dynamic matrix")
   a.resize(2,2);
   TTS_EQUAL(a.rows(), 2ULL);
   TTS_EQUAL(a.cols(), 2ULL);
-
-  for(std::size_t r=0;r<a.rows();++r)
-    for(std::size_t c=0;c<a.cols();++c)
-      TTS_EQUAL(a(r,c),0);
 };
 
-TTS_CASE("Resizing static matrix")
+TTS_CASE_TPL("Resizing static matrix", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
 {
   rotgen::matrix<double, 3, 3> a;
   TTS_EXPECT_NOT_COMPILES(a, { a.resize(4, 5); });
 };
 
-TTS_CASE("Dynamix matrix conservative resizing")
+TTS_CASE_TPL("Dynamix matrix conservative resizing", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
 {
-  rotgen::matrix<double> a(2, 3);
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> a(2, 3);
 
   int i = 1;
   for(std::size_t r=0;r<a.rows();++r)
@@ -74,8 +73,8 @@ TTS_CASE("Dynamix matrix conservative resizing")
   a.conservativeResize(3, 2);
   TTS_EQUAL(a.rows(), 3ULL);
   TTS_EQUAL(a.cols(), 2ULL);
-  int expected[3][2] = { {1, 2}, {4, 5}, {0, 0} };
-  for(std::size_t r = 0; r < 3; ++r)
+  int expected[3][2] = { {1, 2}, {4, 5} };
+  for(std::size_t r = 0; r < 2; ++r)
     for(std::size_t c = 0; c < 2; ++c)
       TTS_EQUAL(a(r, c), expected[r][c]);
 
@@ -104,20 +103,19 @@ TTS_CASE("Dynamix matrix conservative resizing")
   a.conservativeResize(3, 3);
   TTS_EQUAL(a.rows(), 3ULL);
   TTS_EQUAL(a.cols(), 3ULL);
-  for(std::size_t r=0;r<a.rows();++r)
-    for(std::size_t c=0;c<a.cols();++c)
-      TTS_EQUAL(a(r,c),0);
 };
 
-TTS_CASE("Static matrix conservative resizing")
+TTS_CASE_TPL("Static matrix conservative resizing", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
 {
   rotgen::matrix<double, 3, 3> a;
   TTS_EXPECT_NOT_COMPILES(a, { a.conservativeResize(4, 5); });
 };
 
-TTS_CASE("Test coefficient accessors")
+TTS_CASE_TPL("Test coefficient accessors", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
 {
-  rotgen::matrix<double> a(3, 5);
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> a(3, 5);
 
   for(std::size_t r=0;r<a.rows();++r)
     for(std::size_t c=0;c<a.cols();++c)
@@ -132,36 +130,35 @@ TTS_CASE("Test coefficient accessors")
   a(1, 1) = 42;
   TTS_EQUAL(a(1,1), 42);
 
-  double& ref = a(2, 2);
+  T& ref = a(2, 2);
   ref = 17;
   assert(a(2, 2) == 17);
 };
 
-TTS_CASE("Test one index coefficient accessors")
+TTS_CASE_TPL("Test one index coefficient accessors", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
 {
-  rotgen::matrix<double> a(2, 4);
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> a(2, 4);
 
   int i = 1;
-  for(std::size_t c=0;c<a.cols();++c)
+  if constexpr(!O::value)
+  {
+    for(std::size_t c=0;c<a.cols();++c)
+      for(std::size_t r=0;r<a.rows();++r)
+        a(r, c) = i++;
+  }
+  else
+  {
     for(std::size_t r=0;r<a.rows();++r)
-      a(r, c) = i++;
+      for(std::size_t c=0;c<a.cols();++c)
+        a(r, c) = i++;
+  }
 
   i = 1;
   for(std::size_t r=0;r<a.rows()*a.cols();++r)
-    TTS_EQUAL(a(r), i++);
+    TTS_EQUAL(a(r), i++) << a;
 
-  double& ref = a(2);
-  ref = 999.1;
-  assert(a(2) == 999.1);
-
-  rotgen::matrix<double, 2, 3> b (2, 3);
-
-  i = 1;
-  for(std::size_t c=0;c<b.cols()*b.rows(); ++c)
-    a(c) = i++;
-
-  i = 1;
-  for(std::size_t c=0;c<b.cols();++c)
-    for(std::size_t r=0;r<b.rows();++r)
-      TTS_EQUAL(a(r, c), i++);
+  T& ref = a(2);
+  ref = 999.5;
+  TTS_EQUAL(a(2), 999.5);
 };
\ No newline at end of file
diff --git a/test/basic/io.cpp b/test/unit/basic/io.cpp
similarity index 57%
rename from test/basic/io.cpp
rename to test/unit/basic/io.cpp
index e1a0af8..0649ea2 100644
--- a/test/basic/io.cpp
+++ b/test/unit/basic/io.cpp
@@ -5,22 +5,16 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#define TTS_MAIN
+#include "unit/tests.hpp"
 #include <rotgen/matrix.hpp>
 #include <sstream>
-#include "tts.hpp"
 
-TTS_CASE("Sample test")
+TTS_CASE_TPL("Sample test", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
 {
-  rotgen::matrix<double,5,5> x;
+  rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value> x({ {1,2} , {3,4} });
   std::ostringstream os;
   os << x;
 
-  std::string ref = "0 0 0 0 0\n"
-                    "0 0 0 0 0\n"
-                    "0 0 0 0 0\n"
-                    "0 0 0 0 0\n"
-                    "0 0 0 0 0";
-
-  TTS_EQUAL(os.str(), ref);
+  TTS_EQUAL(os.str(), std::string{"1 2\n3 4"});
 };
\ No newline at end of file
diff --git a/test/unit/basic/norms.cpp b/test/unit/basic/norms.cpp
new file mode 100644
index 0000000..189084b
--- /dev/null
+++ b/test/unit/basic/norms.cpp
@@ -0,0 +1,47 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#include "unit/tests.hpp"
+#include <rotgen/matrix.hpp>
+#include <Eigen/Dense>
+
+TTS_CASE_TPL("Matrix norm-related operations", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  std::vector<rotgen::tests::matrix_descriptor> test_matrices =
+  {
+    { 3,  3, [](int r, int c) { return  r + c; }},
+    { 2,  3, [](int  , int  ) { return  T{0}; }},
+    { 2,  7, [](int r, int c) { return  -r*r*r - c*c - T(1.23); }},
+    {17,  3, [](int r, int c) { return  r*c + T(0.98); }},
+    { 1,  1, [](int  , int  ) { return  42; }},
+    {10, 11, [](int r, int c) { return  std::tan(T(r + c)); }},
+    { 1,  5, [](int r, int c) { return  T(-1.5)*r + T(2.56)*c + T(3.33); }},
+    { 7,  1, [](int r, int c) { return  r == c ? T(1) : T(0); }},
+    { 0,  0, [](int r, int c) { return  r == c ? T(1) : T(0); }},
+  };
+
+  for (const auto& [rows, cols, fn] : test_matrices)
+  {
+    rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>  matrix(rows, cols);
+    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic,O::value>   ref(rows, cols);
+
+    for (std::size_t r = 0; r < rows; ++r)
+    {
+      for (std::size_t c = 0; c < cols; ++c)
+      {
+        ref(r, c) =  matrix(r, c)  = fn(r,c);
+      }
+    }
+
+    TTS_EQUAL(matrix.norm(), ref.norm());
+    TTS_EQUAL(matrix.squaredNorm(), ref.squaredNorm());
+    TTS_EQUAL(matrix.template lpNorm<1>()               , ref.template lpNorm<1>());
+    TTS_EQUAL(matrix.template lpNorm<2>()               , ref.template lpNorm<2>());
+    TTS_EQUAL(matrix.template lpNorm<rotgen::Infinity>(), ref.template lpNorm<Eigen::Infinity>());
+  }
+};
\ No newline at end of file
diff --git a/test/unit/basic/operators.cpp b/test/unit/basic/operators.cpp
new file mode 100644
index 0000000..fd1930e
--- /dev/null
+++ b/test/unit/basic/operators.cpp
@@ -0,0 +1,193 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#include "unit/tests.hpp"
+#include <rotgen/matrix.hpp>
+
+template <typename MatrixType>
+void test_matrix_operations(std::size_t rows, std::size_t cols, auto a_init_fn, auto b_init_fn, auto ops, auto self_ops)
+{
+  MatrixType a(rows, cols);
+  MatrixType b(rows, cols);
+  MatrixType ref(rows, cols);
+
+  for (std::size_t r = 0; r < rows; ++r)
+  {
+    for (std::size_t c = 0; c < cols; ++c)
+    {
+      a(r,c) = a_init_fn(r, c);
+      b(r,c) = b_init_fn(r, c);
+      ref(r, c) = ops(a(r,c),b(r,c));
+    }
+  }
+
+  TTS_EQUAL(ops(a, b), ref);
+  self_ops(a,b);
+  TTS_EQUAL(a, ref);
+}
+
+template <typename MatrixType>
+void test_scalar_operations(std::size_t rows, std::size_t cols, auto a_init_fn, auto s, auto ops, auto self_ops)
+{
+  MatrixType a(rows, cols);
+  MatrixType ref(rows, cols);
+
+  for (std::size_t r = 0; r < rows; ++r)
+  {
+    for (std::size_t c = 0; c < cols; ++c)
+    {
+      a(r,c) = a_init_fn(r, c);
+      ref(r, c) = ops(a(r,c),s);
+    }
+  }
+
+  TTS_EQUAL(ops(a, s), ref);
+  self_ops(a,s);
+  TTS_EQUAL(a, ref);
+}
+
+template <typename MatrixType>
+void test_scalar_multiplications(std::size_t rows, std::size_t cols, auto fn, auto s)
+{
+  MatrixType a(rows, cols);
+  MatrixType ref(rows, cols);
+
+  for (std::size_t r = 0; r < rows; ++r)
+  {
+    for (std::size_t c = 0; c < cols; ++c)
+    {
+      a(r,c) = fn(r, c);
+      ref(r, c) = a(r,c) * s;
+    }
+  }
+
+  TTS_EQUAL(a * s, ref);
+  TTS_EQUAL(s * a, ref);
+  a *= s;
+  TTS_EQUAL(a, ref);
+}
+
+template <typename MatrixType>
+void test_matrix_multiplication(std::size_t rows, std::size_t cols, auto a_init_fn, auto b_init_fn)
+{
+  MatrixType a(rows, cols);
+  MatrixType b(cols, rows);
+  MatrixType ref(rows, rows);
+
+  for (std::size_t r = 0; r < a.rows(); ++r)
+    for (std::size_t c = 0; c < a.cols(); ++c)
+      a(r,c) = a_init_fn(r, c);
+
+  for (std::size_t r = 0; r < b.rows(); ++r)
+    for (std::size_t c = 0; c < b.cols(); ++c)
+      b(r,c) = b_init_fn(r, c);
+
+  for (std::size_t i = 0; i < a.rows(); ++i)
+  {
+    for (std::size_t j = 0; j < b.cols(); ++j)
+    {
+      ref(i, j) = 0;
+      for (std::size_t k = 0; k < a.cols(); ++k)
+        ref(i, j) += a(i, k) * b(k, j);
+    }
+  }
+
+  TTS_EQUAL(a * b, ref);
+}
+
+// Basic initializers
+inline constexpr auto init_a = [](auto  r, auto c) { return 9.9*r*r*r - 6*c -12;  };
+inline constexpr auto init_b = [](auto  r, auto c) { return 3.1*r + 4.2*c - 12.3; };
+inline constexpr auto init_0 = [](auto   , auto  ) { return 0; };
+
+TTS_CASE_TPL("Check matrix addition", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  using mat_t = rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>;
+  auto op     = [](auto  a, auto b) { return a + b;   };
+  auto s_op   = [](auto& a, auto b) { return a += b;  };
+
+  test_matrix_operations<mat_t>(1 , 1, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(3 , 5, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(5 , 3, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(5 , 5, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(5 , 5, init_b, init_a, op, s_op);
+  test_matrix_operations<mat_t>(10, 1, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(1 ,10, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(5 , 5, init_0, init_b, op, s_op);
+  test_matrix_operations<mat_t>(5 , 5, init_a, init_0, op, s_op);
+};
+
+TTS_CASE_TPL("Check matrix substraction", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  using mat_t = rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>;
+  auto op     = [](auto  a, auto b) { return a - b;   };
+  auto s_op   = [](auto& a, auto b) { return a -= b;  };
+
+  test_matrix_operations<mat_t>(1 , 1, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(3 , 5, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(5 , 3, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(5 , 5, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(5 , 5, init_b, init_a, op, s_op);
+  test_matrix_operations<mat_t>(10, 1, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(1 ,10, init_a, init_b, op, s_op);
+  test_matrix_operations<mat_t>(5 , 5, init_0, init_b, op, s_op);
+  test_matrix_operations<mat_t>(5 , 5, init_a, init_0, op, s_op);
+};
+
+TTS_CASE_TPL("Check matrix multiplications", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  using mat_t = rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>;
+  auto init_id = [](auto  r, auto c) { return r == c ? 1 : 0; };
+
+  test_matrix_multiplication<mat_t>(1 , 1, init_a , init_b );
+  test_matrix_multiplication<mat_t>(3 , 5, init_a , init_b );
+  test_matrix_multiplication<mat_t>(5 , 3, init_a , init_b );
+  test_matrix_multiplication<mat_t>(5 , 5, init_a , init_b );
+  test_matrix_multiplication<mat_t>(5 , 5, init_b , init_a );
+  test_matrix_multiplication<mat_t>(5 , 5, init_a , init_a );
+  test_matrix_multiplication<mat_t>(5 , 5, init_a , init_id);
+  test_matrix_multiplication<mat_t>(5 , 5, init_id, init_a );
+  test_matrix_multiplication<mat_t>(10, 1, init_a , init_b );
+  test_matrix_multiplication<mat_t>(1 ,10, init_a , init_b );
+  test_matrix_multiplication<mat_t>(5 , 5, init_0 , init_b );
+  test_matrix_multiplication<mat_t>(5 , 5, init_a , init_0 );
+};
+
+TTS_CASE_TPL("Check matrix multiplication with scalar", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  using mat_t = rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>;
+
+  test_scalar_multiplications<mat_t>(1 , 1, init_a, T{ 3.5});
+  test_scalar_multiplications<mat_t>(3 , 5, init_a, T{-2.5});
+  test_scalar_multiplications<mat_t>(5 , 3, init_a, T{ 4. });
+  test_scalar_multiplications<mat_t>(5 , 5, init_a, T{-5. });
+  test_scalar_multiplications<mat_t>(5 , 5, init_a, T{ 1. });
+  test_scalar_multiplications<mat_t>(5 , 5, init_a, T{ 6. });
+  test_scalar_multiplications<mat_t>(10, 1, init_a, T{ 10.});
+  test_scalar_multiplications<mat_t>(1 ,10, init_a, T{-0.5});
+};
+
+TTS_CASE_TPL("Check matrix division with scalar", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  using mat_t = rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic,O::value>;
+  auto op     = [](auto  a, auto b) { return a / b;   };
+  auto s_op   = [](auto& a, auto b) { return a /= b;  };
+
+  test_scalar_operations<mat_t>(1 , 1, init_a, T{ 3.5}, op, s_op);
+  test_scalar_operations<mat_t>(3 , 5, init_a, T{-2.5}, op, s_op);
+  test_scalar_operations<mat_t>(5 , 3, init_a, T{ 4. }, op, s_op);
+  test_scalar_operations<mat_t>(5 , 5, init_a, T{-5. }, op, s_op);
+  test_scalar_operations<mat_t>(5 , 5, init_a, T{ 1. }, op, s_op);
+  test_scalar_operations<mat_t>(5 , 5, init_a, T{ 6. }, op, s_op);
+  test_scalar_operations<mat_t>(10, 1, init_a, T{ 10.}, op, s_op);
+  test_scalar_operations<mat_t>(1 ,10, init_a, T{-0.5}, op, s_op);
+};
diff --git a/test/unit/basic/static_functions.cpp b/test/unit/basic/static_functions.cpp
new file mode 100644
index 0000000..deb95e4
--- /dev/null
+++ b/test/unit/basic/static_functions.cpp
@@ -0,0 +1,81 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#include "unit/tests.hpp"
+#include <rotgen/matrix.hpp>
+
+void test_value(const auto& matrix, std::size_t rows, std::size_t cols, auto constant)
+{
+  for(std::size_t r=0;r<rows;++r)
+    for(std::size_t c=0;c<cols;++c)
+      TTS_EQUAL(matrix(r, c), constant);
+}
+
+void test_random(const auto& matrix, std::size_t rows, std::size_t cols)
+{
+  for(std::size_t r=0;r<rows;++r)
+    for(std::size_t c=0;c<cols;++c)
+    {
+      TTS_GREATER_EQUAL(matrix(r, c), -1.0);
+      TTS_LESS_EQUAL(matrix(r, c), 1.0);
+    }
+}
+
+void test_identity(const auto& matrix, std::size_t rows, std::size_t cols)
+{
+  for(std::size_t r=0;r<rows;++r)
+    for(std::size_t c=0;c<cols;++c)
+      TTS_EQUAL(matrix(r, c), r==c ? 1 : 0);
+}
+
+TTS_CASE_TPL("Test zero", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  test_value(rotgen::matrix<T, 3, 4, O::value>::Zero(), 3, 4, 0);
+  test_value(rotgen::matrix<T, 1, 1, O::value>::Zero(), 1, 1, 0);
+  test_value(rotgen::matrix<T, 10, 10, O::value>::Zero(), 10, 10, 0);
+  test_value(rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic, O::value>::Zero(3, 4), 3, 4, 0);
+  test_value(rotgen::matrix<T, 7, 5, O::value>::Zero(7, 5), 7, 5, 0);
+  test_value(rotgen::matrix<T, 9,rotgen::Dynamic, O::value>::Zero(9, 3), 9, 3, 0);
+  test_value(rotgen::matrix<T,rotgen::Dynamic, 3, O::value>::Zero(2, 3), 2, 3, 0);
+};
+
+TTS_CASE_TPL("Test constant", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  test_value(rotgen::matrix<T, 3, 8, O::value>::Constant(5.12), 3, 8, T(5.12));
+  test_value(rotgen::matrix<T, 1, 1, O::value>::Constant(2.2), 1, 1, T(2.2));
+  test_value(rotgen::matrix<T, 11, 12, O::value>::Constant(13), 11, 12, T(13));
+  test_value(rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic, O::value>::Constant(2, 7, 5.6), 2, 7, T(5.6));
+  test_value(rotgen::matrix<T, 2, 2, O::value>::Constant(2, 2, 2.0), 2, 2, T(2.0));
+  test_value(rotgen::matrix<T, 9,rotgen::Dynamic, O::value>::Constant(9, 3, 1.1), 9, 3, T(1.1));
+  test_value(rotgen::matrix<T,rotgen::Dynamic, 9, O::value>::Constant(5, 9, 42), 5, 9,T(42));
+};
+
+TTS_CASE_TPL("Test random", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  test_random(rotgen::matrix<T, 2, 3, O::value>::Random(), 2, 3);
+  test_random(rotgen::matrix<T, 1, 1, O::value>::Random(), 1, 1);
+  test_random(rotgen::matrix<T, 11, 17, O::value>::Random(), 11, 17);
+  test_random(rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic, O::value>::Random(7, 3), 7, 3);
+  test_random(rotgen::matrix<T, 2, 2, O::value>::Random(2, 2), 2, 2);
+  test_random(rotgen::matrix<T, 4,rotgen::Dynamic, O::value>::Random(4, 3), 4, 3);
+  test_random(rotgen::matrix<T,rotgen::Dynamic, 5, O::value>::Random(5, 5), 5, 5);
+};
+
+TTS_CASE_TPL("Test identity", rotgen::tests::types)
+<typename T, typename O>( tts::type< tts::types<T,O>> )
+{
+  test_identity(rotgen::matrix<T, 4, 5, O::value>::Identity(), 4, 5);
+  test_identity(rotgen::matrix<T, 1, 1, O::value>::Identity(), 1, 1);
+  test_identity(rotgen::matrix<T, 21, 3, O::value>::Identity(), 21, 3);
+  test_identity(rotgen::matrix<T,rotgen::Dynamic,rotgen::Dynamic, O::value>::Identity(2, 7), 2, 7);
+  test_identity(rotgen::matrix<T, 2, 2, O::value>::Identity(2, 2), 2, 2);
+  test_identity(rotgen::matrix<T, 3,rotgen::Dynamic, O::value>::Identity(3, 3), 3, 3);
+  test_identity(rotgen::matrix<T,rotgen::Dynamic, 11, O::value>::Identity(5, 11), 5, 11);
+};
\ No newline at end of file
diff --git a/test/unit/tests.hpp b/test/unit/tests.hpp
new file mode 100644
index 0000000..7798c6d
--- /dev/null
+++ b/test/unit/tests.hpp
@@ -0,0 +1,31 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+#define TTS_MAIN
+#include "tts.hpp"
+#include <rotgen/detail/static_info.hpp>
+#include <functional>
+
+namespace rotgen::tests
+{
+  template<auto N> struct constant { static constexpr auto value = N; };
+
+  using scalar  = tts::types<float,double>;
+  using order   = tts::types<constant<ColMajor>,constant<RowMajor>>;
+
+  using types   = tts::types< tts::types<float ,constant<ColMajor>>
+                            , tts::types<double,constant<ColMajor>>
+                            , tts::types<float ,constant<RowMajor>>
+                            , tts::types<double,constant<RowMajor>>
+                            >;
+
+  struct matrix_descriptor
+  {
+    std::size_t                                     rows, cols;
+    std::function<double(std::size_t,std::size_t)>  init_fn;
+  };
+}
\ No newline at end of file