diff --git a/src/algebra/DenseMatrix.hpp b/src/algebra/DenseMatrix.hpp
index 8f15f125861d0cc9cce0aef716e63959f736d36b..25f607e7bf0bba35d8b919655e63f6e9431724fd 100644
--- a/src/algebra/DenseMatrix.hpp
+++ b/src/algebra/DenseMatrix.hpp
@@ -66,7 +66,7 @@ class DenseMatrix   // 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_columns == x.size());
+    Assert(m_nb_columns == x.size(), "cannot compute matrix-vector product: incompatible sizes");
     const DenseMatrix& A = *this;
     Vector<std::remove_const_t<DataType>> Ax{m_nb_rows};
     for (size_t i = 0; i < m_nb_rows; ++i) {
@@ -85,7 +85,7 @@ class DenseMatrix   // 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_columns == B.numberOfRows());
+    Assert(m_nb_columns == B.numberOfRows(), "cannot compute matrix product: incompatible sizes");
     const DenseMatrix& A = *this;
     DenseMatrix<std::remove_const_t<DataType>> AB{m_nb_rows, B.numberOfColumns()};
 
@@ -124,8 +124,8 @@ class DenseMatrix   // 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());
-    Assert(m_nb_columns == B.numberOfColumns());
+    Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()),
+           "cannot substract matrix: incompatible sizes");
 
     parallel_for(
       m_values.size(), PUGS_LAMBDA(index_type i) { m_values[i] -= B.m_values[i]; });
@@ -138,8 +138,7 @@ class DenseMatrix   // 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());
-    Assert(m_nb_columns == B.numberOfColumns());
+    Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()), "incompatible matrix sizes");
 
     parallel_for(
       m_values.size(), PUGS_LAMBDA(index_type i) { m_values[i] += B.m_values[i]; });
