Merge branch 'develop' into feature/language

src/algebra/BiCGStab.hpp

+#ifndef BICG_STAB_HPP
+#define BICG_STAB_HPP
+
+#include <iomanip>
+#include <iostream>
+
+struct BiCGStab
+{
+  template <typename VectorType, typename MatrixType>
+  BiCGStab(const VectorType& b, const MatrixType& A, VectorType& x, const size_t max_iter, const double epsilon = 1e-6)
+  {
+    std::cout << "- bi-conjugate gradient stabilized\n";
+    std::cout << "  epsilon = " << epsilon << '\n';
+    std::cout << "  maximum number of iterations: " << max_iter << '\n';
+
+    VectorType r_k_1{b.size()};
+
+    r_k_1 = b - A * x;
+
+    double residu = std::sqrt((r_k_1, r_k_1));   // Norm(r_k_1);
+
+    if (residu != 0) {
+      double resid0 = residu;
+
+      VectorType rTilda_0 = copy(r_k_1);
+      VectorType p_k      = copy(r_k_1);
+
+      VectorType s_k{x.size()};
+
+      VectorType Ap_k{x.size()};
+      VectorType As_k{x.size()};
+
+      VectorType r_k{x.size()};
+
+      std::cout << "   initial residu: " << resid0 << '\n';
+      for (size_t i = 1; i <= max_iter; ++i) {
+        std::cout << "  - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0
+                  << "\tabsolute: " << residu << '\n';
+
+        Ap_k = A * p_k;
+
+        const double alpha_k = (r_k_1, rTilda_0) / (Ap_k, rTilda_0);
+
+        s_k  = r_k_1 - alpha_k * Ap_k;
+        As_k = A * s_k;
+
+        const double w_k = (As_k, s_k) / (As_k, As_k);
+
+        x += alpha_k * p_k + w_k * s_k;
+        r_k = s_k - w_k * As_k;
+
+        const double beta_k = (r_k, rTilda_0) / (r_k_1, rTilda_0) * (alpha_k / w_k);
+
+        p_k -= w_k * Ap_k;
+        p_k *= beta_k;
+        p_k += r_k;
+
+        if ((residu = std::sqrt((r_k, r_k))) / resid0 < epsilon) {
+          break;
+        }
+
+        r_k_1 = r_k;
+      }
+
+      if (residu / resid0 > epsilon) {
+        std::cout << "  bi_conjugate gradient stabilized: *NOT CONVERGED*\n";
+        std::cout << "  - epsilon:          " << epsilon << '\n';
+        std::cout << "  - relative residu : " << residu / resid0 << '\n';
+        std::cout << "  - absolute residu : " << residu << '\n';
+      }
+    }
+
+    std::cout << "- bi-conjugate gradient stabilized: done\n";
+  }
+};
+
+#endif   // BICG_STAB_HPP

src/algebra/CRSMatrix.hpp

+#ifndef CRS_MATRIX_HPP
+#define CRS_MATRIX_HPP
+
+#include <Array.hpp>
+#include <Kokkos_StaticCrsGraph.hpp>
+#include <PugsAssert.hpp>
+
+#include <SparseMatrixDescriptor.hpp>
+
+#include <Vector.hpp>
+
+#include <iostream>
+
+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<const DataType> m_values;
+
+ public:
+  PUGS_INLINE
+  size_t
+  numberOfRows() const
+  {
+    return m_host_matrix.numRows();
+  }
+
+  template <typename DataType2>
+  Vector<MutableDataType> operator*(const Vector<DataType2>& x) const
+  {
+    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);
+        MutableDataType sum{0};
+        for (IndexType j = row_begin; j < row_end; ++j) {
+          sum += m_values[j] * x[m_host_matrix.entries(j)];
+        }
+        Ax[i_row] = sum;
+      });
+
+    return Ax;
+  }
+
+  CRSMatrix(const SparseMatrixDescriptor<DataType, IndexType>& M)
+  {
+    m_host_matrix = Kokkos::create_staticcrsgraph<HostMatrix>("connectivity_matrix", M.graphVector());
+    m_values      = M.valueArray();
+  }
+  ~CRSMatrix() = default;
+};
+
+#endif   // CRS_MATRIX_HPP

