test_CG.cpp

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>

#include <algebra/CG.hpp>
#include <algebra/CRSMatrix.hpp>
#include <algebra/CRSMatrixDescriptor.hpp>
#include <algebra/Vector.hpp>

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

TEST_CASE("CG", "[algebra]")
{
  SECTION("no preconditionner")
  {
    Array<int> non_zeros{5};
    non_zeros.fill(3);
    non_zeros[0] = 2;
    non_zeros[4] = 2;
    CRSMatrixDescriptor S{5, 5, non_zeros};
    S(0, 0) = 2;
    S(0, 1) = -1;

    S(1, 0) = -1;
    S(1, 1) = 2;
    S(1, 2) = -1;

    S(2, 1) = -1;
    S(2, 2) = 2;
    S(2, 3) = -1;

    S(3, 2) = -1;
    S(3, 3) = 2;
    S(3, 4) = -1;

    S(4, 3) = -1;
    S(4, 4) = 2;

    CRSMatrix A{S.getCRSMatrix()};

    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 = zero;

    CG{A, x, b, 1e-12, 10, true};
    Vector error = x - x_exact;
    REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x, x)));
  }

  SECTION("singular matrix with RHS in its image")
  {
    Array<int> non_zeros(5);
    non_zeros.fill(0);
    CRSMatrixDescriptor<double> S{5, 5, non_zeros};

    std::ostringstream Sout;
    Sout << S;
    REQUIRE(Sout.str() == R"(0| 0:0
1| 1:0
2| 2:0
3| 3:0
4| 4:0
)");

    CRSMatrix<double> A{S.getCRSMatrix()};

    std::ostringstream Aout;
    Aout << A;
    REQUIRE(Sout.str() == Aout.str());

    Vector<double> b{5};
    b = zero;

    Vector<double> x{5};
    x = zero;

    CG{A, x, b, 1e-12, 10};
    REQUIRE(std::sqrt(dot(x, x)) == 0);
  }

  SECTION("no preconditionner non-converged")
  {
    Array<int> non_zeros(5);
    non_zeros.fill(3);
    non_zeros[0] = 2;
    non_zeros[4] = 2;
    CRSMatrixDescriptor<double> S{5, 5, non_zeros};
    S(0, 0) = 2;
    S(0, 1) = -1;

    S(1, 0) = -1;
    S(1, 1) = 2;
    S(1, 2) = -1;

    S(2, 1) = -1;
    S(2, 2) = 2;
    S(2, 3) = -1;

    S(3, 2) = -1;
    S(3, 3) = 2;
    S(3, 4) = -1;

    S(4, 3) = -1;
    S(4, 4) = 2;

    CRSMatrix A{S.getCRSMatrix()};

    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 = zero;

    CG{A, x, b, 1e-12, 1, false};
    Vector error = x - x_exact;
    REQUIRE(std::sqrt(dot(error, error)) > 1E-10 * std::sqrt(dot(x, x)));
  }
}