diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8fc35b90dd026596316a2a478747cf896f15e7ec..a56fb1f9d0b49108f1770f186cec453758cfd793 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -63,6 +63,7 @@ add_executable (unit_tests test_DiscreteFunctionDescriptorP0Vector.cpp test_DiscreteFunctionType.cpp test_DoWhileProcessor.cpp + test_EigenvalueSolver.cpp test_EmbeddedData.cpp test_EscapedString.cpp test_Exceptions.cpp diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..e779bce5600e32cf26b0fde7339693f8da69bcb1 --- /dev/null +++ b/tests/test_EigenvalueSolver.cpp @@ -0,0 +1,118 @@ +#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 + } + } + } +} diff --git a/tests/test_main.cpp b/tests/test_main.cpp index 89e55c3d7bd024c6ef5f3344283fe0288c292cc5..b5316614d10c21cbe13bb5540cc0a13d7e4b4b24 100644 --- a/tests/test_main.cpp +++ b/tests/test_main.cpp @@ -9,6 +9,7 @@ #include <mesh/SynchronizerManager.hpp> #include <utils/Messenger.hpp> #include <utils/PETScWrapper.hpp> +#include <utils/SLEPcWrapper.hpp> #include <MeshDataBaseForTests.hpp> @@ -19,6 +20,7 @@ main(int argc, char* argv[]) Kokkos::initialize({4, -1, -1, true}); PETScWrapper::initialize(argc, argv); + SLEPcWrapper::initialize(argc, argv); Catch::Session session; int result = session.applyCommandLine(argc, argv); @@ -53,6 +55,7 @@ main(int argc, char* argv[]) } } + SLEPcWrapper::finalize(); PETScWrapper::finalize(); Kokkos::finalize();