src/algebra/PCG.hpp

+#ifndef PCG_HPP
+#define PCG_HPP
+
+#include <iomanip>
+#include <iostream>
+
+struct PCG
+{
+  template <typename VectorType, typename MatrixType, typename PreconditionerType>
+  PCG(const VectorType& f,
+      const MatrixType& A,
+      [[maybe_unused]] const PreconditionerType& C,
+      VectorType& x,
+      const size_t maxiter,
+      const double epsilon = 1e-6)
+  {
+    std::cout << "- conjugate gradient\n";
+    std::cout << "  epsilon = " << epsilon << '\n';
+    std::cout << "  maximum number of iterations: " << maxiter << '\n';
+
+    VectorType h{f.size()};
+    VectorType b = copy(f);
+
+    if (true) {
+      h = A * x;
+      h -= f;
+      std::cout << "- initial *real* residu :   " << (h, h) << '\n';
+    }
+
+    VectorType g{b.size()};
+    VectorType cg = copy(b);
+
+    double gcg  = 0;
+    double gcg0 = 1;
+
+    double relativeEpsilon = epsilon;
+
+    for (size_t i = 1; i <= maxiter; ++i) {
+      if (i == 1) {
+        h = A * x;
+
+        cg -= h;
+
+        g = copy(cg);   // TODO: precond: g = g/C
+
+        gcg = (g, cg);
+
+        h = copy(g);
+      }
+
+      b = A * h;
+
+      double hAh = (h, b);
+
+      if (hAh == 0) {
+        hAh = 1.;
+      }
+      double ro = gcg / hAh;
+      cg -= ro * b;
+
+      // TODO: precond: b <- b/C
+
+      x += ro * h;
+      g -= ro * b;
+
+      double gamma = gcg;
+      gcg          = (g, cg);
+
+      if ((i == 1) && (gcg != 0)) {
+        relativeEpsilon = epsilon * gcg;
+        gcg0            = gcg;
+        std::cout << "  initial residu: " << gcg << '\n';
+      }
+      std::cout << "  - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
+      std::cout << "\tabsolute: " << gcg << '\n';
+
+      if (gcg < relativeEpsilon) {
+        break;
+      }
+
+      gamma = gcg / gamma;
+
+      h *= gamma;
+      h += g;
+    }
+
+    if (gcg > relativeEpsilon) {
+      std::cout << "  conjugate gradient: *NOT CONVERGED*\n";
+      std::cout << "  - epsilon:          " << epsilon << '\n';
+      std::cout << "  - relative residu : " << gcg / gcg0 << '\n';
+      std::cout << "  - absolute residu : " << gcg << '\n';
+    }
+  }
+};
+
+#endif   // PCG_HPP

src/algebra/SparseMatrixDescriptor.hpp

