From 0a24f6aeac3501e32b330be1edc150d76049c71f Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Thu, 1 Aug 2019 14:25:02 +0200
Subject: [PATCH] Add a first simple version of conjugate gradient
No preconditionner is used ATM
---
src/algebra/PCG.hpp | 94 ++++++++++++++++++++++++++++++++++++++++++
src/algebra/Vector.hpp | 19 ++++++++-
src/main.cpp | 12 +++++-
3 files changed, 122 insertions(+), 3 deletions(-)
create mode 100644 src/algebra/PCG.hpp
diff --git a/src/algebra/PCG.hpp b/src/algebra/PCG.hpp
new file mode 100644
index 000000000..357ab3a1d
--- /dev/null
+++ b/src/algebra/PCG.hpp
@@ -0,0 +1,94 @@
+#ifndef PCG_HPP
+#define PCG_HPP
+
+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); // C.computes(cg, g);
+
+ 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;
+
+ // VectorType b2 = b;
+ // C.computes(b2, b);
+
+ 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/Vector.hpp b/src/algebra/Vector.hpp
index bed915b78..5f6cff941 100644
--- a/src/algebra/Vector.hpp
+++ b/src/algebra/Vector.hpp
@@ -20,6 +20,16 @@ class Vector
static_assert(std::is_same_v<typename decltype(m_values)::index_type, index_type>);
public:
+ friend PUGS_INLINE Vector
+ copy(const Vector& x)
+ {
+ auto values = copy(x.m_values);
+
+ Vector x_copy{0};
+ x_copy.m_values = values;
+ return x_copy;
+ }
+
friend std::ostream&
operator<<(std::ostream& os, const Vector<DataType>& x)
{
@@ -29,6 +39,12 @@ class Vector
return os;
}
+ friend Vector operator*(const DataType& a, const Vector& x)
+ {
+ Vector y = copy(x);
+ return y *= a;
+ }
+
PUGS_INLINE
DataType
operator,(const Vector& y)
@@ -37,7 +53,7 @@ class Vector
// 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];
+ sum += m_values[i] * y.m_values[i];
}
return sum;
@@ -102,6 +118,7 @@ class Vector
operator=(const DataType& value) noexcept
{
m_values.fill(value);
+ return *this;
}
template <typename DataType2>
diff --git a/src/main.cpp b/src/main.cpp
index bcd9161bd..760164fcf 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -25,6 +25,7 @@
#include <SynchronizerManager.hpp>
#include <CRSMatrix.hpp>
+#include <PCG.hpp>
#include <iostream>
@@ -44,9 +45,10 @@ main(int argc, char* argv[])
FlexibleMatrix F(matrix_size);
for (size_t i = 0; i < matrix_size; ++i) {
- F(i, i) = 2;
+ F(i, i) = 4;
if (i + 1 < matrix_size) {
- F(i, i + 1) = 1;
+ F(i, i + 1) = -1;
+ F(i + 1, i) = -1;
}
}
std::cout << F << '\n';
@@ -64,6 +66,12 @@ main(int argc, char* argv[])
Vector<double> Ax = A * x;
std::cout << "Ax=\n" << Ax << '\n';
+
+ Vector<double> u{x.size()};
+ u = 0;
+
+ PCG{Ax, A, A, u, 9, 1e-4};
+ std::cout << "u=" << u << '\n';
}
Kokkos::finalize();
--
GitLab