diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp
index d9537f0e488f39bedf8849c299f06221e62a7dc9..f8bbf531304d3952b7e0de052eaf055d56997869 100644
--- a/src/algebra/Vector.hpp
+++ b/src/algebra/Vector.hpp
@@ -9,7 +9,7 @@
 #include <Array.hpp>
 
 template <typename DataType>
-class Vector
+class Vector   // LCOV_EXCL_LINE
 {
  public:
   using data_type  = DataType;
@@ -30,16 +30,8 @@ class Vector
     return x_copy;
   }
 
-  friend std::ostream&
-  operator<<(std::ostream& os, const Vector<DataType>& x)
-  {
-    for (index_type i = 0; i < x.size(); ++i) {
-      os << i << " : " << x[i] << '\n';
-    }
-    return os;
-  }
-
-  friend Vector operator*(const DataType& a, const Vector& x)
+  template <typename DataType2>
+  friend Vector operator*(const DataType2& a, const Vector& x)
   {
     Vector y = copy(x);
     return y *= a;
@@ -59,17 +51,17 @@ class Vector
     return sum;
   }
 
-  PUGS_INLINE
-  Vector&
-  operator/=(const DataType& a)
+  template <typename DataType2>
+  PUGS_INLINE Vector&
+  operator/=(const DataType2& a)
   {
-    const DataType inv_a = 1. / a;
+    const auto inv_a = 1. / a;
     return (*this) *= inv_a;
   }
 
