test_PCG.cpp

#include <catch2/catch.hpp>

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

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

TEST_CASE("PCG", "[algebra]")
{
  SECTION("no preconditionner")
  {
    SparseMatrixDescriptor S{5};
    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};

    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;

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

  SECTION("singular matrix with RHS in its image")
  {
    SparseMatrixDescriptor<double, uint8_t> S{5};

    CRSMatrix<double, uint8_t> A{S};

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

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

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

  SECTION("no preconditionner non-converged")
  {
    SparseMatrixDescriptor S{5};
    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};

    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;

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