diff --git a/src/algebra/SmallMatrix.hpp b/src/algebra/SmallMatrix.hpp
index d0d6ee329f02c41626900384c794d1cfa490089c..b8d3b81b5fad0ac65521a94aa4f94183f4cff333 100644
--- a/src/algebra/SmallMatrix.hpp
+++ b/src/algebra/SmallMatrix.hpp
@@ -38,7 +38,7 @@ class [[nodiscard]] SmallMatrix   // LCOV_EXCL_LINE
 
   friend PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> copy(const SmallMatrix& A) noexcept
   {
-    return {A.m_nb_rows, A.m_nb_columns, copy(A.m_values)};
+    return SmallMatrix<std::remove_const_t<DataType>>{A.m_nb_rows, A.m_nb_columns, copy(A.m_values)};
   }
 
   friend PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> transpose(const SmallMatrix& A)
@@ -130,7 +130,8 @@ class [[nodiscard]] SmallMatrix   // LCOV_EXCL_LINE
   {
     static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
                   "incompatible data types");
-    Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()), "incompatible matrix sizes");
+    Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()),
+           "cannot add matrix: incompatible sizes");
 
     parallel_for(
       m_values.size(), PUGS_LAMBDA(index_type i) { m_values[i] += B.m_values[i]; });
@@ -260,6 +261,12 @@ class [[nodiscard]] SmallMatrix   // LCOV_EXCL_LINE
 
   SmallMatrix(SmallMatrix &&) = default;
 
+  explicit SmallMatrix(size_t nb_rows, size_t nb_columns, const SmallArray<DataType>& values)
+    : m_nb_rows{nb_rows}, m_nb_columns{nb_columns}, m_values{values}
+  {
+    Assert(m_values.size() == m_nb_columns * m_nb_rows, "incompatible array size and matrix dimensions")
+  }
+
   explicit SmallMatrix(size_t nb_rows, size_t nb_columns) noexcept
     : m_nb_rows{nb_rows}, m_nb_columns{nb_columns}, m_values{nb_rows * nb_columns}
   {}
diff --git a/src/algebra/SmallVector.hpp b/src/algebra/SmallVector.hpp
index d20bb48f4cdb700f61a0881898fb5273217693aa..d4a0bf84268bdf1d164cf6216153117adafaab56 100644
--- a/src/algebra/SmallVector.hpp
+++ b/src/algebra/SmallVector.hpp
@@ -30,6 +30,18 @@ class SmallVector   // LCOV_EXCL_LINE
     return os << x.m_values;
   }
 
+  [[nodiscard]] friend std::remove_const_t<DataType>
+  min(const SmallVector& x)
+  {
+    return min(x.m_values);
+  }
+
+  [[nodiscard]] friend std::remove_const_t<DataType>
+  max(const SmallVector& x)
+  {
+    return max(x.m_values);
+  }
+
   friend SmallVector<std::remove_const_t<DataType>>
   copy(const SmallVector& x)
   {
@@ -40,7 +52,18 @@ class SmallVector   // LCOV_EXCL_LINE
     return x_copy;
   }
 
