diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4faf83a054572821d07aa7d29fc08779c3ce8734
--- /dev/null
+++ b/src/algebra/BiCGStab.hpp
@@ -0,0 +1,82 @@
+#ifndef BICG_STAB_HPP
+#define BICG_STAB_HPP
+
+#include <iostream>
+
+class BiCGStab
+{
+ public:
+ template <typename VectorType, typename MatrixType, typename PreconditionerType>
+ BiCGStab(const VectorType& b,
+ const MatrixType& A,
+ [[maybe_unused]] const PreconditionerType& P,
+ 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/PCG.hpp b/src/algebra/PCG.hpp
index 357ab3a1d1275d6f79ceb82627e28f0003f06edd..2a979af2be8ad89c4ae3049f71b5d0b28af1a161 100644
--- a/src/algebra/PCG.hpp
+++ b/src/algebra/PCG.hpp
@@ -38,7 +38,7 @@ struct PCG
cg -= h;
- g = copy(cg); // C.computes(cg, g);
+ g = copy(cg); // TODO: precond: g = g/C
gcg = (g, cg);
@@ -55,8 +55,7 @@ struct PCG
double ro = gcg / hAh;
cg -= ro * b;
- // VectorType b2 = b;
- // C.computes(b2, b);
+ // TODO: precond: b <- b/C
x += ro * h;
g -= ro * b;
diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp
index 3a21d2a534517427fae8fd1afead63baf7bf40e3..d9537f0e488f39bedf8849c299f06221e62a7dc9 100644
--- a/src/algebra/Vector.hpp
+++ b/src/algebra/Vector.hpp
@@ -103,7 +103,7 @@ class Vector
PUGS_INLINE
Vector
- operator+(const Vector& y)
+ operator+(const Vector& y) const
{
Assert(this->size() == y.size());
Vector sum{y.size()};
@@ -116,7 +116,7 @@ class Vector
PUGS_INLINE
Vector
- operator-(const Vector& y)
+ operator-(const Vector& y) const
{
Assert(this->size() == y.size());
Vector sum{y.size()};
diff --git a/src/main.cpp b/src/main.cpp
index 3647754351565d38353ebfabc878aadfeed73292..c2f6b3a32d2d869611c753e69d501c2c29dbb225 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -24,6 +24,7 @@
#include <SynchronizerManager.hpp>
+#include <BiCGStab.hpp>
#include <CRSMatrix.hpp>
#include <PCG.hpp>
@@ -41,41 +42,70 @@ main(int argc, char* argv[])
Kokkos::initialize({4, -1, -1, true});
{
- size_t matrix_size = 10;
- FlexibleMatrix F(matrix_size);
-
- for (size_t i = 0; i < matrix_size; ++i) {
- F(i, i) = 4;
- if (i + 1 < matrix_size) {
- F(i, i + 1) = -1;
- F(i + 1, i) = -1;
+ { // PCG simple test
+ size_t matrix_size = 10;
+ FlexibleMatrix F(matrix_size);
+
+ for (size_t i = 0; i < matrix_size; ++i) {
+ F(i, i) = 4;
+ if (i + 1 < matrix_size) {
+ F(i, i + 1) = -1;
+ F(i + 1, i) = -1;
+ }
}
- }
- std::cout << F << '\n';
- CRSMatrix A{F};
+ CRSMatrix A{F};
+
+ Vector<double> x{10};
+ for (size_t i = 0; i < x.size(); ++i) {
+ x[i] = 2 * i + 1;
+ }
+
+ Vector<double> Ax = A * x;
+
+ Vector<double> u{x.size()};
+ u = 0;
- std::cout << A << '\n';
+ PCG{Ax, A, A, u, 100, 1e-12};
- Vector<double> x{10};
- for (size_t i = 0; i < x.size(); ++i) {
- x[i] = 2 * i + 1;
+ std::cout << "== PCG report ==\n";
+ std::cout << "L2 Error = " << std::sqrt((x - u, x - u)) << '\n';
+ std::cout << "L2 Norm = " << std::sqrt((x, x)) << '\n';
+ std::cout << "L2 RelEr = " << std::sqrt((x - u, x - u)) / std::sqrt((x, x)) << '\n';
+ std::cout << "================\n";
}
- std::cout << "x=\n" << x << '\n';
+ { // BiCGStab simple test
+ size_t matrix_size = 10;
+ FlexibleMatrix F(matrix_size);
- Vector<double> Ax = A * x;
- std::cout << "Ax=\n" << Ax << '\n';
+ for (size_t i = 0; i < matrix_size; ++i) {
+ F(i, i) = 4;
+ if (i + 1 < matrix_size) {
+ F(i, i + 1) = -1;
+ F(i + 1, i) = -0.3;
+ }
+ }
+
+ CRSMatrix A{F};
- Vector<double> u{x.size()};
- u = 0;
+ Vector<double> x{10};
+ for (size_t i = 0; i < x.size(); ++i) {
+ x[i] = 2 * i + 1;
+ }
+
+ Vector<double> Ax = A * x;
+ Vector<double> u{x.size()};
+ u = 0;
- PCG{Ax, A, A, u, static_cast<size_t>(-1), 1e-12};
- std::cout << "u=\n" << u << '\n';
+ BiCGStab{Ax, A, A, u, 100, 1e-12};
- std::cout << "L2 Error = " << std::sqrt((x - u, x - u)) << '\n';
- std::cout << "L2 Norm = " << std::sqrt((x, x)) << '\n';
- std::cout << "L2 RelEr = " << std::sqrt((x - u, x - u)) / std::sqrt((x, x)) << '\n';
+ std::cout << "== BiCGStab report ==\n";
+ std::cout << "L2 Error = " << std::sqrt((x - u, x - u)) << '\n';
+ std::cout << "L2 Norm = " << std::sqrt((x, x)) << '\n';
+ std::cout << "L2 RelEr = " << std::sqrt((x - u, x - u)) / std::sqrt((x, x)) << '\n';
+ std::cout << "=====================\n";
+ }
}
Kokkos::finalize();