Select Git revision
test_EigenvalueSolver.cpp
-
Emmanuel Labourasse authoredEmmanuel Labourasse authored
test_EigenvalueSolver.cpp 8.08 KiB
#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));
}
}
}
}
}