From 6f77ee95d6f09e1124cbc09101934d734efd0fa7 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 6 Aug 2019 18:47:43 +0200
Subject: [PATCH] Add BiCGStab tests and clean-up

---
 src/algebra/BiCGStab.hpp | 13 ++----
 src/main.cpp             |  2 +-
 tests/CMakeLists.txt     |  1 +
 tests/test_BiCGStab.cpp  | 93 ++++++++++++++++++++++++++++++++++++++++
 4 files changed, 99 insertions(+), 10 deletions(-)
 create mode 100644 tests/test_BiCGStab.cpp

diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
index 4faf83a05..fbcfa8103 100644
--- a/src/algebra/BiCGStab.hpp
+++ b/src/algebra/BiCGStab.hpp
@@ -1,18 +1,13 @@
 #ifndef BICG_STAB_HPP
 #define BICG_STAB_HPP
 
+#include <iomanip>
 #include <iostream>
 
-class BiCGStab
+struct 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)
+  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';
diff --git a/src/main.cpp b/src/main.cpp
index e2802e465..de641f550 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -98,7 +98,7 @@ main(int argc, char* argv[])
         Vector<double> u{x.size()};
         u = 0;
 
-        BiCGStab{Ax, A, A, u, 100, 1e-12};
+        BiCGStab{Ax, A, u, 100, 1e-12};
 
         std::cout << "== BiCGStab report ==\n";
         std::cout << "L2 Error = " << std::sqrt((x - u, x - u)) << '\n';
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index a3becc009..f4062e220 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -5,6 +5,7 @@ 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
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
new file mode 100644
index 000000000..8769200f6
--- /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)));
+  }
+}
-- 
GitLab