diff --git a/include/rotgen/impl/matrix_impl64.hpp b/include/rotgen/impl/matrix_impl64.hpp
index e3345a7..1651f71 100644
--- a/include/rotgen/impl/matrix_impl64.hpp
+++ b/include/rotgen/impl/matrix_impl64.hpp
@@ -36,6 +36,22 @@ namespace rotgen
     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&       operator()(std::size_t i, std::size_t j);
     double const& operator()(std::size_t i, std::size_t j) const;
 
@@ -43,8 +59,11 @@ namespace rotgen
     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);
diff --git a/include/rotgen/matrix.hpp b/include/rotgen/matrix.hpp
index a481c23..febe15b 100644
--- a/include/rotgen/matrix.hpp
+++ b/include/rotgen/matrix.hpp
@@ -29,6 +29,8 @@ namespace rotgen
        if constexpr(Cols != -1) assert(c == Cols &&  "Mismatched between dynamic and static column size");
     }
 
+    matrix(parent const& base) : parent(base) {}
+
     matrix(std::initializer_list<std::initializer_list<Scalar>> init) : parent(init)
     {
       if constexpr(Rows != -1) assert(init.size() == Rows &&  "Mismatched between dynamic and static row size");
@@ -58,6 +60,25 @@ namespace rotgen
       parent::conservativeResize(new_rows, new_cols);
     }
 
+    matrix transpose() const
+    {
+      return matrix(static_cast<parent const&>(*this).transpose());
+    }
+
+    matrix conjugate() const
+    {
+      return matrix(static_cast<parent const&>(*this).conjugate());
+    }
+
+    matrix adjoint() const
+    {
+      return matrix(static_cast<parent const&>(*this).adjoint());
+    }
+
+    void transposeInPlace() { parent::transposeInPlace(); }
+
+    void adjointInPlace() { parent::adjointInPlace(); }
+
     friend bool operator==(matrix const& lhs, matrix const& rhs)
     {
       return static_cast<parent const&>(lhs) == static_cast<parent const&>(rhs);
@@ -69,6 +90,17 @@ namespace rotgen
       return *this;
     }
 
+    matrix& operator-=(matrix const& rhs)
+    {
+      static_cast<parent&>(*this) -= static_cast<parent const&>(rhs);
+      return *this;
+    }
+
+    matrix operator-() const
+    {
+      return matrix(static_cast<parent const&>(*this).operator-());
+    }
+
     matrix& operator*=(matrix const& rhs)
     {
       static_cast<parent&>(*this) *= static_cast<parent const&>(rhs);
@@ -80,6 +112,12 @@ namespace rotgen
       static_cast<parent&>(*this) *= rhs;
       return *this;
     }
+
+    matrix& operator/=(double rhs)
+    {
+      static_cast<parent&>(*this) /= rhs;
+      return *this;
+    }
   };
 
   template<typename S, int R, int C, int O, int MR, int MC>
@@ -89,6 +127,13 @@ namespace rotgen
     return that += rhs;
   }
 
+  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, matrix<S,R,C,O,MR,MC> const& 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*(matrix<S,R,C,O,MR,MC> const& lhs, matrix<S,R,C,O,MR,MC> const& rhs)
   {
@@ -108,4 +153,11 @@ namespace rotgen
   {
     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> that(lhs);
+    return that /= rhs;
+  }
 }
\ No newline at end of file
diff --git a/src/matrix_impl64.cpp b/src/matrix_impl64.cpp
index e1fec52..f70259a 100644
--- a/src/matrix_impl64.cpp
+++ b/src/matrix_impl64.cpp
@@ -76,6 +76,50 @@ namespace rotgen
 
   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(); }
+
   //==================================================================================================
   // Operators
   //==================================================================================================
@@ -95,6 +139,19 @@ namespace rotgen
     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;
@@ -106,4 +163,10 @@ namespace rotgen
     storage_->data *= s;
     return *this;
   }
+
+  matrix_impl64& matrix_impl64::operator/=(double s)
+  {
+    storage_->data /= s;
+    return *this;
+  }
 }
\ No newline at end of file
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index e2dfd7b..b9dddab 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -12,7 +12,7 @@ rotgen_setup_test_targets()
 add_library(rotgen_test INTERFACE)
 target_compile_options(rotgen_test INTERFACE -std=c++20 -Werror -Wall -Wextra -Wshadow -Wunused-variable)
 target_include_directories( rotgen_test INTERFACE ${PROJECT_SOURCE_DIR}/test)
-target_link_libraries(rotgen_test INTERFACE rotgen)
+target_link_libraries(rotgen_test INTERFACE rotgen Eigen3::Eigen)
 
 ##======================================================================================================================
 ## Unit Tests registration
diff --git a/test/basic/arithmetic_functions.cpp b/test/basic/arithmetic_functions.cpp
new file mode 100644
index 0000000..c61f505
--- /dev/null
+++ b/test/basic/arithmetic_functions.cpp
@@ -0,0 +1,121 @@
+//==================================================================================================
+/*
+  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/operators.cpp b/test/basic/operators.cpp
index 0734aef..a3051ad 100644
--- a/test/basic/operators.cpp
+++ b/test/basic/operators.cpp
@@ -9,23 +9,257 @@
 #include <rotgen/matrix.hpp>
 #include "tts.hpp"
 
-TTS_CASE("Check operator*")
+template <typename MatrixType>
+void test_matrix_scalar_multiplication(std::size_t rows, std::size_t cols, double scalar, auto init_fn)
 {
-  rotgen::matrix<double> a(2,2);
-  rotgen::matrix<double> ref(2,2);
+  MatrixType a(rows, cols);
+  MatrixType ref(rows, cols);
 
-  for(std::size_t r=0;r<a.rows();r++)
+  for (std::size_t r = 0; r < rows; ++r)
   {
-    for(std::size_t c=0;c<a.cols();c++)
+    for (std::size_t c = 0; c < cols; ++c)
     {
-      a(r,c)   = (1+c) + 10*(1+r);
-      ref(r,c) = ((1+c) + 10*(1+r)) * 10.5;
+      init_fn(a, r, c);
+      ref(r, c) = a(r, c) * scalar;
     }
   }
 
-  TTS_EQUAL(a * 10.5, ref);
-  TTS_EQUAL(10.5 * a, ref);
+  TTS_EQUAL(a * scalar, ref);
+  TTS_EQUAL(scalar * a, ref);
 
-  a *= 10.5;
+  a *= scalar;
   TTS_EQUAL(a, ref);
-};
\ No newline at end of file
+}
+
+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; });
+};