-  PUGS_INLINE
-  Vector&
-  operator*=(const DataType& a)
+  template <typename DataType2>
+  PUGS_INLINE Vector&
+  operator*=(const DataType2& a)
   {
     parallel_for(
       this->size(), PUGS_LAMBDA(const index_type& i) { m_values[i] *= a; });
@@ -165,7 +157,7 @@ class Vector
 
   Vector(Vector&&) = default;
 
-  Vector(const size_t& size) : m_values(size) {}
+  Vector(const size_t& size) : m_values{size} {}
   ~Vector() = default;
 };
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index dad75a12071765ed09cba7816a4a8f71469c18fc..41f931d731e5d4889cf4ff9c72eb53ff36604ab9 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -10,6 +10,7 @@ add_executable (unit_tests
   test_RevisionInfo.cpp
   test_TinyMatrix.cpp
   test_TinyVector.cpp
+  test_Vector.cpp
   )
 
 add_executable (mpi_unit_tests
diff --git a/tests/test_Vector.cpp b/tests/test_Vector.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4384acd0660e1897b4aa0de54926b724c21c08dc
--- /dev/null
+++ b/tests/test_Vector.cpp
@@ -0,0 +1,283 @@
+#include <catch2/catch.hpp>
+
+#include <PugsAssert.hpp>
+#include <Vector.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class Vector<int>;
+
+TEST_CASE("Vector", "[algebra]")
+{
+  SECTION("size")
+  {
+    Vector<int> x{5};
+    REQUIRE(x.size() == 5);
+  }
+
+  SECTION("write access")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    REQUIRE(x[0] == 0);
+    REQUIRE(x[1] == 1);
+    REQUIRE(x[2] == 2);
+    REQUIRE(x[3] == 3);
+    REQUIRE(x[4] == 4);
+  }
+
+  SECTION("fill")
+  {
+    Vector<int> x{5};
+    x = 2;
+
+    REQUIRE(x[0] == 2);
+    REQUIRE(x[1] == 2);
+    REQUIRE(x[2] == 2);
+    REQUIRE(x[3] == 2);
+    REQUIRE(x[4] == 2);
+  }
+
+  SECTION("copy constructor (shallow)")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    const Vector<int> y = x;
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 1);
+    REQUIRE(y[2] == 2);
+    REQUIRE(y[3] == 3);
+    REQUIRE(y[4] == 4);
+  }
+
+  SECTION("copy (deep)")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    const Vector<int> y = copy(x);
+
+    x[0] = 14;
+    x[1] = 13;
+    x[2] = 12;
+    x[3] = 11;
+    x[4] = 10;
+
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 1);
+    REQUIRE(y[2] == 2);
+    REQUIRE(y[3] == 3);
+    REQUIRE(y[4] == 4);
+
+    REQUIRE(x[0] == 14);
+    REQUIRE(x[1] == 13);
+    REQUIRE(x[2] == 12);
+    REQUIRE(x[3] == 11);
+    REQUIRE(x[4] == 10);
+  }
+
+  SECTION("self scalar multiplication")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    x *= 2;
+
+    REQUIRE(x[0] == 0);
+    REQUIRE(x[1] == 2);
+    REQUIRE(x[2] == 4);
+    REQUIRE(x[3] == 6);
+    REQUIRE(x[4] == 8);
+  }
+
+  SECTION("left scalar multiplication")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    const Vector<int> y = 2 * x;
+
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 2);
+    REQUIRE(y[2] == 4);
+    REQUIRE(y[3] == 6);
+    REQUIRE(y[4] == 8);
+  }
+
+  SECTION("dot product")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    Vector<int> y{5};
+    y[0] = 7;
+    y[1] = 3;
+    y[2] = 4;
+    y[3] = 2;
+    y[4] = 8;
+
+    const int dot = (x, y);
+    REQUIRE(dot == 49);
+  }
+
+  SECTION("self scalar division")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    x /= 2;
+
+    REQUIRE(x[0] == 1);
+    REQUIRE(x[1] == 1);
+    REQUIRE(x[2] == 2);
+    REQUIRE(x[3] == 3);
+    REQUIRE(x[4] == 4);
+  }
+
+  SECTION("self minus")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    Vector<int> y{5};
+    y[0] = 3;
+    y[1] = 8;
+    y[2] = 6;
+    y[3] = 2;
+    y[4] = 4;
+
+    x -= y;
+
+    REQUIRE(x[0] == -1);
+    REQUIRE(x[1] == -5);
+    REQUIRE(x[2] == -1);
+    REQUIRE(x[3] == 5);
+    REQUIRE(x[4] == 4);
+  }
+
+  SECTION("self sum")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    Vector<int> y{5};
+    y[0] = 3;
+    y[1] = 8;
+    y[2] = 6;
+    y[3] = 2;
+    y[4] = 4;
+
+    x += y;
+
+    REQUIRE(x[0] == 5);
+    REQUIRE(x[1] == 11);
+    REQUIRE(x[2] == 11);
+    REQUIRE(x[3] == 9);
+    REQUIRE(x[4] == 12);
+  }
+
+  SECTION("sum")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    Vector<int> y{5};
+    y[0] = 3;
+    y[1] = 8;
+    y[2] = 6;
+    y[3] = 2;
+    y[4] = 4;
+
+    Vector z = x + y;
+
+    REQUIRE(z[0] == 5);
+    REQUIRE(z[1] == 11);
+    REQUIRE(z[2] == 11);
+    REQUIRE(z[3] == 9);
+    REQUIRE(z[4] == 12);
+  }
+
+  SECTION("difference")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    Vector<int> y{5};
+    y[0] = 3;
+    y[1] = 8;
+    y[2] = 6;
+    y[3] = 2;
+    y[4] = 4;
+
+    Vector z = x - y;
+
+    REQUIRE(z[0] == -1);
+    REQUIRE(z[1] == -5);
+    REQUIRE(z[2] == -1);
+    REQUIRE(z[3] == 5);
+    REQUIRE(z[4] == 4);
+  }
+
+  //   REQUIRE(l2Norm(TinyVector<2, double>(3, 4)) == Approx(5).epsilon(1E-14));
+
+  //   SECTION("checking for cross product")
+  //   {
+  //     const TinyVector<3, int> a(1, -2, 4);
+  //     const TinyVector<3, int> b(3, 1, 6);
+  //     REQUIRE(crossProduct(a, b) == TinyVector<3, int>(-16, 6, 7));
+  //   }
+
+  // #ifndef NDEBUG
+  //   SECTION("checking for bounds validation")
+  //   {
+  //     REQUIRE_THROWS_AS(x[4] = 0, AssertError);
+  //     const TinyVector<3, int>& const_x = x;
+  //     REQUIRE_THROWS_AS(const_x[-1], AssertError);
+  //   }
+  // #endif   // NDEBUG
+}