diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index 224b09c2109fd41812344d498ad501670c776413..740343cd612f22911e187351b7929eedfbf7dcff 100644
--- a/src/algebra/CRSMatrix.hpp
+++ b/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)];
}
diff --git a/src/algebra/PCG.hpp b/src/algebra/PCG.hpp
index 2a979af2be8ad89c4ae3049f71b5d0b28af1a161..74c17ee40971b96c512c196d630f1edfadfa6ddc 100644
--- a/src/algebra/PCG.hpp
+++ b/src/algebra/PCG.hpp
@@ -1,6 +1,9 @@
#ifndef PCG_HPP
#define PCG_HPP
+#include <iomanip>
+#include <iostream>
+
struct PCG
{
template <typename VectorType, typename MatrixType, typename PreconditionerType>
diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp
index f8bbf531304d3952b7e0de052eaf055d56997869..43865dae1827cdb4875b3963127adafe78a8c6f3 100644
--- a/src/algebra/Vector.hpp
+++ b/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;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 2d0ec377d1b76e1de778ce019f85029fc9098204..a3becc0093004b637b552dc863a00990d9717f81 100644
--- a/tests/CMakeLists.txt
+++ b/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
diff --git a/tests/test_PCG.cpp b/tests/test_PCG.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..11636132cf38612cf662bc9cadfb1883553f966f
--- /dev/null
+++ b/tests/test_PCG.cpp
@@ -0,0 +1,50 @@
+#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)));
+ }
+}
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();