@@ -152,8 +151,9 @@ class DenseMatrix   // 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());
-    Assert(m_nb_columns == B.numberOfColumns());
+    Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()),
+           "cannot compute matrix sum: incompatible sizes");
+
     DenseMatrix<std::remove_const_t<DataType>> sum{B.numberOfRows(), B.numberOfColumns()};
 
     parallel_for(
@@ -168,8 +168,9 @@ class DenseMatrix   // 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());
-    Assert(m_nb_columns == B.numberOfColumns());
+    Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()),
+           "cannot compute matrix difference: incompatible sizes");
+
     DenseMatrix<std::remove_const_t<DataType>> difference{B.numberOfRows(), B.numberOfColumns()};
 
     parallel_for(
@@ -182,7 +183,7 @@ class DenseMatrix   // LCOV_EXCL_LINE
   DataType&
   operator()(index_type i, index_type j) const noexcept(NO_ASSERT)
   {
-    Assert((i < m_nb_rows and j < m_nb_columns), "invalid indices");
+    Assert(i < m_nb_rows and j < m_nb_columns, "cannot access element: invalid indices");
     return m_values[i * m_nb_columns + j];
   }
 
@@ -214,7 +215,7 @@ class DenseMatrix   // LCOV_EXCL_LINE
 
   PUGS_INLINE DenseMatrix& operator=(IdentityType) noexcept(NO_ASSERT)
   {
-    Assert(m_nb_rows == m_nb_columns, "Identity must be a square matrix");
+    Assert(m_nb_rows == m_nb_columns, "identity must be a square matrix");
 
     m_values.fill(0);
     parallel_for(
@@ -286,7 +287,7 @@ class DenseMatrix   // LCOV_EXCL_LINE
   DenseMatrix(size_t nb_rows, size_t nb_columns, const Array<DataType> values) noexcept(NO_ASSERT)
     : m_nb_rows{nb_rows}, m_nb_columns{nb_columns}, m_values{values}
   {
-    Assert(m_values.size() == m_nb_rows * m_nb_columns, "incompatible sizes");
+    Assert(m_values.size() == m_nb_rows * m_nb_columns, "cannot create matrix: incorrect number of values");
   }
 
  public:
diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp
index e9eace62afa018c31c9f727e779a01c9ce39ccd6..2a3e1bd296a4b07e0dd32e9f3017892a2cea290f 100644
--- a/src/algebra/Vector.hpp
+++ b/src/algebra/Vector.hpp
@@ -1,12 +1,11 @@
 #ifndef VECTOR_HPP
 #define VECTOR_HPP
 
+#include <utils/Array.hpp>
+#include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 #include <utils/PugsUtils.hpp>
-
-#include <utils/PugsAssert.hpp>
-
-#include <utils/Array.hpp>
+#include <utils/Types.hpp>
 
 template <typename DataType>
 class Vector   // LCOV_EXCL_LINE
@@ -20,6 +19,7 @@ class Vector   // LCOV_EXCL_LINE
   const bool m_deep_copies;
 
   static_assert(std::is_same_v<typename decltype(m_values)::index_type, index_type>);
+  static_assert(std::is_arithmetic_v<DataType>, "Vector is only defined for arithmetic values");
 
   // Allows const version to access our data
   friend Vector<std::add_const_t<DataType>>;
@@ -55,7 +55,7 @@ class Vector   // LCOV_EXCL_LINE
   PUGS_INLINE friend auto
   dot(const Vector& x, const Vector<DataType2>& y)
   {
-    Assert(x.size() == y.size());
+    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};
@@ -94,7 +94,7 @@ class Vector   // LCOV_EXCL_LINE
   PUGS_INLINE Vector&
   operator-=(const Vector<DataType2>& y)
   {
-    Assert(this->size() == y.size());
+    Assert(this->size() == y.size(), "cannot substract vector: incompatible sizes");
 
     parallel_for(
       this->size(), PUGS_LAMBDA(index_type i) { m_values[i] -= y[i]; });
@@ -106,7 +106,7 @@ class Vector   // LCOV_EXCL_LINE
   PUGS_INLINE Vector&
   operator+=(const Vector<DataType2>& y)
   {
-    Assert(this->size() == y.size());
+    Assert(this->size() == y.size(), "cannot add vector: incompatible sizes");
 
     parallel_for(
       this->size(), PUGS_LAMBDA(index_type i) { m_values[i] += y[i]; });
@@ -118,7 +118,7 @@ class Vector   // LCOV_EXCL_LINE
   PUGS_INLINE Vector<std::remove_const_t<DataType>>
   operator+(const Vector<DataType2>& y) const
   {
-    Assert(this->size() == y.size());
+    Assert(this->size() == y.size(), "cannot compute vector sum: incompatible sizes");
     Vector<std::remove_const_t<DataType>> sum{y.size()};
 
     parallel_for(
@@ -131,7 +131,7 @@ class Vector   // LCOV_EXCL_LINE
   PUGS_INLINE Vector<std::remove_const_t<DataType>>
   operator-(const Vector<DataType2>& y) const
   {
-    Assert(this->size() == y.size());
+    Assert(this->size() == y.size(), "cannot compute vector difference: incompatible sizes");
     Vector<std::remove_const_t<DataType>> sum{y.size()};
 
     parallel_for(
@@ -155,28 +155,34 @@ class Vector   // LCOV_EXCL_LINE
   }
 
   PUGS_INLINE Vector&
-  operator=(const DataType& value) noexcept
+  fill(const DataType& value) noexcept
   {
     m_values.fill(value);
     return *this;
   }
 
+  PUGS_INLINE Vector& operator=(const DataType&) = delete;
+
+  PUGS_INLINE Vector&
+  operator=(const ZeroType&) noexcept
+  {
+    m_values.fill(0);
+    return *this;
+  }
+
   template <typename DataType2>
   PUGS_INLINE Vector&
-  operator=(const Vector<DataType2>& x) noexcept
+  operator=(const Vector<DataType2>& x)
   {
     // ensures that DataType is the same as source DataType2
     static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
                   "Cannot assign Vector of different type");
-    // ensures that const is not lost through copy
-    static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
-                  "Cannot assign Vector of const to Vector of non-const");
+    // ensures that const is not lost through copy. If this point is
+    // reached data types only differ through their constness
+    static_assert((not std::is_const_v<DataType2>), "Cannot assign Vector of const to Vector of non-const");
 
-    if (m_deep_copies) {
-      copy_to(x.m_values, m_values);
-    } else {
-      m_values = x.m_values;
-    }
+    Assert(not m_deep_copies, "cannot assign value to a Vector of const that required deep copies");
+    m_values = x.m_values;
     return *this;
   }
 
@@ -185,7 +191,10 @@ class Vector   // LCOV_EXCL_LINE
   operator=(const Vector& x)
   {
     if (m_deep_copies) {
-      copy_to(x.m_values, m_values);
+      Assert(not std::is_const_v<DataType>, "cannot assign value to a Vector of const that required deep copies");
+      if constexpr (not std::is_const_v<DataType>) {
+        copy_to(x.m_values, m_values);
+      }
     } else {
       m_values = x.m_values;
     }
@@ -211,22 +220,20 @@ class Vector   // LCOV_EXCL_LINE
     static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
                   "Cannot assign Vector of different type");
     // ensures that const is not lost through copy
-    static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
+    static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
                   "Cannot assign Vector of const to Vector of non-const");
 
     m_values = x.m_values;
   }
 
-  explicit Vector(const Array<DataType>& values) : m_deep_copies{true}
-  {
-    m_values = values;
-  }
+  explicit Vector(const Array<DataType>& values) : m_values{values}, m_deep_copies{true} {}
 
   Vector(const Vector&) = default;
-
-  Vector(Vector&&) = default;
+  Vector(Vector&&)      = default;
 
   Vector(size_t size) : m_values{size}, m_deep_copies{false} {}
+  Vector() : m_deep_copies{false} {}
+
   ~Vector() = default;
 };
 
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 5b4d9cb57cdea7a374b68447fc7a79ba81562e84..e5171badb11fa509a020025ad38e7e7bbd96281d 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -241,14 +241,6 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   {
     BoundaryConditionList bc_list;
 
-    constexpr ItemType FaceType = [] {
-      if constexpr (Dimension > 1) {
-        return ItemType::face;
-      } else {
-        return ItemType::node;
-      }
-    }();
-
     for (const auto& bc_descriptor : bc_descriptor_list) {
       bool is_valid_boundary_condition = true;
 
@@ -288,6 +280,14 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
 
             bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
           } else {
+            constexpr ItemType FaceType = [] {
+              if constexpr (Dimension > 1) {
+                return ItemType::face;
+              } else {
+                return ItemType::node;
+              }
+            }();
+
             MeshFaceBoundary<Dimension> mesh_face_boundary =
               getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
 
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index e47e56575fe725f67f6c6815c34c42eac4ce8d44..e2a4e5c71c6938e96a2d2cf63318b1f7ed9c8f4c 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -60,7 +60,7 @@ class [[nodiscard]] Array
   PUGS_INLINE
   void fill(const DataType& data) const
   {
-    static_assert(not std::is_const<DataType>(), "Cannot modify Array of const");
+    static_assert(not std::is_const_v<DataType>, "Cannot modify Array of const");
 
     Kokkos::deep_copy(m_values, data);
   }
@@ -72,7 +72,7 @@ class [[nodiscard]] Array
     static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
                   "Cannot assign Array of different type");
     // ensures that const is not lost through copy
-    static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
+    static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
                   "Cannot assign Array of const to Array of non-const");
     m_values = array.m_values;
     return *this;
diff --git a/src/utils/PugsAssert.hpp b/src/utils/PugsAssert.hpp
index 7ac521a51f73df6c5129ddd6ba346bd06453c07f..7793a9fc6a4ef537545b7691510b3553477d1bdf 100644
--- a/src/utils/PugsAssert.hpp
+++ b/src/utils/PugsAssert.hpp
@@ -9,7 +9,7 @@
 #include <string>
 #include <tuple>
 
-class AssertError
+class AssertError : public std::runtime_error
 {
  private:
   const std::string m_file;
@@ -41,7 +41,8 @@ class AssertError
               const std::string& function,
               const std::tuple<Args...>& tuple_args,
               const std::string_view args_string)
-    : m_file{filename},
+    : std::runtime_error(""),
+      m_file{filename},
       m_line{line},
       m_function{function},
       m_test{[&] {
@@ -64,7 +65,11 @@ class AssertError
         }
       }()}
   {
-    ;
+    if (m_message.empty()) {
+      std::runtime_error::operator=(std::runtime_error(m_test.c_str()));
+    } else {
+      std::runtime_error::operator=(std::runtime_error(m_message.c_str()));
+    }
   }
 };
 
diff --git a/src/utils/RandomEngine.cpp b/src/utils/RandomEngine.cpp
index 508f77a39123492376a5141662f4b3daa41476f9..9b7efc162c4d2d7b0ff912689b9f718eae39af5b 100644
--- a/src/utils/RandomEngine.cpp
+++ b/src/utils/RandomEngine.cpp
@@ -41,13 +41,7 @@ RandomEngine::destroy()
 
 RandomEngine::RandomEngine()
 {
-  uint64_t random_seed = std::random_device{}();
-  parallel::broadcast(random_seed, 0);
-
-  m_random_engine = std::default_random_engine(random_seed);
-
-  std::cout << " * setting " << rang::fgB::green << "random seed" << rang::style::reset << " to " << rang::fgB::yellow
-            << random_seed << rang::style::reset << '\n';
+  this->resetRandomSeed();
 }
 
 void
@@ -62,5 +56,8 @@ RandomEngine::setRandomSeed(const uint64_t random_seed)
 void
 RandomEngine::resetRandomSeed()
 {
-  m_instance = std::unique_ptr<RandomEngine>(new RandomEngine);
+  uint64_t random_seed = std::random_device{}();
+  parallel::broadcast(random_seed, 0);
+
+  this->setRandomSeed(random_seed);
 }
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 7904eedc03c7a203667e5c604096a5cded6e136b..06967be54d53c7ff8e88d6ec38d675381cd8ad07 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -86,6 +86,7 @@ add_executable (unit_tests
   test_NameProcessor.cpp
   test_OStreamProcessor.cpp
   test_ParseError.cpp
+  test_PETScUtils.cpp
   test_PugsAssert.cpp
   test_PugsFunctionAdapter.cpp
   test_PugsUtils.cpp
@@ -113,6 +114,7 @@ add_executable (mpi_unit_tests
   test_ItemValueUtils.cpp
   test_Messenger.cpp
   test_Partitioner.cpp
+  test_RandomEngine.cpp
   test_SubItemValuePerItem.cpp
   test_SubItemArrayPerItem.cpp
   )
diff --git a/tests/mpi_test_main.cpp b/tests/mpi_test_main.cpp
index 51a01b0cfa99b36995b8849fc1fdb94360d971ad..4c94634173799a17519b5d417bfc914a99e2e740 100644
--- a/tests/mpi_test_main.cpp
+++ b/tests/mpi_test_main.cpp
@@ -9,6 +9,7 @@
 #include <mesh/SynchronizerManager.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PETScWrapper.hpp>
+#include <utils/RandomEngine.hpp>
 #include <utils/pugs_config.hpp>
 
 #include <MeshDataBaseForTests.hpp>
@@ -57,6 +58,7 @@ main(int argc, char* argv[])
       std::cout.setstate(std::ios::badbit);
 
       SynchronizerManager::create();
+      RandomEngine::create();
       MeshDataManager::create();
       DiamondDualConnectivityManager::create();
       DiamondDualMeshManager::create();
@@ -89,6 +91,7 @@ main(int argc, char* argv[])
       DiamondDualMeshManager::destroy();
       DiamondDualConnectivityManager::destroy();
       MeshDataManager::destroy();
+      RandomEngine::destroy();
       SynchronizerManager::destroy();
     }
   }
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
index 231783090756e52187b46da02a97d44cf2279ea9..904b59974def0ba927e8bb4ddc4e8f4f291625c5 100644
--- a/tests/test_BiCGStab.cpp
+++ b/tests/test_BiCGStab.cpp
@@ -51,7 +51,7 @@ TEST_CASE("BiCGStab", "[algebra]")
     Vector<double> b = A * x_exact;
 
     Vector<double> x{5};
