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

#include <algebra/CRSMatrix.hpp>
#include <algebra/CRSMatrixDescriptor.hpp>
#include <algebra/EigenvalueSolver.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/TinyMatrix.hpp>
#include <scheme/HyperelasticSolver.hpp>

#include <utils/pugs_config.hpp>

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

TEST_CASE("EigenvalueSolver", "[algebra]")
{
  SECTION("symmetric system")
  {
    SECTION("SLEPc")
    {
      EigenvalueSolverOptions options;
      options.library() = ESLibrary::slepsc;

      SECTION("Sparse Matrices")
      {
        Array<int> non_zeros(3);
        non_zeros.fill(3);
        CRSMatrixDescriptor<double> S{3, 3, non_zeros};

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

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

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

        CRSMatrix A{S.getCRSMatrix()};

        SECTION("eigenvalues")
        {
          SmallArray<double> eigenvalues;

#ifdef PUGS_HAS_SLEPC
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
          REQUIRE(eigenvalues[0] == Catch::Approx(-1));
          REQUIRE(eigenvalues[1] == Catch::Approx(-1));
          REQUIRE(eigenvalues[2] == Catch::Approx(8));
#else    //  PUGS_HAS_SLEPC
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
                              "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_SLEPC
        }

        SECTION("eigenvalues and eigenvectors")
        {
          SmallArray<double> eigenvalue_list;
          std::vector<SmallVector<double>> eigenvector_list;

#ifdef PUGS_HAS_SLEPC
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> P{3};
          SmallMatrix<double> PT{3};
          SmallMatrix<double> D{3};
          D = zero;

          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              P(i, j)  = eigenvector_list[j][i];
              PT(i, j) = eigenvector_list[i][j];
            }
            D(i, i) = eigenvalue_list[i];
          }

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_SLEPC
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
                              "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_SLEPC
        }

        SECTION("eigenvalues and transition matrix")
        {
          SmallArray<double> eigenvalue_list;
          SmallMatrix<double> P{};

#ifdef PUGS_HAS_SLEPC
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> D{3};
          D = zero;
          for (size_t i = 0; i < 3; ++i) {
            D(i, i) = eigenvalue_list[i];
          }
          SmallMatrix PT = transpose(P);

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_SLEPC
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
                              "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_SLEPC
        }
      }

      SECTION("SmallMatrix")
      {
        SmallMatrix<double> A{3, 3};
        A.fill(0);
        A(0, 0) = 3;
        A(0, 1) = 2;
        A(0, 2) = 4;

        A(1, 0) = 2;
        A(1, 1) = 0;
        A(1, 2) = 2;

        A(2, 0) = 4;
        A(2, 1) = 2;
        A(2, 2) = 3;

        SECTION("eigenvalues")
        {
          SmallArray<double> eigenvalues;

#ifdef PUGS_HAS_SLEPC
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
          REQUIRE(eigenvalues[0] == Catch::Approx(-1));
          REQUIRE(eigenvalues[1] == Catch::Approx(-1));
          REQUIRE(eigenvalues[2] == Catch::Approx(8));
#else    //  PUGS_HAS_SLEPC
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
                              "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_SLEPC
        }

        SECTION("eigenvalues and eigenvectors")
        {
          SmallArray<double> eigenvalue_list;
          std::vector<SmallVector<double>> eigenvector_list;

#ifdef PUGS_HAS_SLEPC
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> P{3};
          SmallMatrix<double> PT{3};
          SmallMatrix<double> D{3};
          D = zero;

          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              P(i, j)  = eigenvector_list[j][i];
              PT(i, j) = eigenvector_list[i][j];
            }
            D(i, i) = eigenvalue_list[i];
          }

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_SLEPC
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
                              "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_SLEPC
        }

        SECTION("eigenvalues and transition matrix")
        {
          SmallArray<double> eigenvalue_list;
          SmallMatrix<double> P{};

#ifdef PUGS_HAS_SLEPC
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> D{3};
          D = zero;
          for (size_t i = 0; i < 3; ++i) {
            D(i, i) = eigenvalue_list[i];
          }
          SmallMatrix PT = transpose(P);

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_SLEPC
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
                              "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_SLEPC
        }
      }

      SECTION("TinyMatrix")
      {
        TinyMatrix<3> A = zero;

        A(0, 0) = 3;
        A(0, 1) = 2;
        A(0, 2) = 4;

        A(1, 0) = 2;
        A(1, 1) = 0;
        A(1, 2) = 2;

        A(2, 0) = 4;
        A(2, 1) = 2;
        A(2, 2) = 3;

        SECTION("eigenvalues")
        {
          SmallArray<double> eigenvalues;

#ifdef PUGS_HAS_SLEPC
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
          REQUIRE(eigenvalues[0] == Catch::Approx(-1));
          REQUIRE(eigenvalues[1] == Catch::Approx(-1));
          REQUIRE(eigenvalues[2] == Catch::Approx(8));
#else    //  PUGS_HAS_SLEPC
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
                              "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_SLEPC
        }

        SECTION("eigenvalues and eigenvectors")
        {
          SmallArray<double> eigenvalue_list;
          std::vector<SmallVector<double>> eigenvector_list;

#ifdef PUGS_HAS_SLEPC
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> P{3};
          SmallMatrix<double> PT{3};
          SmallMatrix<double> D{3};
          D = zero;

          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              P(i, j)  = eigenvector_list[j][i];
              PT(i, j) = eigenvector_list[i][j];
            }
            D(i, i) = eigenvalue_list[i];
          }

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_SLEPC
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
                              "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_SLEPC
        }

        SECTION("eigenvalues and transition matrix")
        {
          SmallArray<double> eigenvalue_list;
          SmallMatrix<double> P{};

#ifdef PUGS_HAS_SLEPC
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> D{3};
          D = zero;
          for (size_t i = 0; i < 3; ++i) {
            D(i, i) = eigenvalue_list[i];
          }
          SmallMatrix PT = transpose(P);

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_SLEPC
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
                              "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_SLEPC
        }
      }
    }

    SECTION("Eigen3")
    {
      EigenvalueSolverOptions options;
      options.library() = ESLibrary::eigen3;

      SECTION("Sparse Matrices")
      {
        Array<int> non_zeros(3);
        non_zeros.fill(3);
        CRSMatrixDescriptor<double> S{3, 3, non_zeros};

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

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

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

        CRSMatrix A{S.getCRSMatrix()};

        SECTION("eigenvalues")
        {
          SmallArray<double> eigenvalues;

#ifdef PUGS_HAS_EIGEN3
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
          REQUIRE(eigenvalues[0] == Catch::Approx(-1));
          REQUIRE(eigenvalues[1] == Catch::Approx(-1));
          REQUIRE(eigenvalues[2] == Catch::Approx(8));
#else    //  PUGS_HAS_EIGEN3
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
                              "error: Eigen3 is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_EIGEN3
        }

        SECTION("eigenvalues and eigenvectors")
        {
          SmallArray<double> eigenvalue_list;
          std::vector<SmallVector<double>> eigenvector_list;

#ifdef PUGS_HAS_EIGEN3
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> P{3};
          SmallMatrix<double> PT{3};
          SmallMatrix<double> D{3};
          D = zero;

          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              P(i, j)  = eigenvector_list[j][i];
              PT(i, j) = eigenvector_list[i][j];
            }
            D(i, i) = eigenvalue_list[i];
          }

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_EIGEN3
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
                              "error: Eigen3 is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_EIGEN3
        }

        SECTION("eigenvalues and transition matrix")
        {
          SmallArray<double> eigenvalue_list;
          SmallMatrix<double> P{};

#ifdef PUGS_HAS_EIGEN3
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> D{3};
          D = zero;
          for (size_t i = 0; i < 3; ++i) {
            D(i, i) = eigenvalue_list[i];
          }
          SmallMatrix PT = transpose(P);

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_EIGEN3
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
                              "error: Eigen3 is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_EIGEN3
        }
      }

      SECTION("SmallMatrix")
      {
        SmallMatrix<double> A{3, 3};
        A.fill(0);
        A(0, 0) = 3;
        A(0, 1) = 2;
        A(0, 2) = 4;

        A(1, 0) = 2;
        A(1, 1) = 0;
        A(1, 2) = 2;

        A(2, 0) = 4;
        A(2, 1) = 2;
        A(2, 2) = 3;

        SECTION("eigenvalues")
        {
          SmallArray<double> eigenvalues;

#ifdef PUGS_HAS_EIGEN3
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
          REQUIRE(eigenvalues[0] == Catch::Approx(-1));
          REQUIRE(eigenvalues[1] == Catch::Approx(-1));
          REQUIRE(eigenvalues[2] == Catch::Approx(8));
#else    //  PUGS_HAS_EIGEN3
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
                              "error: Eigen3 is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_EIGEN3
        }

        SECTION("eigenvalues and eigenvectors")
        {
          SmallArray<double> eigenvalue_list;
          std::vector<SmallVector<double>> eigenvector_list;

#ifdef PUGS_HAS_EIGEN3
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> P{3};
          SmallMatrix<double> PT{3};
          SmallMatrix<double> D{3};
          D = zero;

          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              P(i, j)  = eigenvector_list[j][i];
              PT(i, j) = eigenvector_list[i][j];
            }
            D(i, i) = eigenvalue_list[i];
          }

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_EIGEN3
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
                              "error: Eigen3 is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_EIGEN3
        }

        SECTION("eigenvalues and transition matrix")
        {
          SmallArray<double> eigenvalue_list;
          SmallMatrix<double> P{};

#ifdef PUGS_HAS_EIGEN3
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> D{3};
          D = zero;
          for (size_t i = 0; i < 3; ++i) {
            D(i, i) = eigenvalue_list[i];
          }
          SmallMatrix PT = transpose(P);

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_EIGEN3
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
                              "error: Eigen3 is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_EIGEN3
        }
      }

      SECTION("TinyMatrix")
      {
        TinyMatrix<3> A = zero;

        A(0, 0) = 3;
        A(0, 1) = 2;
        A(0, 2) = 4;

        A(1, 0) = 2;
        A(1, 1) = 0;
        A(1, 2) = 2;

        A(2, 0) = 4;
        A(2, 1) = 2;
        A(2, 2) = 3;

        SECTION("eigenvalues")
        {
          SmallArray<double> eigenvalues;

#ifdef PUGS_HAS_EIGEN3
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
          REQUIRE(eigenvalues[0] == Catch::Approx(-1));
          REQUIRE(eigenvalues[1] == Catch::Approx(-1));
          REQUIRE(eigenvalues[2] == Catch::Approx(8));
#else    //  PUGS_HAS_EIGEN3
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
                              "error: Eigen3 is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_EIGEN3
        }

        SECTION("eigenvalues and eigenvectors")
        {
          SmallArray<double> eigenvalue_list;
          std::vector<SmallVector<double>> eigenvector_list;

#ifdef PUGS_HAS_EIGEN3
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> P{3};
          SmallMatrix<double> PT{3};
          SmallMatrix<double> D{3};
          D = zero;

          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              P(i, j)  = eigenvector_list[j][i];
              PT(i, j) = eigenvector_list[i][j];
            }
            D(i, i) = eigenvalue_list[i];
          }

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_EIGEN3
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
                              "error: Eigen3 is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_EIGEN3
        }

        SECTION("eigenvalues and transition matrix")
        {
          SmallArray<double> eigenvalue_list;
          SmallMatrix<double> P{};

#ifdef PUGS_HAS_EIGEN3
          EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);

          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

          SmallMatrix<double> D{3};
          D = zero;
          for (size_t i = 0; i < 3; ++i) {
            D(i, i) = eigenvalue_list[i];
          }
          SmallMatrix PT = transpose(P);

          SmallMatrix PDPT = P * D * PT;
          for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
              REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
            }
          }
