Skip to content
Snippets Groups Projects
Commit 6f77ee95 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Add BiCGStab tests and clean-up

parent dda63dda
No related branches found
No related tags found
1 merge request!26Feature/linear systems
#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';
......
......@@ -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';
......
......@@ -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
......
#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)));
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment