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