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