-    x = 0;
+    x = zero;
 
     BiCGStab{A, x, b, 1e-12, 10, false};
     Vector error = x - x_exact;
@@ -96,7 +96,7 @@ TEST_CASE("BiCGStab", "[algebra]")
     Vector<double> b = A * x_exact;
 
     Vector<double> x{5};
-    x = 0;
+    x = zero;
 
     BiCGStab{A, x, b, 1e-12, 1, true};
     Vector error = x - x_exact;
diff --git a/tests/test_CG.cpp b/tests/test_CG.cpp
index dbc1c35952f57043c6ad0b23bf29883f4c678bdb..f0c14f24f9c416048aa05905123621deb011219f 100644
--- a/tests/test_CG.cpp
+++ b/tests/test_CG.cpp
@@ -49,7 +49,7 @@ TEST_CASE("CG", "[algebra]")
     Vector<double> b = A * x_exact;
 
     Vector<double> x{5};
-    x = 0;
+    x = zero;
 
     CG{A, x, b, 1e-12, 10, true};
     Vector error = x - x_exact;
@@ -78,10 +78,10 @@ TEST_CASE("CG", "[algebra]")
     REQUIRE(Sout.str() == Aout.str());
 
     Vector<double> b{5};
-    b = 0;
+    b = zero;
 
     Vector<double> x{5};
