Skip to content
Snippets Groups Projects
Commit 548029ae authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Add minimalist EigenvalueSolver tests

parent ef90e8c1
No related branches found
No related tags found
1 merge request!93Do not initializa Kokkos Arrays anymore
...@@ -63,6 +63,7 @@ add_executable (unit_tests ...@@ -63,6 +63,7 @@ add_executable (unit_tests
test_DiscreteFunctionDescriptorP0Vector.cpp test_DiscreteFunctionDescriptorP0Vector.cpp
test_DiscreteFunctionType.cpp test_DiscreteFunctionType.cpp
test_DoWhileProcessor.cpp test_DoWhileProcessor.cpp
test_EigenvalueSolver.cpp
test_EmbeddedData.cpp test_EmbeddedData.cpp
test_EscapedString.cpp test_EscapedString.cpp
test_Exceptions.cpp test_Exceptions.cpp
......
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
#include <algebra/CRSMatrix.hpp>
#include <algebra/DenseMatrix.hpp>
#include <algebra/EigenvalueSolver.hpp>
#include <algebra/SparseMatrixDescriptor.hpp>
#include <utils/pugs_config.hpp>
// clazy:excludeall=non-pod-global-static
TEST_CASE("EigenvalueSolver", "[algebra]")
{
SECTION("Sparse Matrices")
{
SECTION("symmetric system")
{
SparseMatrixDescriptor S{3};
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};
SECTION("eigenvalues")
{
Array<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")
{
Array<double> eigenvalue_list;
std::vector<Vector<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));
DenseMatrix<double> P{3};
DenseMatrix<double> PT{3};
DenseMatrix<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];
}
DenseMatrix 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")
{
Array<double> eigenvalue_list;
DenseMatrix<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));
DenseMatrix<double> D{3};
D = zero;
for (size_t i = 0; i < 3; ++i) {
D(i, i) = eigenvalue_list[i];
}
DenseMatrix PT = transpose(P);
DenseMatrix 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
}
}
}
}
...@@ -9,6 +9,7 @@ ...@@ -9,6 +9,7 @@
#include <mesh/SynchronizerManager.hpp> #include <mesh/SynchronizerManager.hpp>
#include <utils/Messenger.hpp> #include <utils/Messenger.hpp>
#include <utils/PETScWrapper.hpp> #include <utils/PETScWrapper.hpp>
#include <utils/SLEPcWrapper.hpp>
#include <MeshDataBaseForTests.hpp> #include <MeshDataBaseForTests.hpp>
...@@ -19,6 +20,7 @@ main(int argc, char* argv[]) ...@@ -19,6 +20,7 @@ main(int argc, char* argv[])
Kokkos::initialize({4, -1, -1, true}); Kokkos::initialize({4, -1, -1, true});
PETScWrapper::initialize(argc, argv); PETScWrapper::initialize(argc, argv);
SLEPcWrapper::initialize(argc, argv);
Catch::Session session; Catch::Session session;
int result = session.applyCommandLine(argc, argv); int result = session.applyCommandLine(argc, argv);
...@@ -53,6 +55,7 @@ main(int argc, char* argv[]) ...@@ -53,6 +55,7 @@ main(int argc, char* argv[])
} }
} }
SLEPcWrapper::finalize();
PETScWrapper::finalize(); PETScWrapper::finalize();
Kokkos::finalize(); Kokkos::finalize();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment