#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 <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 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); 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, expP), "not implemented yet: SLEPc is required to solve eigenvalue problems"); #endif // PUGS_HAS_SLEPC } } } }