-    x = 0;
+    x = zero;
 
     CG{A, x, b, 1e-12, 10};
     REQUIRE(std::sqrt(dot(x, x)) == 0);
@@ -127,7 +127,7 @@ TEST_CASE("CG", "[algebra]")
     Vector<double> b = A * x_exact;
 
     Vector<double> x{5};
-    x = 0;
+    x = zero;
 
     CG{A, x, b, 1e-12, 1, false};
     Vector error = x - x_exact;
diff --git a/tests/test_CRSMatrix.cpp b/tests/test_CRSMatrix.cpp
index e0aa27b57c1a9fd773c99d7d96c6d5e47d79f74e..d0d9476ab0e47ff90874391c63f39a6321adbc3b 100644
--- a/tests/test_CRSMatrix.cpp
+++ b/tests/test_CRSMatrix.cpp
@@ -49,7 +49,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
     CRSMatrix<int> A{S.getCRSMatrix()};
 
     Vector<int> x(A.numberOfColumns());
-    x    = 0;
+    x    = zero;
     x[0] = 1;
 
     Vector<int> y = A * x;
@@ -59,7 +59,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
     REQUIRE(y[3] == 5);
     REQUIRE(y[4] == 0);
 
-    x    = 0;
+    x    = zero;
     x[1] = 2;
 
     y = A * x;
@@ -69,7 +69,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
     REQUIRE(y[3] == -6);
     REQUIRE(y[4] == 2);
 
-    x    = 0;
+    x    = zero;
     x[2] = -1;
 
     y = A * x;
@@ -79,7 +79,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
     REQUIRE(y[3] == 0);
     REQUIRE(y[4] == 0);
 
-    x    = 0;
+    x    = zero;
     x[3] = 3;
 
     y = A * x;
@@ -89,7 +89,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
     REQUIRE(y[3] == 0);
     REQUIRE(y[4] == 6);
 
-    x    = 0;
+    x    = zero;
     x[4] = 1;
 
     y = A * x;
diff --git a/tests/test_DenseMatrix.cpp b/tests/test_DenseMatrix.cpp
index 6e2470fa0f325fab1e17eb0f51c2ad34a61c318b..938ebefbda40aae6ff6d7a3c33f11355c51b23c1 100644
--- a/tests/test_DenseMatrix.cpp
+++ b/tests/test_DenseMatrix.cpp
@@ -20,6 +20,15 @@ TEST_CASE("DenseMatrix", "[algebra]")
     REQUIRE(A.numberOfColumns() == 3);
   }
 
+  SECTION("is square")
+  {
+    REQUIRE(not DenseMatrix<int>{2, 3}.isSquare());
+    REQUIRE(not DenseMatrix<int>{4, 3}.isSquare());
+    REQUIRE(DenseMatrix<int>{1, 1}.isSquare());
+    REQUIRE(DenseMatrix<int>{2, 2}.isSquare());
+    REQUIRE(DenseMatrix<int>{3, 3}.isSquare());
+  }
+
   SECTION("write access")
   {
     DenseMatrix<int> A{2, 3};
@@ -56,6 +65,28 @@ TEST_CASE("DenseMatrix", "[algebra]")
     REQUIRE(A(1, 2) == 2);
   }
 
+  SECTION("zero matrix")
+  {
+    DenseMatrix<int> A{2, 3};
+    A = zero;
+    for (size_t i = 0; i < A.numberOfRows(); ++i) {
+      for (size_t j = 0; j < A.numberOfColumns(); ++j) {
+        REQUIRE(A(i, j) == 0);
+      }
+    }
+  }
+
+  SECTION("identity matrix")
+  {
+    DenseMatrix<int> A{3, 3};
+    A = identity;
+    for (size_t i = 0; i < A.numberOfRows(); ++i) {
+      for (size_t j = 0; j < A.numberOfColumns(); ++j) {
+        REQUIRE(A(i, j) == (i == j));
+      }
+    }
+  }
+
   SECTION("copy constructor (shallow)")
   {
     DenseMatrix<int> A{2, 3};
@@ -341,4 +372,62 @@ TEST_CASE("DenseMatrix", "[algebra]")
     REQUIRE(C(1, 0) == 64);
     REQUIRE(C(1, 1) == 137);
   }
