diff --git a/src/algebra/CRSMatrixDescriptor.hpp b/src/algebra/CRSMatrixDescriptor.hpp
index cbf82a904e85dcb8aa735e978612a4baec0affa7..8f123896d486b7d6b67be1350a96183423b4235e 100644
--- a/src/algebra/CRSMatrixDescriptor.hpp
+++ b/src/algebra/CRSMatrixDescriptor.hpp
@@ -35,6 +35,7 @@ class CRSMatrixDescriptor
     Array<IndexType> row_map(m_nb_rows + 1);
     row_map[0] = 0;
     for (IndexType i = 0; i < m_nb_rows; ++i) {
+      Assert(non_zeros[i] >= 0);
       row_map[i + 1] = row_map[i] + non_zeros[i];
     }
 
@@ -83,7 +84,7 @@ class CRSMatrixDescriptor
       const auto& overflow_row = A.m_overflow_per_row[i];
       os << i << "|";
 
-      if (A.m_nb_values_per_row[i] > 0) {
+      if (A.m_nb_values_per_row[i] + overflow_row.size() > 0) {
         IndexType j     = 0;
         auto i_overflow = overflow_row.begin();
         while (j < A.m_nb_values_per_row[i] or i_overflow != overflow_row.end()) {
@@ -207,7 +208,7 @@ class CRSMatrixDescriptor
       for (IndexType i = 0; i < m_nb_rows; ++i) {
         const auto& overflow_row = m_overflow_per_row[i];
 
-        if (m_nb_values_per_row[i] > 0) {
+        if (m_nb_values_per_row[i] + overflow_row.size() > 0) {
           IndexType j     = 0;
           auto i_overflow = overflow_row.begin();
           while (j < m_nb_values_per_row[i] or i_overflow != overflow_row.end()) {
@@ -257,6 +258,13 @@ class CRSMatrixDescriptor
   {
     m_nb_values_per_row.fill(0);
     m_values.fill(0);
+
+    // Diagonal is always set to fulfill PETSc's requirements
+    if (m_nb_columns == m_nb_rows) {
+      for (IndexType i = 0; i < m_nb_rows; ++i) {
+        this->operator()(i, i) = DataType{0};
+      }
+    }
   }
 
   ~CRSMatrixDescriptor() = default;
diff --git a/tests/test_CG.cpp b/tests/test_CG.cpp
index b76d90a047f2aab41e0576ab6f69ed5b133b474e..dbc1c35952f57043c6ad0b23bf29883f4c678bdb 100644
--- a/tests/test_CG.cpp
+++ b/tests/test_CG.cpp
@@ -62,8 +62,21 @@ TEST_CASE("CG", "[algebra]")
     non_zeros.fill(0);
     CRSMatrixDescriptor<double> S{5, 5, non_zeros};
 
+    std::ostringstream Sout;
+    Sout << S;
+    REQUIRE(Sout.str() == R"(0| 0:0
+1| 1:0
+2| 2:0
+3| 3:0
+4| 4:0
+)");
+
     CRSMatrix<double> A{S.getCRSMatrix()};
 
+    std::ostringstream Aout;
+    Aout << A;
+    REQUIRE(Sout.str() == Aout.str());
+
     Vector<double> b{5};
     b = 0;
 
diff --git a/tests/test_CRSMatrixDescriptor.cpp b/tests/test_CRSMatrixDescriptor.cpp
index 2bf59af22201ee9a8f164d2d6344b16357c56303..2b311d84b06dff839b3376d8c72c4439bc6855e0 100644
--- a/tests/test_CRSMatrixDescriptor.cpp
+++ b/tests/test_CRSMatrixDescriptor.cpp
@@ -162,5 +162,13 @@ TEST_CASE("CRSMatrixDescriptor", "[algebra]")
     REQUIRE_THROWS_AS(S(0, 5) = 1, AssertError);
   }
 
+  SECTION("negative number of non-zeros")
+  {
+    Array<int> non_zeros(2);
+    non_zeros[0] = 1;
+    non_zeros[1] = -2;
+    REQUIRE_THROWS_AS((CRSMatrixDescriptor<int>{2, 4, non_zeros}), AssertError);
+  }
+
 #endif   // NDEBUG
 }