+#ifndef SPARSE_MATRIX_DESCRIPTOR_HPP
+#define SPARSE_MATRIX_DESCRIPTOR_HPP
+
+#include <Array.hpp>
+
+#include <map>
+#include <type_traits>
+
+template <typename DataType = double, typename IndexType = size_t>
+class SparseMatrixDescriptor
+{
+  static_assert(std::is_integral_v<IndexType>, "Index type must be an integral type");
+  static_assert(std::is_unsigned_v<IndexType>, "Index type must be unsigned");
+
+ public:
+  using data_type  = DataType;
+  using index_type = IndexType;
+
+  class SparseRowDescriptor
+  {
+    friend class SparseMatrixDescriptor;
+    std::map<IndexType, DataType> m_id_value_map;
+
+   public:
+    IndexType
+    numberOfValues() const noexcept
+    {
+      return m_id_value_map.size();
+    }
+
+    const DataType& operator[](const IndexType& j) const
+    {
+      auto i_value = m_id_value_map.find(j);
+      Assert(i_value != m_id_value_map.end());
+      return i_value->second;
+    }
+
+    DataType& operator[](const IndexType& j)
+    {
+      auto i_value = m_id_value_map.find(j);
+      if (i_value != m_id_value_map.end()) {
+        return i_value->second;
+      } else {
+        auto [i_inserted, success] = m_id_value_map.insert(std::make_pair(j, 0));
+        Assert(success);
+        return i_inserted->second;
+      }
+    }
+
+    friend std::ostream&
+    operator<<(std::ostream& os, const SparseRowDescriptor& row)
+    {
+      for (auto [j, value] : row.m_id_value_map) {
+        os << ' ' << j << ':' << value;
+      }
+      return os;
+    }
+
+    SparseRowDescriptor() = default;
+  };
+
+ private:
+  Array<SparseRowDescriptor> m_row_array;
+
+ public:
+  PUGS_INLINE
+  size_t
+  numberOfRows() const noexcept
+  {
+    return m_row_array.size();
+  }
+
+  SparseRowDescriptor&
+  row(const IndexType i)
+  {
+    return m_row_array[i];
+  }
+
+  const SparseRowDescriptor&
+  row(const IndexType i) const
+  {
+    return m_row_array[i];
+  }
+
+  DataType&
+  operator()(const IndexType& i, const IndexType& j)
+  {
+    return m_row_array[i][j];
+  }
+
+  const DataType&
+  operator()(const IndexType& i, const IndexType& j) const
+  {
+    const auto& r = m_row_array[i];   // split to ensure const-ness of call
+    return r[j];
+  }
+
+  friend std::ostream&
+  operator<<(std::ostream& os, const SparseMatrixDescriptor& M)
+  {
+    for (IndexType i = 0; i < M.m_row_array.size(); ++i) {
+      os << i << " |" << M.m_row_array[i] << '\n';
+    }
+    return os;
+  }
+
+  auto
+  graphVector() const
+  {
+    std::vector<std::vector<IndexType>> graph_vector(m_row_array.size());
+    for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
+      const SparseRowDescriptor& row = m_row_array[i_row];
+      for (auto [id, value] : row.m_id_value_map) {
+        graph_vector[i_row].push_back(id);
+      }
+    }
+    return graph_vector;
+  }
+
+  Array<DataType>
+  valueArray() const
+  {
+    IndexType size = 0;
+    for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
+      size += m_row_array[i_row].m_id_value_map.size();
+    }
+
+    Array<DataType> values(size);
+
+    IndexType cpt = 0;
+    for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
+      const SparseRowDescriptor& row = m_row_array[i_row];
+      for (auto [id, value] : row.m_id_value_map) {
+        values[cpt++] = value;
+      }
+    }
+    return values;
+  }
+
+  SparseMatrixDescriptor(size_t nb_row) : m_row_array{nb_row} {}
+
+  ~SparseMatrixDescriptor() = default;
+};
+
+#endif   // SPARSE_MATRIX_DESCRIPTOR_HPP

src/algebra/Vector.hpp