+
+#ifndef NDEBUG
+  SECTION("non square identity matrix")
+  {
+    DenseMatrix<int> A(2, 3);
+
+    REQUIRE_THROWS_WITH(A = identity, "identity must be a square matrix");
+  }
+
+  SECTION("invalid matrix vector product")
+  {
+    DenseMatrix<int> A(2, 3);
+    Vector<int> u(2);
+
+    REQUIRE_THROWS_WITH(A * u, "cannot compute matrix-vector product: incompatible sizes");
+  }
+
+  SECTION("invalid matrix product")
+  {
+    DenseMatrix<int> A(2, 3);
+    DenseMatrix<int> B(2, 3);
+
+    REQUIRE_THROWS_WITH(A * B, "cannot compute matrix product: incompatible sizes");
+  }
+
+  SECTION("invalid matrix substract")
+  {
+    DenseMatrix<int> A(2, 3);
+    DenseMatrix<int> B(3, 2);
+
+    REQUIRE_THROWS_WITH(A -= B, "cannot substract matrix: incompatible sizes");
+  }
+
+  SECTION("invalid matrix sum")
+  {
+    DenseMatrix<int> A(2, 3);
+    DenseMatrix<int> B(3, 2);
+
+    REQUIRE_THROWS_WITH(A + B, "cannot compute matrix sum: incompatible sizes");
+  }
+
+  SECTION("invalid matrix difference")
+  {
+    DenseMatrix<int> A(2, 3);
+    DenseMatrix<int> B(3, 2);
+
+    REQUIRE_THROWS_WITH(A - B, "cannot compute matrix difference: incompatible sizes");
+  }
+
+  SECTION("invalid indices")
+  {
+    DenseMatrix<int> A(2, 3);
+
+    REQUIRE_THROWS_WITH(A(2, 0), "cannot access element: invalid indices");
+    REQUIRE_THROWS_WITH(A(0, 3), "cannot access element: invalid indices");
+  }
+
+#endif   // NDEBUG
 }
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index a9b4cab1fe49efafa6a8ec4d38e34b7b947c6eff..f8bcd9972ca5b52a6ea5fe48d70eb1b8c2064107 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -195,7 +195,7 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.precond() = LSPrecond::none;
 
             Vector<double> x{5};
-            x = 0;
+            x = zero;
 
             LinearSolver solver{options};
 
@@ -222,7 +222,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::none;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -236,7 +236,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::diagonal;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -250,7 +250,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::incomplete_choleski;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -264,7 +264,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::amg;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -282,7 +282,7 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.precond() = LSPrecond::none;
 
             Vector<double> x{5};
-            x = 0;
+            x = zero;
 
             LinearSolver solver{options};
 
@@ -355,7 +355,7 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.precond() = LSPrecond::none;
 
             Vector<double> x{5};
-            x = 0;
+            x = zero;
 
             LinearSolver solver{options};
 
@@ -382,7 +382,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::none;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -396,7 +396,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::diagonal;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -410,7 +410,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::incomplete_LU;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -432,7 +432,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::none;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -446,7 +446,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::diagonal;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -468,7 +468,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::none;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -482,7 +482,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::diagonal;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -496,7 +496,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::incomplete_LU;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -514,7 +514,7 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.precond() = LSPrecond::none;
 
             Vector<double> x{5};
-            x = 0;
+            x = zero;
 
             LinearSolver solver{options};
 
@@ -588,7 +588,7 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.precond() = LSPrecond::none;
 
             Vector<double> x{5};
-            x = 0;
+            x = zero;
 
             LinearSolver solver{options};
 
@@ -615,7 +615,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::none;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -629,7 +629,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::diagonal;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -643,7 +643,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::incomplete_choleski;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -657,7 +657,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::amg;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -675,7 +675,7 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.precond() = LSPrecond::none;
 
             Vector<double> x{5};
-            x = 0;
+            x = zero;
 
             LinearSolver solver{options};
 
@@ -743,7 +743,7 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.precond() = LSPrecond::none;
 
             Vector<double> x{5};
-            x = 0;
+            x = zero;
 
             LinearSolver solver{options};
 
@@ -770,7 +770,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::none;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -784,7 +784,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::diagonal;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -798,7 +798,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::incomplete_LU;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -820,7 +820,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::none;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -834,7 +834,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::diagonal;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -856,7 +856,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::none;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -870,7 +870,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::diagonal;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -884,7 +884,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               options.precond() = LSPrecond::incomplete_LU;
 
               Vector<double> x{5};
-              x = 0;
+              x = zero;
 
               LinearSolver solver{options};
 
@@ -902,7 +902,7 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.precond() = LSPrecond::none;
 
             Vector<double> x{5};
-            x = 0;
+            x = zero;
 
             LinearSolver solver{options};
 
diff --git a/tests/test_Messenger.cpp b/tests/test_Messenger.cpp
index eebd5dc939c23c7875c3b92ff697cb1f00afca63..4f76e250ddab60d81e1d49eacb1bcfc23425e317 100644
--- a/tests/test_Messenger.cpp
+++ b/tests/test_Messenger.cpp
@@ -1,6 +1,7 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_all.hpp>
 
+#include <algebra/TinyVector.hpp>
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
 
@@ -143,6 +144,14 @@ TEST_CASE("Messenger", "[mpi]")
 
     const bool or_value_2 = parallel::allReduceOr(parallel::rank() > 0);
     REQUIRE(or_value_2 == (parallel::size() > 1));