-  friend SmallVector operator*(const DataType& a, const SmallVector& x)
+  PUGS_INLINE
+  SmallVector<std::remove_const_t<DataType>>
+  operator-() const
+  {
+    SmallVector<std::remove_const_t<DataType>> opposite(this->size());
+    parallel_for(
+      opposite.size(), PUGS_LAMBDA(index_type i) { opposite.m_values[i] = -m_values[i]; });
+    return opposite;
+  }
+
+  friend SmallVector
+  operator*(const DataType& a, const SmallVector& x)
   {
     SmallVector<std::remove_const_t<DataType>> y = copy(x);
     return y *= a;
@@ -51,14 +74,8 @@ class SmallVector   // LCOV_EXCL_LINE
   dot(const SmallVector& x, const SmallVector<DataType2>& y)
   {
     Assert(x.size() == y.size(), "cannot compute dot product: incompatible vector sizes");
-    // Quite ugly, TODO: fix me in C++20
-    auto promoted = [] {
-      DataType a{0};
-      DataType2 b{0};
-      return a * b;
-    }();
 
-    decltype(promoted) sum = 0;
+    decltype(std::declval<DataType>() * std::declval<DataType2>()) sum = 0;
 
     // Does not use parallel_for to preserve sum order
     for (index_type i = 0; i < x.size(); ++i) {
@@ -136,7 +153,8 @@ class SmallVector   // LCOV_EXCL_LINE
   }
 
   PUGS_INLINE
-  DataType& operator[](index_type i) const noexcept(NO_ASSERT)
+  DataType&
+  operator[](index_type i) const noexcept(NO_ASSERT)
   {
     return m_values[i];
   }
diff --git a/src/algebra/TinyVector.hpp b/src/algebra/TinyVector.hpp
index ae616a4849d0e681bf09e613ad4d19263d0f035a..e39dcec6e9834e5cb5eb80e1e845de1054beb1b2 100644
--- a/src/algebra/TinyVector.hpp
+++ b/src/algebra/TinyVector.hpp
@@ -42,14 +42,12 @@ class [[nodiscard]] TinyVector
     return opposite;
   }
 
-  PUGS_INLINE
-  constexpr size_t dimension() const
+  [[nodiscard]] PUGS_INLINE constexpr size_t dimension() const
   {
     return N;
   }
 
-  PUGS_INLINE
-  constexpr bool operator==(const TinyVector& v) const
+  [[nodiscard]] PUGS_INLINE constexpr bool operator==(const TinyVector& v) const
   {
     for (size_t i = 0; i < N; ++i) {
       if (m_values[i] != v.m_values[i])
@@ -58,14 +56,12 @@ class [[nodiscard]] TinyVector
     return true;
   }
 
-  PUGS_INLINE
-  constexpr bool operator!=(const TinyVector& v) const
+  [[nodiscard]] PUGS_INLINE constexpr bool operator!=(const TinyVector& v) const
   {
     return not this->operator==(v);
   }
 
-  PUGS_INLINE
-  constexpr friend T dot(const TinyVector& u, const TinyVector& v)
+  [[nodiscard]] PUGS_INLINE constexpr friend T dot(const TinyVector& u, const TinyVector& v)
   {
     T t = u.m_values[0] * v.m_values[0];
     for (size_t i = 1; i < N; ++i) {
@@ -242,7 +238,7 @@ class [[nodiscard]] TinyVector
 };
 
 template <size_t N, typename T>
-PUGS_INLINE constexpr T
+[[nodiscard]] PUGS_INLINE constexpr T
 l2Norm(const TinyVector<N, T>& x)
 {
   static_assert(std::is_arithmetic<T>(), "Cannot compute L2 norm for non-arithmetic types");
@@ -250,13 +246,39 @@ l2Norm(const TinyVector<N, T>& x)
   return std::sqrt(dot(x, x));
 }
 
-// Cross product is only defined for K^3 vectors
+template <size_t N, typename T>
+[[nodiscard]] PUGS_INLINE constexpr T
+min(const TinyVector<N, T>& x)
+{
+  T m = x[0];
+  for (size_t i = 1; i < N; ++i) {
+    if (x[i] < m) {
+      m = x[i];
+    }
+  }
+  return m;
+}
+
+template <size_t N, typename T>
+[[nodiscard]] PUGS_INLINE constexpr T
+max(const TinyVector<N, T>& x)
+{
+  T m = x[0];
+  for (size_t i = 1; i < N; ++i) {
+    if (x[i] > m) {
+      m = x[i];
+    }
+  }
+  return m;
+}
+
+// Cross product is only defined for dimension 3 vectors
 template <typename T>
-PUGS_INLINE constexpr TinyVector<3, T>
+[[nodiscard]] PUGS_INLINE constexpr TinyVector<3, T>
 crossProduct(const TinyVector<3, T>& u, const TinyVector<3, T>& v)
 {
   TinyVector<3, T> cross_product(u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]);
   return cross_product;
 }
 
-#endif   // TINYVECTOR_HPP
+#endif   // TINY_VECTOR_HPP
diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp
index c918ed409746a3c212720ad8ae6769cfb0882ba9..749a795a489f12a671f4ecfd452104702ba68147 100644
--- a/src/algebra/Vector.hpp
+++ b/src/algebra/Vector.hpp
@@ -31,6 +31,18 @@ class Vector   // LCOV_EXCL_LINE
     return os << x.m_values;
   }
 
+  [[nodiscard]] friend std::remove_const_t<DataType>
+  min(const Vector& x)
+  {
+    return min(x.m_values);
+  }
+
+  [[nodiscard]] friend std::remove_const_t<DataType>
+  max(const Vector& x)
+  {
+    return max(x.m_values);
+  }
+
   friend Vector<std::remove_const_t<DataType>>
   copy(const Vector& x)
   {
@@ -41,6 +53,16 @@ class Vector   // LCOV_EXCL_LINE
     return x_copy;
   }
 
+  PUGS_INLINE
+  Vector<std::remove_const_t<DataType>>
+  operator-() const
+  {
+    Vector<std::remove_const_t<DataType>> opposite(this->size());
+    parallel_for(
+      opposite.size(), PUGS_LAMBDA(index_type i) { opposite.m_values[i] = -m_values[i]; });
+    return opposite;
+  }
+
   friend Vector
   operator*(const DataType& a, const Vector& x)
   {
@@ -53,14 +75,7 @@ class Vector   // LCOV_EXCL_LINE
   dot(const Vector& x, const Vector<DataType2>& y)
   {
     Assert(x.size() == y.size(), "cannot compute dot product: incompatible vector sizes");
-    // Quite ugly, TODO: fix me in C++20
-    auto promoted = [] {
-      DataType a{0};
-      DataType2 b{0};
-      return a * b;
-    }();
-
-    decltype(promoted) sum = 0;
+    decltype(std::declval<DataType>() * std::declval<DataType2>()) sum = 0;
 
     // Does not use parallel_for to preserve sum order
     for (index_type i = 0; i < x.size(); ++i) {
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 0e585acac5ba1e5b7656ffa29411567e0f0d54f4..b5301f36ae05dd0ad9e034a4953f7d7fa95d17ec 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -131,6 +131,7 @@ add_executable (unit_tests
   test_RefItemList.cpp
   test_RevisionInfo.cpp
   test_SmallArray.cpp
+  test_SmallMatrix.cpp
   test_SmallVector.cpp
   test_Socket.cpp
   test_SocketModule.cpp
diff --git a/tests/test_SmallMatrix.cpp b/tests/test_SmallMatrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9712e13e330ee84adc55afeee1bcf81150312bc9
--- /dev/null
+++ b/tests/test_SmallMatrix.cpp
@@ -0,0 +1,416 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <Kokkos_Core.hpp>
+
+#include <utils/PugsAssert.hpp>
+#include <utils/Types.hpp>
+
+#include <algebra/SmallMatrix.hpp>
+
+#include <sstream>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SmallMatrix", "[algebra]")
+{
+  SmallMatrix<int> A{3, 4};
+
+  REQUIRE(A.numberOfRows() == 3);
+  REQUIRE(A.numberOfColumns() == 4);
+
+  for (size_t i = 0; i < 3; ++i) {
+    for (size_t j = 0; j < 4; ++j) {
+      A(i, j) = 1 + 4 * i + j;
+    }
+  }
+
+  REQUIRE(not A.isSquare());
+
+  REQUIRE(((A(0, 0) == 1) and (A(0, 1) == 2) and (A(0, 2) == 3) and (A(0, 3) == 4) and   //
+           (A(1, 0) == 5) and (A(1, 1) == 6) and (A(1, 2) == 7) and (A(1, 3) == 8) and   //
+           (A(2, 0) == 9) and (A(2, 1) == 10) and (A(2, 2) == 11) and (A(2, 3) == 12)));
+
+  SmallMatrix<int> B{3, 4};
+  B(0, 0) = -6;
+  B(0, 1) = 5;
+  B(0, 2) = 3;
+  B(0, 3) = 8;
+
+  B(1, 0) = -4;
+  B(1, 1) = 6;
+  B(1, 2) = 5;
+  B(1, 3) = 3;
+
+  B(2, 0) = 7;
+  B(2, 1) = 1;
+  B(2, 2) = 3;
+  B(2, 3) = 6;
+
+  SECTION("copy and shallow copy")
+  {
+    SmallMatrix<const int> affected_A;
+
+    REQUIRE(affected_A.numberOfRows() == 0);
+    REQUIRE(affected_A.numberOfColumns() == 0);
+
+    affected_A = A;
+    REQUIRE(&A(0, 0) == &affected_A(0, 0));
+
+    REQUIRE(affected_A.numberOfRows() == A.numberOfRows());
+    REQUIRE(affected_A.numberOfColumns() == A.numberOfColumns());
+
+    const SmallMatrix<int> copy_A = copy(A);
+
+    REQUIRE(copy_A.numberOfRows() == A.numberOfRows());
+    REQUIRE(copy_A.numberOfColumns() == A.numberOfColumns());
+
+    bool is_same = true;
+    for (size_t i = 0; i < A.numberOfRows(); ++i) {
+      for (size_t j = 0; j < A.numberOfColumns(); ++j) {
+        is_same &= (copy_A(i, j) == A(i, j));
+      }
+    }
+    REQUIRE(is_same);
+
+    REQUIRE(&A(0, 0) != &copy_A(0, 0));
+  }
+
+  SECTION("fill and special affectations")
+  {
+    SmallMatrix<int> M(5);
+    REQUIRE(M.isSquare());
+    M.fill(1);
+    for (size_t i = 0; i < M.numberOfRows(); ++i) {
+      for (size_t j = 0; j < M.numberOfColumns(); ++j) {
+        REQUIRE(M(i, j) == 1);
+      }
+    }
+
+    M.fill(-4);
+    for (size_t i = 0; i < M.numberOfRows(); ++i) {
+      for (size_t j = 0; j < M.numberOfColumns(); ++j) {
+        REQUIRE(M(i, j) == -4);
+      }
+    }
+
+    M = zero;
+    for (size_t i = 0; i < M.numberOfRows(); ++i) {
+      for (size_t j = 0; j < M.numberOfColumns(); ++j) {
+        REQUIRE(M(i, j) == 0);
+      }
+    }
+
+    M = identity;
+    for (size_t i = 0; i < M.numberOfRows(); ++i) {
+      for (size_t j = 0; j < M.numberOfColumns(); ++j) {
+        REQUIRE(M(i, j) == (i == j));
+      }
+    }
+  }
+
+  SECTION("build from TinyMatrix")
+  {
+    TinyMatrix<3, 4> T{1, 2, 3, 4, 0.5, 1, 1.5, 2, 0.3, 0.6, 0.9, 1.2};
+    SmallMatrix S{T};
+
+    REQUIRE(S.numberOfRows() == T.numberOfRows());
+    REQUIRE(S.numberOfColumns() == T.numberOfColumns());
+
+    for (size_t i = 0; i < S.numberOfRows(); ++i) {
+      for (size_t j = 0; j < S.numberOfColumns(); ++j) {
+        REQUIRE(S(i, j) == T(i, j));
+      }
+    }
+  }
+
+  SmallMatrix AT = transpose(A);
+
+  SECTION("transposed")
+  {
+    REQUIRE(AT.numberOfColumns() == A.numberOfRows());
+    REQUIRE(AT.numberOfRows() == A.numberOfColumns());
+
+    bool is_same = true;
+    for (size_t i = 0; i < A.numberOfRows(); ++i) {
+      for (size_t j = 0; j < A.numberOfColumns(); ++j) {
+        is_same &= (AT(j, i) == A(i, j));
+      }
+    }
+    REQUIRE(is_same);
+  }
+
+  SECTION("sum")
+  {
+    SmallMatrix AaddB = copy(A);
+    AaddB += B;
+
+    REQUIRE(AaddB(0, 0) == -5);
+    REQUIRE(AaddB(0, 1) == 7);
+    REQUIRE(AaddB(0, 2) == 6);
+    REQUIRE(AaddB(0, 3) == 12);
+
+    REQUIRE(AaddB(1, 0) == 1);
+    REQUIRE(AaddB(1, 1) == 12);
+    REQUIRE(AaddB(1, 2) == 12);
+    REQUIRE(AaddB(1, 3) == 11);
+
+    REQUIRE(AaddB(2, 0) == 16);
+    REQUIRE(AaddB(2, 1) == 11);
+    REQUIRE(AaddB(2, 2) == 14);
+    REQUIRE(AaddB(2, 3) == 18);
+
+    SmallMatrix ApB = A + B;
+    for (size_t i = 0; i < ApB.numberOfRows(); ++i) {
+      for (size_t j = 0; j < ApB.numberOfColumns(); ++j) {
+        REQUIRE(AaddB(i, j) == ApB(i, j));
+      }
+    }
+  }
+
+  SECTION("difference")
+  {
+    SmallMatrix AsubsB = copy(A);
+    AsubsB -= B;
+
+    REQUIRE(AsubsB(0, 0) == 7);
+    REQUIRE(AsubsB(0, 1) == -3);
+    REQUIRE(AsubsB(0, 2) == 0);
+    REQUIRE(AsubsB(0, 3) == -4);
+
+    REQUIRE(AsubsB(1, 0) == 9);
+    REQUIRE(AsubsB(1, 1) == 0);
+    REQUIRE(AsubsB(1, 2) == 2);
+    REQUIRE(AsubsB(1, 3) == 5);
+
+    REQUIRE(AsubsB(2, 0) == 2);
+    REQUIRE(AsubsB(2, 1) == 9);
+    REQUIRE(AsubsB(2, 2) == 8);
+    REQUIRE(AsubsB(2, 3) == 6);
+
+    SmallMatrix AmB = A - B;
+    for (size_t i = 0; i < AmB.numberOfRows(); ++i) {
+      for (size_t j = 0; j < AmB.numberOfColumns(); ++j) {
+        REQUIRE(AsubsB(i, j) == AmB(i, j));
+      }
+    }
+  }
+
+  SECTION("left product by a scalar")
+  {
+    int a = 2;
+
+    SmallMatrix aA = a * A;
+
+    {
+      bool is_same = true;
+      for (size_t i = 0; i < A.numberOfRows(); ++i) {
+        for (size_t j = 0; j < A.numberOfColumns(); ++j) {
+          is_same &= (aA(i, j) == a * A(i, j));
+        }
+      }
+      REQUIRE(is_same);
+    }
+
+    SmallMatrix copy_A = copy(A);
+    copy_A *= a;
+
+    {
+      bool is_same = true;
+      for (size_t i = 0; i < A.numberOfRows(); ++i) {
+        for (size_t j = 0; j < A.numberOfColumns(); ++j) {
+          is_same &= (aA(i, j) == copy_A(i, j));
+        }
+      }
+      REQUIRE(is_same);
+    }
+  }
+
+  SECTION("divide by a scalar")
+  {
+    SmallMatrix<double> M(3, 5);
+    M(0, 0) = 2;
+    M(0, 1) = 4;
+    M(0, 2) = 6;
+    M(0, 3) = 8;
+    M(0, 4) = 10;
+
+    M(1, 0) = 12;
+    M(1, 1) = 14;
+    M(1, 2) = 16;
+    M(1, 3) = 18;
+    M(1, 4) = 20;
+
+    M(2, 0) = 22;
+    M(2, 1) = 24;
+    M(2, 2) = 26;
+    M(2, 3) = 28;
+    M(2, 4) = 30;
+
+    M /= 2;
+
+    for (size_t i = 0; i < M.numberOfRows(); ++i) {
+      for (size_t j = 0; j < M.numberOfColumns(); ++j) {
+        REQUIRE(M(i, j) == 1 + i * M.numberOfColumns() + j);
+      }
+    }
+  }
+
+  SECTION("matrix-vector product")
+  {
+    SmallVector<int> u(4);
+
+    u[0] = 1;
+    u[1] = 3;
+    u[2] = -2;
+    u[3] = -1;
+
+    SmallVector v = B * u;
+    REQUIRE(((v[0] == -5) and (v[1] == 1) and (v[2] == -2)));
+  }
+
+  SECTION("matrix-matrix product")
+  {
+    SmallMatrix<int> BAT = B * AT;
+    REQUIRE(BAT.isSquare());
+    REQUIRE(BAT.numberOfRows() == 3);
+    REQUIRE(BAT.numberOfColumns() == 3);
+
+    REQUIRE(BAT(0, 0) == 45);
+    REQUIRE(BAT(0, 1) == 85);
+    REQUIRE(BAT(0, 2) == 125);
+
+    REQUIRE(BAT(1, 0) == 35);
+    REQUIRE(BAT(1, 1) == 75);
+    REQUIRE(BAT(1, 2) == 115);
+
+    REQUIRE(BAT(2, 0) == 42);
+    REQUIRE(BAT(2, 1) == 110);
+    REQUIRE(BAT(2, 2) == 178);
+
+    SmallMatrix<int> ATB = AT * B;
+    REQUIRE(ATB.isSquare());
+    REQUIRE(ATB.numberOfRows() == 4);
+    REQUIRE(ATB.numberOfColumns() == 4);
+
+    REQUIRE(ATB(0, 0) == 37);
+    REQUIRE(ATB(0, 1) == 44);
+    REQUIRE(ATB(0, 2) == 55);
+    REQUIRE(ATB(0, 3) == 77);
+
+    REQUIRE(ATB(1, 0) == 34);
+    REQUIRE(ATB(1, 1) == 56);
+    REQUIRE(ATB(1, 2) == 66);
+    REQUIRE(ATB(1, 3) == 94);
+
+    REQUIRE(ATB(2, 0) == 31);
+    REQUIRE(ATB(2, 1) == 68);
+    REQUIRE(ATB(2, 2) == 77);
+    REQUIRE(ATB(2, 3) == 111);
+
+    REQUIRE(ATB(3, 0) == 28);
+    REQUIRE(ATB(3, 1) == 80);
+    REQUIRE(ATB(3, 2) == 88);
+    REQUIRE(ATB(3, 3) == 128);
+  }
+
+  SECTION("output")
+  {
+    std::ostringstream out;
+    out << A;
+
+    std::ostringstream expected;
+    expected << "0| 0:1 1:2 2:3 3:4\n";
+    expected << "1| 0:5 1:6 2:7 3:8\n";
+    expected << "2| 0:9 1:10 2:11 3:12\n";
+
+    REQUIRE(out.str() == expected.str());
+  }
+
+#ifndef NDEBUG
+  SECTION("matrix-vector compatibility")
+  {
+    SmallVector<int> v(5);
+    REQUIRE_THROWS_WITH(A * v, "cannot compute matrix-vector product: incompatible sizes");
+  }
+
+  SECTION("matrix-matrix compatibility")
+  {
+    SmallMatrix C = copy(A);
+    REQUIRE_THROWS_WITH(C * A, "cannot compute matrix product: incompatible sizes");
+    REQUIRE_THROWS_WITH(C -= AT, "cannot substract matrix: incompatible sizes");
+    REQUIRE_THROWS_WITH(C += AT, "cannot add matrix: incompatible sizes");
+    REQUIRE_THROWS_WITH(C + AT, "cannot compute matrix sum: incompatible sizes");
+    REQUIRE_THROWS_WITH(C - AT, "cannot compute matrix difference: incompatible sizes");
+  }
+
+  SECTION("invalid element")
+  {
+    REQUIRE_THROWS_WITH(A(A.numberOfRows(), 0), "cannot access element: invalid indices");
+    REQUIRE_THROWS_WITH(A(0, A.numberOfColumns()), "cannot access element: invalid indices");
+  }
+
+  SECTION("invalid identity affectation")
+  {
+    SmallMatrix<double> C(3, 2);
+    REQUIRE_THROWS_WITH((C = identity), "identity must be a square matrix");
+  }
+
+  SECTION("invalid constructor")
+  {
+    REQUIRE_THROWS_WITH((SmallMatrix{2, 3, SmallArray<double>{7}}), "incompatible array size and matrix dimensions");
+  }
+
+  SECTION("output with signaling NaN")
+  {
+    SmallMatrix<double> A(2, 3);
+    A(0, 0) = 1;
+    A(0, 2) = 3;
+    A(1, 0) = 2;
+    std::ostringstream A_ost;
+    A_ost << A;
+
+    std::ostringstream ref_ost;
+    ref_ost << "0| 0:1 1:nan 2:3\n";
+    ref_ost << "1| 0:2 1:nan 2:nan\n";
+
+    REQUIRE(A_ost.str() == ref_ost.str());
+  }
+
+  // SECTION("checking for bounds violation")
+  // {
+  //   REQUIRE_THROWS_AS(A(3, 0), AssertError);
+  //   REQUIRE_THROWS_AS(A(0, 4), AssertError);
+
+  //   REQUIRE_THROWS_AS(getMinor(A, 3, 0), AssertError);
+  //   REQUIRE_THROWS_AS(getMinor(A, 0, 4), AssertError);
+
+  //   const TinyMatrix<3, 4, int>& constA = A;
+  //   REQUIRE_THROWS_AS(constA(3, 0), AssertError);
+  //   REQUIRE_THROWS_AS(constA(0, 4), AssertError);
+  // }
+
+  // SECTION("checking for nan initialization")
+  // {
+  //   TinyMatrix<3, 4, double> B;
+
+  //   for (size_t i = 0; i < B.numberOfRows(); ++i) {
+  //     for (size_t j = 0; j < B.numberOfColumns(); ++j) {
+  //       REQUIRE(std::isnan(B(i, j)));
+  //     }
+  //   }
+  // }
+
+  // SECTION("checking for bad initialization")
+  // {
+  //   TinyMatrix<3, 4, int> B;
+
+  //   for (size_t i = 0; i < B.numberOfRows(); ++i) {
+  //     for (size_t j = 0; j < B.numberOfColumns(); ++j) {
+  //       REQUIRE(B(i, j) == std::numeric_limits<int>::max() / 2);
+  //     }
+  //   }
+  // }
+#endif   // NDEBUG
+}
diff --git a/tests/test_SmallVector.cpp b/tests/test_SmallVector.cpp
index 3b499c3ab3e984247b7f20f96b626f8ee413763a..b5ffc6819736b2e55d33706aeb18f1e1d28dabca 100644
--- a/tests/test_SmallVector.cpp
+++ b/tests/test_SmallVector.cpp
@@ -360,6 +360,60 @@ TEST_CASE("SmallVector", "[algebra]")
     REQUIRE(z[4] == 4);
   }
 
+  SECTION("unary minus")
+  {
+    SmallVector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    SmallVector<const int> y = x;
+
+    SmallVector<int> z = -x;
+
+    REQUIRE(z[0] == -2);
+    REQUIRE(z[1] == -3);
+    REQUIRE(z[2] == -5);
+    REQUIRE(z[3] == -7);
+    REQUIRE(z[4] == -8);
+
+    z = -y;
+    REQUIRE(z[0] == -2);
+    REQUIRE(z[1] == -3);
+    REQUIRE(z[2] == -5);
+    REQUIRE(z[3] == -7);
+    REQUIRE(z[4] == -8);
+  }
+
+  SECTION("min/max")
+  {
+    SmallVector<int> u(5);
+    u[0] = 1;
+    u[1] = 7;
+    u[2] = 6;
+    u[3] = 2;
+    u[4] = 9;
+
+    REQUIRE(min(u) == 1);
+    REQUIRE(max(u) == 9);
+    REQUIRE(min(-u) == -9);
+    REQUIRE(max(-u) == -1);
+
+    SmallVector<int> v(5);
+    v[0] = 1;
+    v[1] = 11;
+    v[2] = 6;
+    v[3] = -2;
+    v[4] = 9;
+
+    REQUIRE(min(v) == -2);
+    REQUIRE(max(v) == 11);
+    REQUIRE(min(-v) == -11);
+    REQUIRE(max(-v) == 2);
+  }
+
   SECTION("output")
   {
     SmallVector<int> x{5};
diff --git a/tests/test_TinyVector.cpp b/tests/test_TinyVector.cpp
index 9c4dc2e5a2c4b185904fcf04a706fa70b572a197..6a3942e84b1c95680d8c63080b5c8b37330f139b 100644
--- a/tests/test_TinyVector.cpp
+++ b/tests/test_TinyVector.cpp
@@ -87,6 +87,23 @@ TEST_CASE("TinyVector", "[algebra]")
     REQUIRE(crossProduct(a, b) == TinyVector<3, int>(-16, 6, 7));
   }
 
+  SECTION("min/max")
+  {
+    TinyVector<5, int> u{1, 7, 6, 2, 9};
+
+    REQUIRE(min(u) == 1);
+    REQUIRE(max(u) == 9);
+    REQUIRE(min(-u) == -9);
+    REQUIRE(max(-u) == -1);
+
+    TinyVector<5, int> v{1, 11, 6, -2, 9};
+
+    REQUIRE(min(v) == -2);
+    REQUIRE(max(v) == 11);
+    REQUIRE(min(-v) == -11);
+    REQUIRE(max(-v) == 2);
+  }
+
 #ifndef NDEBUG
   SECTION("output with signaling NaN")
   {
diff --git a/tests/test_Vector.cpp b/tests/test_Vector.cpp
index e31bf166713b70f69ca108d7dace6dc07179dc5b..0d8d0e898321070f74b3cf99f8fb0e04686c9b0d 100644
--- a/tests/test_Vector.cpp
+++ b/tests/test_Vector.cpp
@@ -359,6 +359,60 @@ TEST_CASE("Vector", "[algebra]")
     REQUIRE(z[4] == 4);
   }
 
+  SECTION("unary minus")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    Vector<const int> y = x;
+
+    Vector<int> z = -x;
+
+    REQUIRE(z[0] == -2);
+    REQUIRE(z[1] == -3);
+    REQUIRE(z[2] == -5);
+    REQUIRE(z[3] == -7);
+    REQUIRE(z[4] == -8);
+
+    z = -y;
+    REQUIRE(z[0] == -2);
+    REQUIRE(z[1] == -3);
+    REQUIRE(z[2] == -5);
+    REQUIRE(z[3] == -7);
+    REQUIRE(z[4] == -8);
+  }
+
+  SECTION("min/max")
+  {
+    Vector<int> u(5);
+    u[0] = 1;
+    u[1] = 7;
+    u[2] = 6;
+    u[3] = 2;
+    u[4] = 9;
+
+    REQUIRE(min(u) == 1);
+    REQUIRE(max(u) == 9);
+    REQUIRE(min(-u) == -9);
+    REQUIRE(max(-u) == -1);
+
+    Vector<int> v(5);
+    v[0] = 1;
+    v[1] = 11;
+    v[2] = 6;
+    v[3] = -2;
+    v[4] = 9;
+
+    REQUIRE(min(v) == -2);
+    REQUIRE(max(v) == 11);
+    REQUIRE(min(-v) == -11);
+    REQUIRE(max(-v) == 2);
+  }
+
   SECTION("output")
   {
     Vector<int> x{5};