+#ifndef VECTOR_HPP
+#define VECTOR_HPP
+
+#include <PugsMacros.hpp>
+#include <PugsUtils.hpp>
+
+#include <PugsAssert.hpp>
+
+#include <Array.hpp>
+
+template <typename DataType>
+class Vector   // LCOV_EXCL_LINE
+{
+ public:
+  using data_type  = DataType;
+  using index_type = size_t;
+
+ private:
+  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 Vector<std::remove_const_t<DataType>>
+  copy(const Vector& x)
+  {
+    auto values = copy(x.m_values);
+
+    Vector<std::remove_const_t<DataType>> x_copy{0};
+    x_copy.m_values = values;
+    return x_copy;
+  }
+
+  friend Vector operator*(const DataType& a, const Vector& x)
+  {
+    Vector<std::remove_const_t<DataType>> y = copy(x);
+    return y *= a;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE
+  auto
+  operator,(const Vector<DataType2>& y)
+  {
+    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) {
+      sum += m_values[i] * y.m_values[i];
+    }
+
+    return sum;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE Vector&
+  operator/=(const DataType2& a)
+  {
+    const auto inv_a = 1. / a;
+    return (*this) *= inv_a;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE Vector&
+  operator*=(const DataType2& a)
+  {
+    parallel_for(
+      this->size(), PUGS_LAMBDA(const index_type& i) { m_values[i] *= a; });
+    return *this;
+  }
+
+  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[i]; });
+
+    return *this;
+  }
+
+  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[i]; });
+
+    return *this;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE Vector<std::remove_const_t<DataType>>
+  operator+(const Vector<DataType2>& y) const
+  {
+    Assert(this->size() == 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[i]; });
+
+    return sum;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE Vector<std::remove_const_t<DataType>>
+  operator-(const Vector<DataType2>& y) const
+  {
+    Assert(this->size() == 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[i]; });
+
+    return sum;
+  }
+
+  PUGS_INLINE
+  DataType& operator[](const index_type& i) const noexcept(NO_ASSERT)
+  {
+    return m_values[i];
+  }
+
+  PUGS_INLINE
+  size_t
+  size() const noexcept
+  {
+    return m_values.size();
+  }
+
+  PUGS_INLINE Vector&
+  operator=(const DataType& value) noexcept
+  {
+    m_values.fill(value);
+    return *this;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE Vector&
+  operator=(const Vector<DataType2>& vector) noexcept
+  {
+    m_values = vector.m_values;
+    return *this;
+  }
+
+  PUGS_INLINE
+  Vector& operator=(const Vector&) = default;
+
+  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;
+
+  Vector(const size_t& size) : m_values{size} {}
+  ~Vector() = default;
+};
+
+#endif   // VECTOR_HPP

tests/CMakeLists.txt