+
+    const size_t sum_value = parallel::allReduceSum(parallel::rank() + 1);
+    REQUIRE(sum_value == parallel::size() * (parallel::size() + 1) / 2);
+
+    const TinyVector<2, size_t> sum_tiny_vector =
+      parallel::allReduceSum(TinyVector<2, size_t>(parallel::rank() + 1, 1));
+    REQUIRE(
+      (sum_tiny_vector == TinyVector<2, size_t>{parallel::size() * (parallel::size() + 1) / 2, parallel::size()}));
   }
 
   SECTION("all to all")
diff --git a/tests/test_PETScUtils.cpp b/tests/test_PETScUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2e0daa8519241fde2f7fd1749b40398977b839de
--- /dev/null
+++ b/tests/test_PETScUtils.cpp
@@ -0,0 +1,117 @@
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_PETSC
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <algebra/PETScUtils.hpp>
+
+#include <algebra/CRSMatrixDescriptor.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PETScUtils", "[algebra]")
+{
+  SECTION("PETScAijMatrixEmbedder")
+  {
+    SECTION("from TinyMatrix")
+    {
+      TinyMatrix<3> A{1, 2, 3, 4, 5, 6, 7, 8, 9};
+      PETScAijMatrixEmbedder petsc_A{A};
+
+      for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          double value = 0;
+          MatGetValue(petsc_A, i, j, &value);
+          REQUIRE(value == 1 + i * 3 + j);
+        }
+      }
+    }
+
+    SECTION("from DenseMatrix")
+    {
+      DenseMatrix<double> A(3, 3);
+      for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          A(i, j) = 1 + i * 3 + j;
+        }
+      }
+
+      PETScAijMatrixEmbedder petsc_A{A};
+
+      REQUIRE(petsc_A.numberOfRows() == 3);
+      REQUIRE(petsc_A.numberOfColumns() == 3);
+
+      for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          double value;
+          MatGetValue(petsc_A, i, j, &value);
+          REQUIRE(value == 1 + i * 3 + j);
+        }
+      }
+    }
+
+    SECTION("from DenseMatrix [non-square]")
+    {
+      DenseMatrix<double> A(4, 3);
+      for (size_t i = 0; i < 4; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          A(i, j) = 1 + i * 3 + j;
+        }
+      }
+
+      PETScAijMatrixEmbedder petsc_A{A};
+
+      REQUIRE(petsc_A.numberOfRows() == 4);
+      REQUIRE(petsc_A.numberOfColumns() == 3);
+
+      for (size_t i = 0; i < 4; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          double value = 0;
+          MatGetValue(petsc_A, i, j, &value);
+          REQUIRE(value == 1 + i * 3 + j);
+        }
+      }
+    }
+  }
+
+  SECTION("from CRSMatrix")
+  {
+    Array<int> non_zeros(4);
+    non_zeros[0] = 1;
+    non_zeros[1] = 2;
+    non_zeros[2] = 3;
+    non_zeros[3] = 2;
+
+    CRSMatrixDescriptor<double, int> A(4, 3, non_zeros);
+
+    A(0, 0) = 1;
+    A(1, 0) = 2;
+    A(1, 2) = 3;
+    A(2, 0) = 4;
+    A(2, 1) = 5;
+    A(2, 2) = 6;
+    A(3, 1) = 7;
+    A(3, 2) = 8;
+
+    PETScAijMatrixEmbedder petsc_A{A.getCRSMatrix()};
+
+    auto get_value = [&](int i, int j) {
+      double value;
+      MatGetValue(petsc_A, i, j, &value);
+      return value;
+    };
+
+    REQUIRE(get_value(0, 0) == 1);
+    REQUIRE(get_value(1, 0) == 2);
+    REQUIRE(get_value(1, 2) == 3);
+    REQUIRE(get_value(2, 0) == 4);
+    REQUIRE(get_value(2, 1) == 5);
+    REQUIRE(get_value(2, 2) == 6);
+    REQUIRE(get_value(3, 1) == 7);
+    REQUIRE(get_value(3, 2) == 8);
+  }
+}
+
+#endif   // PUGS_HAS_PETSC
diff --git a/tests/test_RandomEngine.cpp b/tests/test_RandomEngine.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..241b40dbe958b68ccf631fce2291affca9e1e29c
--- /dev/null
+++ b/tests/test_RandomEngine.cpp
@@ -0,0 +1,52 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <utils/Messenger.hpp>
+#include <utils/RandomEngine.hpp>
+
+#include <utils/pugs_config.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("RandomEngine", "[random]")
+{
+  SECTION("current state")
+  {
+    RandomEngine& random_engine = RandomEngine::instance();
+    REQUIRE(isSynchronized(random_engine));
+  }
+
+  SECTION("set seed")
+  {
+    RandomEngine& random_engine = RandomEngine::instance();
+    random_engine.setRandomSeed(1402339680);
+
+    REQUIRE(isSynchronized(random_engine));
+
+    REQUIRE(random_engine.getCurrentSeed() == 1402339680);
+  }
+
+  SECTION("reset seed")
+  {
+    RandomEngine& random_engine = RandomEngine::instance();
+    random_engine.resetRandomSeed();
+
+    REQUIRE(isSynchronized(random_engine));
+  }
+
+  SECTION("de-synchronize seeds")
+  {
+    RandomEngine& random_engine = RandomEngine::instance();
+    random_engine.resetRandomSeed();
+    auto& engine = random_engine.engine();
+
+    for (size_t i = 0; i < parallel::rank(); ++i) {
+      engine();
+    }
+
+    REQUIRE(((parallel::size() == 1) or not isSynchronized(random_engine)));
+
+    random_engine.resetRandomSeed();
+    REQUIRE(isSynchronized(random_engine));
+  }
+}
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index 3a03b47f71e9e5491485053870f93ab45349d942..31ded1e245bc026a0a31405ddccca0413717c402 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -24,6 +24,7 @@ TEST_CASE("TinyMatrix", "[algebra]")
            (A(1, 2) == 6) and (A(2, 0) == 7) and (A(2, 1) == 8) and (A(2, 2) == 9)));
 
   TinyMatrix<3, int> B(6, 5, 3, 8, 34, 6, 35, 6, 7);
