Add a bunch of improvements and test for PCG

- Vector<const T> should now be supported
- CRSMatrix data is constant

src/algebra/CRSMatrix.hpp

@@ -14,11 +14,13 @@
 template <typename DataType = double, typename IndexType = size_t>
 class CRSMatrix
 {
+  using MutableDataType = std::remove_const_t<DataType>;
+
  private:
   using HostMatrix = Kokkos::StaticCrsGraph<IndexType, Kokkos::HostSpace>;
 
   HostMatrix m_host_matrix;
-  Array<DataType> m_values;
+  Array<const DataType> m_values;
 
  public:
   PUGS_INLINE
@@ -28,16 +30,22 @@ class CRSMatrix
     return m_host_matrix.numRows();
   }
 
-  Vector<DataType> operator*(const Vector<DataType>& x) const
+  template <typename DataType2>
+  Vector<MutableDataType> operator*(const Vector<DataType2>& x) const
   {
-    Vector<DataType> Ax{m_host_matrix.numRows()};
+    static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+                  "Cannot multiply matrix and vector of different type");
+
+    Assert(x.size() == m_host_matrix.numRows());
+
+    Vector<MutableDataType> Ax{m_host_matrix.numRows()};
     auto host_row_map = m_host_matrix.row_map;
 
     parallel_for(
       m_host_matrix.numRows(), PUGS_LAMBDA(const IndexType& i_row) {
         const auto row_begin = host_row_map(i_row);
         const auto row_end   = host_row_map(i_row + 1);
-        DataType sum{0};
+        MutableDataType sum{0};
         for (IndexType j = row_begin; j < row_end; ++j) {
           sum += m_values[j] * x[m_host_matrix.entries(j)];
         }

src/algebra/PCG.hpp

 #ifndef PCG_HPP
 #define PCG_HPP
 
+#include <iomanip>
+#include <iostream>
+
 struct PCG
 {
   template <typename VectorType, typename MatrixType, typename PreconditionerType>

src/algebra/Vector.hpp

@@ -19,29 +19,40 @@ class Vector   // LCOV_EXCL_LINE
   Array<DataType> m_values;
   static_assert(std::is_same_v<typename decltype(m_values)::index_type, index_type>);
 
+  // Allows const version to access our data
+  friend Vector<std::add_const_t<DataType>>;
+
  public:
-  friend PUGS_INLINE Vector
+  friend Vector<std::remove_const_t<DataType>>
   copy(const Vector& x)
   {
     auto values = copy(x.m_values);
 
-    Vector x_copy{0};
+    Vector<std::remove_const_t<DataType>> x_copy{0};
     x_copy.m_values = values;
     return x_copy;
   }
 
-  template <typename DataType2>
-  friend Vector operator*(const DataType2& a, const Vector& x)
+  friend Vector operator*(const DataType& a, const Vector& x)
   {
-    Vector y = copy(x);
+    Vector<std::remove_const_t<DataType>> y = copy(x);
     return y *= a;
   }
 
+  template <typename DataType2>
   PUGS_INLINE
-  DataType
-  operator,(const Vector& y)
+  auto
+  operator,(const Vector<DataType2>& y)
   {
-    DataType sum = 0;
+    Assert(this->size() == y.size());
+    // Quite ugly, TODO: fix me in C++20
+    auto promoted = [] {
+      DataType a{0};
+      DataType2 b{0};
+      return a * b;
+    }();
+
+    decltype(promoted) sum = 0;
 
     // Does not use parallel_for to preserve sum order
     for (index_type i = 0; i < this->size(); ++i) {
@@ -65,56 +76,55 @@ class Vector   // LCOV_EXCL_LINE
   {
     parallel_for(
       this->size(), PUGS_LAMBDA(const index_type& i) { m_values[i] *= a; });
-
     return *this;
   }
 
-  PUGS_INLINE
-  Vector&
-  operator-=(const Vector& y)
+  template <typename DataType2>
+  PUGS_INLINE Vector&
+  operator-=(const Vector<DataType2>& y)
   {
     Assert(this->size() == y.size());
 
     parallel_for(
-      this->size(), PUGS_LAMBDA(const index_type& i) { m_values[i] -= y.m_values[i]; });
+      this->size(), PUGS_LAMBDA(const index_type& i) { m_values[i] -= y[i]; });
 
     return *this;
   }
 
-  PUGS_INLINE
-  Vector&
-  operator+=(const Vector& y)
+  template <typename DataType2>
+  PUGS_INLINE Vector&
+  operator+=(const Vector<DataType2>& y)
   {
     Assert(this->size() == y.size());
 
     parallel_for(
-      this->size(), PUGS_LAMBDA(const index_type& i) { m_values[i] += y.m_values[i]; });
+      this->size(), PUGS_LAMBDA(const index_type& i) { m_values[i] += y[i]; });
 
     return *this;
   }
 
-  PUGS_INLINE
-  Vector
-  operator+(const Vector& y) const
+  template <typename DataType2>
+  PUGS_INLINE Vector<std::remove_const_t<DataType>>
+  operator+(const Vector<DataType2>& y) const
   {
     Assert(this->size() == y.size());
-    Vector sum{y.size()};
+    Vector<std::remove_const_t<DataType>> sum{y.size()};
 
     parallel_for(
-      this->size(), PUGS_LAMBDA(const index_type& i) { sum.m_values[i] = m_values[i] + y.m_values[i]; });
+      this->size(), PUGS_LAMBDA(const index_type& i) { sum.m_values[i] = m_values[i] + y[i]; });
 
     return sum;
   }
 
-  PUGS_INLINE
-  Vector
-  operator-(const Vector& y) const
+  template <typename DataType2>
+  PUGS_INLINE Vector<std::remove_const_t<DataType>>
+  operator-(const Vector<DataType2>& y) const
   {
     Assert(this->size() == y.size());
-    Vector sum{y.size()};
+    Vector<std::remove_const_t<DataType>> sum{y.size()};
 
     parallel_for(
-      this->size(), PUGS_LAMBDA(const index_type& i) { sum.m_values[i] = m_values[i] - y.m_values[i]; });
+      this->size(), PUGS_LAMBDA(const index_type& i) { sum.m_values[i] = m_values[i] - y[i]; });
 
     return sum;
   }
@@ -153,6 +163,19 @@ class Vector   // LCOV_EXCL_LINE
   PUGS_INLINE
   Vector& operator=(Vector&&) = default;
 
+  template <typename DataType2>
+  Vector(const Vector<DataType2>& x)
+  {
+    // ensures that DataType is the same as source DataType2
+    static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+                  "Cannot assign Vector of different type");
+    // ensures that const is not lost through copy
+    static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
+                  "Cannot assign Vector of const to Vector of non-const");
+
+    this->operator=(x);
+  }
+
   Vector(const Vector&) = default;
 
   Vector(Vector&&) = default;