@@ -18,12 +18,17 @@ add_executable (unit_tests
   test_ASTPrinter.cpp
   test_ASTSymbolTableBuilder.cpp
   test_ASTSymbolInitializationChecker.cpp
+  test_BiCGStab.cpp
+  test_CRSMatrix.cpp
   test_ItemType.cpp
+  test_PCG.cpp
   test_PugsAssert.cpp
   test_RevisionInfo.cpp
+  test_SparseMatrixDescriptor.cpp
   test_SymbolTable.cpp
   test_TinyMatrix.cpp
   test_TinyVector.cpp
+  test_Vector.cpp
   )
 
 add_executable (mpi_unit_tests

tests/test_BiCGStab.cpp

+#include <catch2/catch.hpp>
+
+#include <BiCGStab.hpp>
+#include <CRSMatrix.hpp>
+
+TEST_CASE("BiCGStab", "[algebra]")
+{
+  SECTION("no preconditionner")
+  {
+    SparseMatrixDescriptor S{5};
+    S(0, 0) = 2;
+    S(0, 1) = -1;
+
+    S(1, 0) = -0.2;
+    S(1, 1) = 2;
+    S(1, 2) = -1;
+
+    S(2, 1) = -1;
+    S(2, 2) = 4;
+    S(2, 3) = -2;
+
+    S(3, 2) = -1;
+    S(3, 3) = 2;
+    S(3, 4) = -0.1;
+
+    S(4, 3) = 1;
+    S(4, 4) = 3;
+
+    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;
+
+    BiCGStab{b, A, x, 10, 1e-12};
+    Vector error = x - x_exact;
+    REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
+  }
+
+  SECTION("no preconditionner non-converged")
+  {
+    SparseMatrixDescriptor S{5};
+    S(0, 0) = 2;
+    S(0, 1) = -1;
+
+    S(1, 0) = -0.2;
+    S(1, 1) = 2;
+    S(1, 2) = -1;
+
+    S(2, 1) = -1;
+    S(2, 2) = 4;
+    S(2, 3) = -2;
+
+    S(3, 2) = -1;
+    S(3, 3) = 2;
+    S(3, 4) = -0.1;
+
+    S(4, 3) = 1;
+    S(4, 4) = 3;
+
+    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;
+
+    BiCGStab{b, A, x, 1, 1e-12};
+    Vector error = x - x_exact;
+    REQUIRE(std::sqrt((error, error)) > 1E-5 * std::sqrt((x, x)));
+  }
+}

tests/test_CRSMatrix.cpp

+#include <catch2/catch.hpp>
+
+#include <CRSMatrix.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class SparseMatrixDescriptor<int, uint8_t>;
+
+TEST_CASE("CRSMatrix", "[algebra]")
+{
+  SECTION("matrix size")
+  {
+    SparseMatrixDescriptor<int, uint8_t> S{5};
+    S(0, 2) = 3;
+    S(1, 2) = 11;
+    S(1, 1) = 1;
+    S(3, 0) = 5;
+    S(3, 1) = -3;
+    S(4, 1) = 1;
+    S(2, 2) = 4;
+    S(4, 3) = 2;
+    S(4, 4) = -2;
+
+    CRSMatrix<int, uint8_t> A{S};
+
+    REQUIRE(A.numberOfRows() == S.numberOfRows());
+  }
+
+  SECTION("matrix vector product (simple)")
+  {
+    SparseMatrixDescriptor<int, uint8_t> S{5};
+    S(0, 2) = 3;
+    S(1, 2) = 11;
+    S(1, 1) = 1;
+    S(3, 0) = 5;
+    S(3, 1) = -3;
+    S(4, 1) = 1;
+    S(2, 2) = 4;
+    S(4, 3) = 2;
+    S(4, 4) = -2;
+
+    CRSMatrix<int, uint8_t> A{S};
+
+    Vector<int> x{A.numberOfRows()};
+    x    = 0;
+    x[0] = 1;
+
+    Vector<int> y = A * x;
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 0);
+    REQUIRE(y[2] == 0);
+    REQUIRE(y[3] == 5);
+    REQUIRE(y[4] == 0);
+
+    x    = 0;
+    x[1] = 2;
+
+    y = A * x;
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 2);
+    REQUIRE(y[2] == 0);
+    REQUIRE(y[3] == -6);
+    REQUIRE(y[4] == 2);
+
+    x    = 0;
+    x[2] = -1;
+
+    y = A * x;
+    REQUIRE(y[0] == -3);
+    REQUIRE(y[1] == -11);
+    REQUIRE(y[2] == -4);
+    REQUIRE(y[3] == 0);
+    REQUIRE(y[4] == 0);
+
+    x    = 0;
+    x[3] = 3;
+
+    y = A * x;
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 0);
+    REQUIRE(y[2] == 0);
+    REQUIRE(y[3] == 0);
+    REQUIRE(y[4] == 6);
+
+    x    = 0;
+    x[4] = 1;
+
+    y = A * x;
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 0);
+    REQUIRE(y[2] == 0);
+    REQUIRE(y[3] == 0);
+    REQUIRE(y[4] == -2);
+  }
+
+  SECTION("matrix vector product (complet)")
+  {
+    SparseMatrixDescriptor<int, uint8_t> S{4};
+    S(0, 0) = 1;
+    S(0, 1) = 2;
+    S(0, 2) = 3;
+    S(0, 3) = 4;
+    S(1, 0) = 5;
+    S(1, 1) = 6;
+    S(1, 2) = 7;
+    S(1, 3) = 8;
+    S(2, 0) = 9;
+    S(2, 1) = 10;
+    S(2, 2) = 11;
+    S(2, 3) = 12;
+    S(3, 0) = 13;
+    S(3, 1) = 14;
+    S(3, 2) = 15;
+    S(3, 3) = 16;
+
+    CRSMatrix<int, uint8_t> A{S};
+
+    Vector<int> x{A.numberOfRows()};
+    x[0] = 1;
+    x[1] = 2;
+    x[2] = 3;
+    x[3] = 4;
+
+    Vector<int> y = A * x;
+    REQUIRE(y[0] == 30);
+    REQUIRE(y[1] == 70);
+    REQUIRE(y[2] == 110);
+    REQUIRE(y[3] == 150);
+  }
+
+#ifndef NDEBUG
+  SECTION("incompatible runtime matrix/vector product")
+  {
+    CRSMatrix<int, uint8_t> A{SparseMatrixDescriptor<int, uint8_t>{4}};
+    Vector<int> x{2};
+    REQUIRE_THROWS_AS(A * x, AssertError);
+  }
+#endif   // NDEBUG
+}

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)));
+  }
+
+  SECTION("singular matrix with RHS in its image")
+  {
+    SparseMatrixDescriptor<double, uint8_t> S{5};
+
+    CRSMatrix<double, uint8_t> A{S};
+
+    Vector<double> b{5};
+    b = 0;
+
+    Vector<double> x{5};
+    x = 0;
+
+    PCG{b, A, A, x, 10, 1e-12};
+    REQUIRE(std::sqrt((x, x)) == 0);
+  }
+
+  SECTION("no preconditionner non-converged")
+  {
+    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, 1, 1e-12};
+    Vector error = x - x_exact;
+    REQUIRE(std::sqrt((error, error)) > 1E-10 * std::sqrt((x, x)));
+  }
+}

