#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}; 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)); } } 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)); } } } } }