Fix the case of incorrect 0 number of values prediction in rows

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;