tests/test_SparseMatrixDescriptor.cpp

+#include <catch2/catch.hpp>
+
+#include <PugsAssert.hpp>
+#include <SparseMatrixDescriptor.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class SparseMatrixDescriptor<int, uint8_t>;
+
+TEST_CASE("SparseMatrixDescriptor", "[algebra]")
+{
+  SECTION("SparseRowDescriptor subclass")
+  {
+    SECTION("number of values")
+    {
+      SparseMatrixDescriptor<int, uint8_t>::SparseRowDescriptor r;
+      REQUIRE(r.numberOfValues() == 0);
+
+      r[0] = 3;
+      r[0] += 2;
+
+      REQUIRE(r.numberOfValues() == 1);
+
+      r[1] = -2;
+
+      REQUIRE(r.numberOfValues() == 2);
+    }
+
+    SECTION("values access")
+    {
+      SparseMatrixDescriptor<int, uint8_t>::SparseRowDescriptor r;
+
+      r[0] = 3;
+      r[0] += 2;
+
+      r[3] = 8;
+
+      REQUIRE(r[0] == 5);
+      REQUIRE(r[3] == 8);
+
+      const auto& const_r = r;
+      REQUIRE(const_r[0] == 5);
+      REQUIRE(const_r[3] == 8);
+
+#ifndef NDEBUG
+      REQUIRE_THROWS_AS(const_r[2], AssertError);
+#endif   // NDEBUG
+    }
+  }
+
+  SECTION("number of columns")
+  {
+    SparseMatrixDescriptor<int, uint8_t> S{5};
+    REQUIRE(S.numberOfRows() == 5);
+  }
+
+  SECTION("location operators")
+  {
+    SparseMatrixDescriptor<int, uint8_t> S{5};
+    S(0, 2) = 3;
+    S(1, 1) = 1;
+    S(1, 2) = 11;
+    S(0, 2) += 2;
+    S(2, 2) = 4;
+    S(3, 1) = 5;
+    S(4, 1) = 1;
+    S(4, 4) = -2;
+
+    REQUIRE(S.row(0).numberOfValues() == 1);
+    REQUIRE(S.row(1).numberOfValues() == 2);
+    REQUIRE(S.row(2).numberOfValues() == 1);
+    REQUIRE(S.row(3).numberOfValues() == 1);
+    REQUIRE(S.row(4).numberOfValues() == 2);
+
+    const auto& const_S = S;
+
+    REQUIRE(const_S.row(0).numberOfValues() == 1);
+    REQUIRE(const_S.row(1).numberOfValues() == 2);
+    REQUIRE(const_S.row(2).numberOfValues() == 1);
+    REQUIRE(const_S.row(3).numberOfValues() == 1);
+    REQUIRE(const_S.row(4).numberOfValues() == 2);
+
+#ifndef NDEBUG
+    REQUIRE_THROWS_AS(S.row(5), AssertError);
+    REQUIRE_THROWS_AS(const_S.row(5), AssertError);
+#endif   // NDEBUG
+
+    REQUIRE(S(0, 2) == 5);
+    REQUIRE(S(1, 1) == 1);
+    REQUIRE(S(1, 2) == 11);
+    REQUIRE(S(2, 2) == 4);
+    REQUIRE(S(3, 1) == 5);
+    REQUIRE(S(4, 1) == 1);
+    REQUIRE(S(4, 4) == -2);
+
+    REQUIRE(const_S(0, 2) == 5);
+    REQUIRE(const_S(1, 1) == 1);
+    REQUIRE(const_S(1, 2) == 11);
+    REQUIRE(const_S(2, 2) == 4);
+    REQUIRE(const_S(3, 1) == 5);
+    REQUIRE(const_S(4, 1) == 1);
+    REQUIRE(const_S(4, 4) == -2);
+
+#ifndef NDEBUG
+    REQUIRE_THROWS_AS(S(5, 0), AssertError);
+    REQUIRE_THROWS_AS(const_S(6, 1), AssertError);
+    REQUIRE_THROWS_AS(const_S(0, 1), AssertError);
+#endif   // NDEBUG
+  }
+
+  SECTION("vector-graph")
+  {
+    SparseMatrixDescriptor<int, uint8_t> S{5};
+    S(0, 2) = 3;
+    S(1, 2) = 11;
+    S(1, 1) = 1;
+    S(0, 2) += 2;
+    S(3, 3) = 5;
+    S(3, 1) = 5;
+    S(4, 1) = 1;
+    S(2, 2) = 4;
+    S(4, 4) = -2;
+
+    const auto graph = S.graphVector();
+
+    REQUIRE(graph.size() == S.numberOfRows());
+    REQUIRE(graph[0].size() == 1);
+    REQUIRE(graph[1].size() == 2);
+    REQUIRE(graph[2].size() == 1);
+    REQUIRE(graph[3].size() == 2);
+    REQUIRE(graph[4].size() == 2);
+
+    REQUIRE(graph[0][0] == 2);
+    REQUIRE(graph[1][0] == 1);
+    REQUIRE(graph[1][1] == 2);
+    REQUIRE(graph[2][0] == 2);
+    REQUIRE(graph[3][0] == 1);
+    REQUIRE(graph[3][1] == 3);
+    REQUIRE(graph[4][0] == 1);
+    REQUIRE(graph[4][1] == 4);
+  }
+
+  SECTION("value array")
+  {
+    SparseMatrixDescriptor<int, uint8_t> S{5};
+    S(0, 2) = 3;
+    S(1, 2) = 11;
+    S(1, 1) = 1;
+    S(0, 2) += 2;
+    S(3, 3) = 5;
+    S(3, 1) = -3;
+    S(4, 1) = 1;
+    S(2, 2) = 4;
+    S(4, 4) = -2;
+
+    const auto value_array = S.valueArray();
+
+    REQUIRE(value_array.size() == 8);
+
+    REQUIRE(value_array[0] == 5);
+    REQUIRE(value_array[1] == 1);
+    REQUIRE(value_array[2] == 11);
+    REQUIRE(value_array[3] == 4);
+    REQUIRE(value_array[4] == -3);
+    REQUIRE(value_array[5] == 5);
+    REQUIRE(value_array[6] == 1);
+    REQUIRE(value_array[7] == -2);
+  }
+}