+
   SECTION("checking for opposed matrix")
   {
     const TinyMatrix<3, int> minus_A = -A;
@@ -31,6 +32,7 @@ TEST_CASE("TinyMatrix", "[algebra]")
              (minus_A(1, 1) == -5) and (minus_A(1, 2) == -6) and (minus_A(2, 0) == -7) and (minus_A(2, 1) == -8) and
              (minus_A(2, 2) == -9)));
   }
+
   SECTION("checking for equality and difference tests")
   {
     const TinyMatrix<3, int> copy_A = A;
@@ -209,14 +211,24 @@ TEST_CASE("TinyMatrix", "[algebra]")
     REQUIRE(TinyMatrix<1>{}.numberOfRows() == 1);
     REQUIRE(TinyMatrix<1>{}.numberOfColumns() == 1);
     REQUIRE(TinyMatrix<1>{}.dimension() == 1);
+    REQUIRE(TinyMatrix<1>{}.numberOfValues() == 1);
 
     REQUIRE(TinyMatrix<2>{}.numberOfRows() == 2);
     REQUIRE(TinyMatrix<2>{}.numberOfColumns() == 2);
-    REQUIRE(TinyMatrix<2>{}.dimension() == 4);
+    REQUIRE(TinyMatrix<2>{}.numberOfValues() == 4);
 
     REQUIRE(TinyMatrix<3>{}.numberOfRows() == 3);
     REQUIRE(TinyMatrix<3>{}.numberOfColumns() == 3);
     REQUIRE(TinyMatrix<3>{}.dimension() == 9);
+    REQUIRE(TinyMatrix<3>{}.numberOfValues() == 9);
+  }
+
+  SECTION("is square matrix")
+  {
+    REQUIRE(TinyMatrix<1>{}.isSquare());
+    REQUIRE(TinyMatrix<2>{}.isSquare());
+    REQUIRE(TinyMatrix<3>{}.isSquare());
+    REQUIRE(TinyMatrix<4>{}.isSquare());
   }
 
   SECTION("checking for matrices output")
diff --git a/tests/test_Vector.cpp b/tests/test_Vector.cpp
index 89fedcec868e026cfdf569eb8667b57b730a6e16..5a8f2f77ced534e9f1e2029e228dd2ed6fc372fd 100644
--- a/tests/test_Vector.cpp
+++ b/tests/test_Vector.cpp
@@ -37,7 +37,7 @@ TEST_CASE("Vector", "[algebra]")
   SECTION("fill")
   {
     Vector<int> x{5};
-    x = 2;
+    x.fill(2);
 
     REQUIRE(x[0] == 2);
     REQUIRE(x[1] == 2);
@@ -63,6 +63,59 @@ TEST_CASE("Vector", "[algebra]")
     REQUIRE(y[4] == 4);
   }
 
+  SECTION("copy constructor (move)")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    const Vector<int> y = std::move(x);
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 1);
+    REQUIRE(y[2] == 2);
+    REQUIRE(y[3] == 3);
+    REQUIRE(y[4] == 4);
+  }
+
+  SECTION("copy (shallow)")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    Vector<int> y{5};
+    y = x;
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 1);
+    REQUIRE(y[2] == 2);
+    REQUIRE(y[3] == 3);
+    REQUIRE(y[4] == 4);
+
+    x[0] = 14;
+    x[1] = 13;
+    x[2] = 12;
+    x[3] = 11;
+    x[4] = 10;
+
+    REQUIRE(x[0] == 14);
+    REQUIRE(x[1] == 13);
+    REQUIRE(x[2] == 12);
+    REQUIRE(x[3] == 11);
+    REQUIRE(x[4] == 10);
+
+    REQUIRE(y[0] == 14);
+    REQUIRE(y[1] == 13);
+    REQUIRE(y[2] == 12);
+    REQUIRE(y[3] == 11);
+    REQUIRE(y[4] == 10);
+  }
+
   SECTION("copy (deep)")
   {
     Vector<int> x{5};
@@ -72,7 +125,8 @@ TEST_CASE("Vector", "[algebra]")
     x[3] = 3;
     x[4] = 4;
 
-    const Vector<int> y = copy(x);
+    Vector<int> y{5};
+    y = copy(x);
 
     x[0] = 14;
     x[1] = 13;
@@ -93,6 +147,44 @@ TEST_CASE("Vector", "[algebra]")
     REQUIRE(x[4] == 10);
   }
 
+  SECTION("copy to const (shallow)")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    Vector<const int> y;
+    y = x;
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 1);
+    REQUIRE(y[2] == 2);
+    REQUIRE(y[3] == 3);
+    REQUIRE(y[4] == 4);
+
+    Vector<int> z{5};
+    z[0] = 14;
+    z[1] = 13;
+    z[2] = 12;
+    z[3] = 11;
+    z[4] = 10;
+
+    y = z;
+    REQUIRE(z[0] == 14);
+    REQUIRE(z[1] == 13);
+    REQUIRE(z[2] == 12);
+    REQUIRE(z[3] == 11);
+    REQUIRE(z[4] == 10);
+
+    REQUIRE(y[0] == 14);
+    REQUIRE(y[1] == 13);
+    REQUIRE(y[2] == 12);
+    REQUIRE(y[3] == 11);
+    REQUIRE(y[4] == 10);
+  }
+
   SECTION("self scalar multiplication")
   {
     Vector<int> x{5};
@@ -266,4 +358,141 @@ TEST_CASE("Vector", "[algebra]")
     REQUIRE(z[3] == 5);
     REQUIRE(z[4] == 4);
   }
