diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp new file mode 100644 index 0000000000000000000000000000000000000000..fbcfa8103e3e929ed2042b0d92b733c239850333 --- /dev/null +++ b/src/algebra/BiCGStab.hpp @@ -0,0 +1,77 @@ +#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 diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp new file mode 100644 index 0000000000000000000000000000000000000000..740343cd612f22911e187351b7929eedfbf7dcff --- /dev/null +++ b/src/algebra/CRSMatrix.hpp @@ -0,0 +1,66 @@ +#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 diff --git a/src/algebra/PCG.hpp b/src/algebra/PCG.hpp new file mode 100644 index 0000000000000000000000000000000000000000..74c17ee40971b96c512c196d630f1edfadfa6ddc --- /dev/null +++ b/src/algebra/PCG.hpp @@ -0,0 +1,96 @@ +#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 diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp new file mode 100644 index 0000000000000000000000000000000000000000..badfb96305c72dbe5777c689a0861dfed299ff51 --- /dev/null +++ b/src/algebra/SparseMatrixDescriptor.hpp @@ -0,0 +1,145 @@ +#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 diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp new file mode 100644 index 0000000000000000000000000000000000000000..43865dae1827cdb4875b3963127adafe78a8c6f3 --- /dev/null +++ b/src/algebra/Vector.hpp @@ -0,0 +1,187 @@ +#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 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index dad75a12071765ed09cba7816a4a8f71469c18fc..f4062e220422fb115235955b6c5d01da09f0f16d 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -5,11 +5,16 @@ add_executable (unit_tests test_main.cpp test_Array.cpp test_ArrayUtils.cpp + test_BiCGStab.cpp + test_CRSMatrix.cpp test_ItemType.cpp + test_PCG.cpp test_PugsAssert.cpp test_RevisionInfo.cpp + test_SparseMatrixDescriptor.cpp test_TinyMatrix.cpp test_TinyVector.cpp + test_Vector.cpp ) add_executable (mpi_unit_tests diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8769200f65f7b0f393a1fde3838a68307e798730 --- /dev/null +++ b/tests/test_BiCGStab.cpp @@ -0,0 +1,93 @@ +#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))); + } +} diff --git a/tests/test_CRSMatrix.cpp b/tests/test_CRSMatrix.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4d555b2810467edc761a1758299887c2d851a1b7 --- /dev/null +++ b/tests/test_CRSMatrix.cpp @@ -0,0 +1,138 @@ +#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 +} diff --git a/tests/test_PCG.cpp b/tests/test_PCG.cpp new file mode 100644 index 0000000000000000000000000000000000000000..97bbe1ce25076d05f24b3e250d779d799012beab --- /dev/null +++ b/tests/test_PCG.cpp @@ -0,0 +1,109 @@ +#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))); + } +} diff --git a/tests/test_SparseMatrixDescriptor.cpp b/tests/test_SparseMatrixDescriptor.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c34352f61cfdf81b7e9c2b3e3753562b29a94de0 --- /dev/null +++ b/tests/test_SparseMatrixDescriptor.cpp @@ -0,0 +1,168 @@ +#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); + } +} diff --git a/tests/test_Vector.cpp b/tests/test_Vector.cpp new file mode 100644 index 0000000000000000000000000000000000000000..15e02d0eb2c7f9c0ad28af2cdb365ed790bd563f --- /dev/null +++ b/tests/test_Vector.cpp @@ -0,0 +1,265 @@ +#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); + } +} diff --git a/tests/test_main.cpp b/tests/test_main.cpp index e0ed7b5fe279e1b1b4b21d10a4145701ab619fc8..c7355e990f3322d22de1ab0398d7b5ed0fab3c06 100644 --- a/tests/test_main.cpp +++ b/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();