diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp index 4faf83a054572821d07aa7d29fc08779c3ce8734..fbcfa8103e3e929ed2042b0d92b733c239850333 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 e2802e4654f6024ce88dd082c81113f1441982f2..de641f550fb44ab2a83055bd275ca34df4418020 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 a3becc0093004b637b552dc863a00990d9717f81..f4062e220422fb115235955b6c5d01da09f0f16d 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 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))); + } +}