#else    //  PUGS_HAS_EIGEN3
          REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
                              "error: Eigen3 is not linked to pugs. Cannot use it!");
#endif   // PUGS_HAS_EIGEN3
        }
      }
    }

    SECTION("symmetric tiny matrix")
    {
#ifdef PUGS_HAS_SLEPC
      TinyMatrix<3> TestA{3e10, 2e10, 4e10, 2e10, 0, 2e10, 4e10, 2e10, 3e10};
      TinyMatrix<3> TestA2{4, 2, 3, 2, 0, 2, 3, 2, 4};
      TinyMatrix<3> TestB{1, -1, 0, -1, 1, 0, 0, 0, 3};
      TinyMatrix<3> TestC{0, 0, 0, 0, 0, 0, 0, 0, 0};
      TinyMatrix<3> TestD{1e10, 0, 0, 0, 1, 0, 0, 0, 1};
      TinyMatrix<3> TestE{3e-10, 2e-10, 4e-10, 2e-10, 0, 2e-10, 4e-10, 2e-10, 3e-10};

      TinyMatrix<3> expA2;

      auto [eigenvalues, eigenmatrix] = EigenvalueSolver{}.findEigen(TestA);
      TinyMatrix<3> Diag              = zero;
      for (size_t i = 0; i < 3; ++i) {
        Diag(i, i) = eigenvalues[i];
      }
      TinyMatrix<3> B = eigenmatrix * Diag * transpose(eigenmatrix);
      for (size_t i = 0; i < 3; ++i) {
        for (size_t j = 0; j < 3; ++j) {
          REQUIRE((B(i, j) - TestA(i, j)) / frobeniusNorm(TestA) == Catch::Approx(0).margin(1E-8));
        }
      }
      auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestA2);
      Diag                              = zero;
      for (size_t i = 0; i < 3; ++i) {
        Diag(i, i) = eigenvalues2[i];
      }
      B = eigenmatrix2 * Diag * transpose(eigenmatrix2);
      for (size_t i = 0; i < 3; ++i) {
        for (size_t j = 0; j < 3; ++j) {
          REQUIRE((B(i, j) - TestA2(i, j)) / frobeniusNorm(TestA2) == Catch::Approx(0).margin(1E-8));
        }
      }
      expA2 = EigenvalueSolver{}.computeExpForTinyMatrix(TestA2);
      for (size_t i = 0; i < 3; ++i) {
        Diag(i, i) = std::exp(eigenvalues2[i]);
      }
      B = eigenmatrix2 * Diag * transpose(eigenmatrix2);

      for (size_t i = 0; i < 3; ++i) {
        for (size_t j = 0; j < 3; ++j) {
          REQUIRE(expA2(i, j) - B(i, j) == Catch::Approx(0).margin(1E-12));
        }
      }

      auto [eigenvalues3, eigenmatrix3] = EigenvalueSolver{}.findEigen(TestB);
      Diag                              = zero;
      for (size_t i = 0; i < 3; ++i) {
        Diag(i, i) = eigenvalues3[i];
      }
      B = eigenmatrix3 * Diag * transpose(eigenmatrix3);
      for (size_t i = 0; i < 3; ++i) {
        for (size_t j = 0; j < 3; ++j) {
          REQUIRE((B(i, j) - TestB(i, j)) / frobeniusNorm(TestB) == Catch::Approx(0).margin(1E-8));
        }
      }

      auto [eigenvalues4, eigenmatrix4] = EigenvalueSolver{}.findEigen(TestC);
      Diag                              = zero;
      for (size_t i = 0; i < 3; ++i) {
        Diag(i, i) = eigenvalues4[i];
      }
      B = eigenmatrix4 * Diag * transpose(eigenmatrix4);
      for (size_t i = 0; i < 3; ++i) {
        for (size_t j = 0; j < 3; ++j) {
          REQUIRE((B(i, j) - TestC(i, j)) == Catch::Approx(0).margin(1E-8));
        }
      }

      auto [eigenvalues5, eigenmatrix5] = EigenvalueSolver{}.findEigen(TestD);
      Diag                              = zero;
      for (size_t i = 0; i < 3; ++i) {
        Diag(i, i) = eigenvalues5[i];
      }
      B = eigenmatrix5 * Diag * transpose(eigenmatrix5);
      for (size_t i = 0; i < 3; ++i) {
        for (size_t j = 0; j < 3; ++j) {
          REQUIRE((B(i, j) - TestD(i, j)) / frobeniusNorm(TestD) == Catch::Approx(0).margin(1E-8));
        }
      }

      auto [eigenvalues6, eigenmatrix6] = EigenvalueSolver{}.findEigen(TestE);
      Diag                              = zero;
      for (size_t i = 0; i < 3; ++i) {
        Diag(i, i) = eigenvalues6[i];
      }
      B = eigenmatrix6 * Diag * transpose(eigenmatrix6);
      for (size_t i = 0; i < 3; ++i) {
        for (size_t j = 0; j < 3; ++j) {
          REQUIRE((B(i, j) - TestE(i, j)) / frobeniusNorm(TestE) == Catch::Approx(0).margin(1E-8));
        }
      }
    }
  }
#endif
}