#include <catch2/catch.hpp>

#include <algebra/BiCGStab.hpp>
#include <algebra/CRSMatrix.hpp>

// clazy:excludeall=non-pod-global-static

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<false>{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)));
  }
}