+
+  SECTION("output")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    std::ostringstream vector_ost;
+    vector_ost << x;
+    std::ostringstream ref_ost;
+    for (size_t i = 0; i < x.size(); ++i) {
+      ref_ost << ' ' << x[i];
+    }
+    REQUIRE(vector_ost.str() == ref_ost.str());
+  }
+
+  SECTION("Vector view from Array")
+  {
+    Array<double> values{5};
+    for (size_t i = 0; i < values.size(); ++i) {
+      values[i] = 2 * i + 1;
+    }
+
+    Vector vector{values};
+
+    SECTION("Constructor")
+    {
+      REQUIRE(vector[0] == 1);
+      REQUIRE(vector[1] == 3);
+      REQUIRE(vector[2] == 5);
+      REQUIRE(vector[3] == 7);
+      REQUIRE(vector[4] == 9);
+    }
+
+    SECTION("Affectation (implicit deep copy)")
+    {
+      Vector<double> x{vector.size()};
+      for (size_t i = 0; i < x.size(); ++i) {
+        x[i] = i;
+      }
+
+      vector = x;
+      REQUIRE(values[0] == 0);
+      REQUIRE(values[1] == 1);
+      REQUIRE(values[2] == 2);
+      REQUIRE(values[3] == 3);
+      REQUIRE(values[4] == 4);
+
+      for (size_t i = 0; i < x.size(); ++i) {
+        x[i] = x.size() - i;
+      }
+
+      REQUIRE(values[0] == 0);
+      REQUIRE(values[1] == 1);
+      REQUIRE(values[2] == 2);
+      REQUIRE(values[3] == 3);
+      REQUIRE(values[4] == 4);
+
+      vector = std::move(x);
+      x.fill(2);
+
+      REQUIRE(values[0] == 5);
+      REQUIRE(values[1] == 4);
+      REQUIRE(values[2] == 3);
+      REQUIRE(values[3] == 2);
+      REQUIRE(values[4] == 1);
+    }
+  }
+
+#ifndef NDEBUG
+  SECTION("invalid dot product")
+  {
+    Vector<int> x{5};
+    Vector<int> y{4};
+
+    REQUIRE_THROWS_WITH(dot(x, y), "cannot compute dot product: incompatible vector sizes");
+  }
+
+  SECTION("invalid substract")
+  {
+    Vector<int> x{5};
+    Vector<int> y{4};
+
+    REQUIRE_THROWS_WITH(x -= y, "cannot substract vector: incompatible sizes");
+  }
+
+  SECTION("invalid add")
+  {
+    Vector<int> x{5};
+    Vector<int> y{4};
+
+    REQUIRE_THROWS_WITH(x += y, "cannot add vector: incompatible sizes");
+  }
+
+  SECTION("invalid difference")
+  {
+    Vector<int> x{5};
+    Vector<int> y{4};
+
+    REQUIRE_THROWS_WITH(x - y, "cannot compute vector difference: incompatible sizes");
+  }
+
+  SECTION("invalid sum")
+  {
+    Vector<int> x{5};
+    Vector<int> y{4};
+
+    REQUIRE_THROWS_WITH(x + y, "cannot compute vector sum: incompatible sizes");
+  }
+
+  SECTION("Affectation to Vector of const values built from Array")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    Array<int> y_values{5};
+    y_values[0] = 0;
+    y_values[1] = 1;
+    y_values[2] = 2;
+    y_values[3] = 3;
+    y_values[4] = 4;
+    Vector<const int> y{y_values};
+
+    REQUIRE_THROWS_WITH(y = x, "cannot assign value to a Vector of const that required deep copies");
+
+    Vector<const int> z = x;
+
+    REQUIRE_THROWS_WITH(y = z, "cannot assign value to a Vector of const that required deep copies");
+  }
+#endif   // NDEBUG
 }
diff --git a/tests/test_main.cpp b/tests/test_main.cpp
index b5316614d10c21cbe13bb5540cc0a13d7e4b4b24..9c03ba9213ece8947688036ca7d5554c8e38a279 100644
--- a/tests/test_main.cpp
+++ b/tests/test_main.cpp
@@ -9,6 +9,7 @@
 #include <mesh/SynchronizerManager.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PETScWrapper.hpp>
+#include <utils/RandomEngine.hpp>
 #include <utils/SLEPcWrapper.hpp>
 
 #include <MeshDataBaseForTests.hpp>
@@ -34,6 +35,7 @@ main(int argc, char* argv[])
       std::cout.setstate(std::ios::badbit);
 
       SynchronizerManager::create();
+      RandomEngine::create();
       MeshDataManager::create();
       DiamondDualConnectivityManager::create();
       DiamondDualMeshManager::create();
@@ -51,6 +53,7 @@ main(int argc, char* argv[])
       DiamondDualMeshManager::destroy();
       DiamondDualConnectivityManager::destroy();
       MeshDataManager::destroy();
+      RandomEngine::destroy();
       SynchronizerManager::destroy();
     }
   }