Merge branch 'feature/linear-algebra' into 'develop'

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

See merge request !95

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;

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;
 

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
 }