tests/CMakeLists.txt

@@ -7,6 +7,7 @@ add_executable (unit_tests
   test_ArrayUtils.cpp
   test_CRSMatrix.cpp
   test_ItemType.cpp
+  test_PCG.cpp
   test_PugsAssert.cpp
   test_RevisionInfo.cpp
   test_SparseMatrixDescriptor.cpp

tests/test_PCG.cpp

+#include <catch2/catch.hpp>
+
+#include <CRSMatrix.hpp>
+#include <PCG.hpp>
+
+TEST_CASE("PCG", "[algebra]")
+{
+  SECTION("no preconditionner")
+  {
+    SparseMatrixDescriptor S{5};
+    S(0, 0) = 2;
+    S(0, 1) = -1;
+
+    S(1, 0) = -1;
+    S(1, 1) = 2;
+    S(1, 2) = -1;
+
+    S(2, 1) = -1;
+    S(2, 2) = 2;
+    S(2, 3) = -1;
+
+    S(3, 2) = -1;
+    S(3, 3) = 2;
+    S(3, 4) = -1;
+
+    S(4, 3) = -1;
+    S(4, 4) = 2;
+
+    CRSMatrix A{S};
+
+    Vector<const double> x_exact = [] {
+      Vector<double> y{5};
+      y[0] = 1;
+      y[1] = 3;
+      y[2] = 2;
+      y[3] = 4;
+      y[4] = 5;
+      return y;
+    }();
+
+    Vector<double> b = A * x_exact;
+
+    Vector<double> x{5};
+    x = 0;
+
+    PCG{b, A, A, x, 10, 1e-12};
+    Vector error = x - x_exact;
+    REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
+  }
+}

tests/test_main.cpp

@@ -8,6 +8,9 @@ main(int argc, char* argv[])
 {
   Kokkos::initialize({4, -1, -1, true});
 
+  // Disable outputs from tested classes to the standard output
+  std::cout.setstate(std::ios::badbit);
+
   int result = Catch::Session().run(argc, argv);
 
   Kokkos::finalize();