tests/test_Vector.cpp

+#include <catch2/catch.hpp>
+
+#include <PugsAssert.hpp>
+#include <Vector.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class Vector<int>;
+
+TEST_CASE("Vector", "[algebra]")
+{
+  SECTION("size")
+  {
+    Vector<int> x{5};
+    REQUIRE(x.size() == 5);
+  }
+
+  SECTION("write access")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    REQUIRE(x[0] == 0);
+    REQUIRE(x[1] == 1);
+    REQUIRE(x[2] == 2);
+    REQUIRE(x[3] == 3);
+    REQUIRE(x[4] == 4);
+  }
+
+  SECTION("fill")
+  {
+    Vector<int> x{5};
+    x = 2;
+
+    REQUIRE(x[0] == 2);
+    REQUIRE(x[1] == 2);
+    REQUIRE(x[2] == 2);
+    REQUIRE(x[3] == 2);
+    REQUIRE(x[4] == 2);
+  }
+
+  SECTION("copy constructor (shallow)")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    const Vector<int> y = x;
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 1);
+    REQUIRE(y[2] == 2);
+    REQUIRE(y[3] == 3);
+    REQUIRE(y[4] == 4);
+  }
+
+  SECTION("copy (deep)")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    const Vector<int> y = copy(x);
+
+    x[0] = 14;
+    x[1] = 13;
+    x[2] = 12;
+    x[3] = 11;
+    x[4] = 10;
+
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 1);
+    REQUIRE(y[2] == 2);
+    REQUIRE(y[3] == 3);
+    REQUIRE(y[4] == 4);
+
+    REQUIRE(x[0] == 14);
+    REQUIRE(x[1] == 13);
+    REQUIRE(x[2] == 12);
+    REQUIRE(x[3] == 11);
+    REQUIRE(x[4] == 10);
+  }
+
+  SECTION("self scalar multiplication")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    x *= 2;
+
+    REQUIRE(x[0] == 0);
+    REQUIRE(x[1] == 2);
+    REQUIRE(x[2] == 4);
+    REQUIRE(x[3] == 6);
+    REQUIRE(x[4] == 8);
+  }
+
+  SECTION("left scalar multiplication")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    const Vector<int> y = 2 * x;
+
+    REQUIRE(y[0] == 0);
+    REQUIRE(y[1] == 2);
+    REQUIRE(y[2] == 4);
+    REQUIRE(y[3] == 6);
+    REQUIRE(y[4] == 8);
+  }
+
+  SECTION("dot product")
+  {
+    Vector<int> x{5};
+    x[0] = 0;
+    x[1] = 1;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 4;
+
+    Vector<int> y{5};
+    y[0] = 7;
+    y[1] = 3;
+    y[2] = 4;
+    y[3] = 2;
+    y[4] = 8;
+
+    const int dot = (x, y);
+    REQUIRE(dot == 49);
+  }
+
+  SECTION("self scalar division")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    x /= 2;
+
+    REQUIRE(x[0] == 1);
+    REQUIRE(x[1] == 1);
+    REQUIRE(x[2] == 2);
+    REQUIRE(x[3] == 3);
+    REQUIRE(x[4] == 4);
+  }
+
+  SECTION("self minus")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    Vector<int> y{5};
+    y[0] = 3;
+    y[1] = 8;
+    y[2] = 6;
+    y[3] = 2;
+    y[4] = 4;
+
+    x -= y;
+
+    REQUIRE(x[0] == -1);
+    REQUIRE(x[1] == -5);
+    REQUIRE(x[2] == -1);
+    REQUIRE(x[3] == 5);
+    REQUIRE(x[4] == 4);
+  }
+
+  SECTION("self sum")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    Vector<int> y{5};
+    y[0] = 3;
+    y[1] = 8;
+    y[2] = 6;
+    y[3] = 2;
+    y[4] = 4;
+
+    x += y;
+
+    REQUIRE(x[0] == 5);
+    REQUIRE(x[1] == 11);
+    REQUIRE(x[2] == 11);
+    REQUIRE(x[3] == 9);
+    REQUIRE(x[4] == 12);
+  }
+
+  SECTION("sum")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    Vector<int> y{5};
+    y[0] = 3;
+    y[1] = 8;
+    y[2] = 6;
+    y[3] = 2;
+    y[4] = 4;
+
+    Vector z = x + y;
+
+    REQUIRE(z[0] == 5);
+    REQUIRE(z[1] == 11);
+    REQUIRE(z[2] == 11);
+    REQUIRE(z[3] == 9);
+    REQUIRE(z[4] == 12);
+  }
+
+  SECTION("difference")
+  {
+    Vector<int> x{5};
+    x[0] = 2;
+    x[1] = 3;
+    x[2] = 5;
+    x[3] = 7;
+    x[4] = 8;
+
+    Vector<int> y{5};
+    y[0] = 3;
+    y[1] = 8;
+    y[2] = 6;
+    y[3] = 2;
+    y[4] = 4;
+
+    Vector z = x - y;
+
+    REQUIRE(z[0] == -1);
+    REQUIRE(z[1] == -5);
+    REQUIRE(z[2] == -1);
+    REQUIRE(z[3] == 5);
+    REQUIRE(z[4] == 4);
+  }
+}