#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("Sparse Matrices")
  {
    SECTION("symmetric system")
    {
      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{}.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{}.computeForSymmetricMatrix(A, eigenvalues),
                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
      }

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

#ifdef PUGS_HAS_SLEPC
        EigenvalueSolver{}.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<double> eigen{3};
        eigen = zero;

        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{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
      }

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

#ifdef PUGS_HAS_SLEPC
        EigenvalueSolver{}.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);

        TinyMatrix<3, 3, double> TinyA;
        for (size_t i = 0; i < 3; ++i) {
          for (size_t j = 0; j < 3; ++j) {
            TinyA(i, j) = S(i, j);
          }
        }

        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{}.computeForSymmetricMatrix(A, eigenvalue_list, P),
                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
      }
      SECTION("exponential of a matrix")
      {
        SmallMatrix<double> expA{};
        SmallMatrix<double> expS{3};

        expS(0, 0) = 1325.074593930307812;
        expS(0, 1) = 662.353357244568285;
        expS(0, 2) = 1324.70671448913637;
        expS(1, 0) = expS(0, 1);
        expS(1, 1) = 331.544558063455535;
        expS(1, 2) = 662.353357244568185;
        expS(2, 0) = expS(0, 2);
        expS(2, 1) = expS(1, 2);
        expS(2, 2) = expS(0, 0);

#ifdef PUGS_HAS_SLEPC
        EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA);

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

#else    //  PUGS_HAS_SLEPC
        REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA),
                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
      }
    }
    SECTION("symmetric tiny matrix")
    {
      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> eigenvectors0{1, 0, 0, 0, 1, 0, 0, 0, 1};
      TinyVector<3> eigenvalues0{0, 0, 0};
      auto [eigenvalues, eigenmatrix] = EigenvalueSolver{}.findEigen(TestA);
    }
  }
}