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/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_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_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
 }