diff --git a/src/algebra/CMakeLists.txt b/src/algebra/CMakeLists.txt index d215b09b50d9933f84aab02855388bb4bd25dced..78c041c5cee2e297258d299152d5cc3a4880f599 100644 --- a/src/algebra/CMakeLists.txt +++ b/src/algebra/CMakeLists.txt @@ -3,6 +3,7 @@ add_library( PugsAlgebra EigenvalueSolver.cpp + EigenvalueSolverOptions.cpp LinearSolver.cpp LinearSolverOptions.cpp PETScUtils.cpp) diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp index c71be58ddcb1c37498a172f30a1eb39b89f43aab..5888d1c0b3ecba96e908984f6d6d5f022e8be548 100644 --- a/src/algebra/EigenvalueSolver.cpp +++ b/src/algebra/EigenvalueSolver.cpp @@ -2,10 +2,13 @@ #include <utils/pugs_config.hpp> #ifdef PUGS_HAS_SLEPC +#include <algebra/PETScUtils.hpp> #include <slepc.h> +#endif // PUGS_HAS_SLEPC struct EigenvalueSolver::Internals { +#ifdef PUGS_HAS_SLEPC static PetscReal computeSmallestRealEigenvalueOfSymmetricMatrix(EPS& eps) { @@ -65,93 +68,210 @@ struct EigenvalueSolver::Internals computeAllEigenvaluesOfSymmetricMatrixInInterval(eps, left_bound - 0.01 * std::abs(left_bound), right_bound + 0.01 * std::abs(right_bound)); } -}; -void -EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues) -{ - EPS eps; + static void + slepscComputeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues) + { + EPS eps; - EPSCreate(PETSC_COMM_SELF, &eps); - EPSSetOperators(eps, A, nullptr); - EPSSetProblemType(eps, EPS_HEP); + EPSCreate(PETSC_COMM_SELF, &eps); + EPSSetOperators(eps, A, nullptr); + EPSSetProblemType(eps, EPS_HEP); - Internals::computeAllEigenvaluesOfSymmetricMatrix(eps); + Internals::computeAllEigenvaluesOfSymmetricMatrix(eps); - PetscInt nb_eigenvalues; - EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr); + PetscInt nb_eigenvalues; + EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr); - eigenvalues = SmallArray<double>(nb_eigenvalues); - for (PetscInt i = 0; i < nb_eigenvalues; ++i) { - EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, nullptr, nullptr); + eigenvalues = SmallArray<double>(nb_eigenvalues); + for (PetscInt i = 0; i < nb_eigenvalues; ++i) { + EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, nullptr, nullptr); + } + + EPSDestroy(&eps); } - EPSDestroy(&eps); + static void + slepscComputeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, + SmallArray<double>& eigenvalues, + std::vector<SmallVector<double>>& eigenvector_list) + { + EPS eps; + + EPSCreate(PETSC_COMM_SELF, &eps); + EPSSetOperators(eps, A, nullptr); + EPSSetProblemType(eps, EPS_HEP); + + Internals::computeAllEigenvaluesOfSymmetricMatrix(eps); + + PetscInt nb_eigenvalues; + EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr); + + eigenvalues = SmallArray<double>(nb_eigenvalues); + eigenvector_list.reserve(nb_eigenvalues); + for (PetscInt i = 0; i < nb_eigenvalues; ++i) { + Vec Vr; + SmallVector<double> eigenvector{A.numberOfRows()}; + VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr); + EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr); + VecDestroy(&Vr); + eigenvector_list.push_back(eigenvector); + } + + EPSDestroy(&eps); + } + + static void + slepscComputeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, + SmallArray<double>& eigenvalues, + SmallMatrix<double>& P) + { + EPS eps; + + EPSCreate(PETSC_COMM_SELF, &eps); + EPSSetOperators(eps, A, nullptr); + EPSSetProblemType(eps, EPS_HEP); + + Internals::computeAllEigenvaluesOfSymmetricMatrix(eps); + + PetscInt nb_eigenvalues; + EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr); + + eigenvalues = SmallArray<double>(nb_eigenvalues); + P = SmallMatrix<double>(nb_eigenvalues, nb_eigenvalues); + + Array<double> eigenvector(nb_eigenvalues); + for (PetscInt i = 0; i < nb_eigenvalues; ++i) { + Vec Vr; + VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr); + EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr); + VecDestroy(&Vr); + for (size_t j = 0; j < eigenvector.size(); ++j) { + P(j, i) = eigenvector[j]; + } + } + + EPSDestroy(&eps); + } +#endif // PUGS_HAS_SLEPC +}; + +#ifdef PUGS_HAS_SLEPC +void +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues) +{ + Internals::slepscComputeForSymmetricMatrix(A, eigenvalues); } void -EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, - SmallArray<double>& eigenvalues, - std::vector<SmallVector<double>>& eigenvector_list) +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues) { - EPS eps; - - EPSCreate(PETSC_COMM_SELF, &eps); - EPSSetOperators(eps, A, nullptr); - EPSSetProblemType(eps, EPS_HEP); - - Internals::computeAllEigenvaluesOfSymmetricMatrix(eps); - - PetscInt nb_eigenvalues; - EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr); - - eigenvalues = SmallArray<double>(nb_eigenvalues); - eigenvector_list.reserve(nb_eigenvalues); - for (PetscInt i = 0; i < nb_eigenvalues; ++i) { - Vec Vr; - SmallVector<double> eigenvector{A.numberOfRows()}; - VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr); - EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr); - VecDestroy(&Vr); - eigenvector_list.push_back(eigenvector); - } + Internals::slepscComputeForSymmetricMatrix(A, eigenvalues); +} - EPSDestroy(&eps); +void +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, + SmallArray<double>& eigenvalues, + std::vector<SmallVector<double>>& eigenvectors) +{ + Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors); } void -EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, - SmallArray<double>& eigenvalues, - SmallMatrix<double>& P) +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, + SmallArray<double>& eigenvalues, + std::vector<SmallVector<double>>& eigenvectors) { - EPS eps; + Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors); +} - EPSCreate(PETSC_COMM_SELF, &eps); - EPSSetOperators(eps, A, nullptr); - EPSSetProblemType(eps, EPS_HEP); +void +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, + SmallArray<double>& eigenvalues, + SmallMatrix<double>& P) +{ + Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, P); +} - Internals::computeAllEigenvaluesOfSymmetricMatrix(eps); +void +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, + SmallArray<double>& eigenvalues, + SmallMatrix<double>& P) +{ + Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, P); +} - PetscInt nb_eigenvalues; - EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr); +#else // PUGS_HAS_SLEPC - eigenvalues = SmallArray<double>(nb_eigenvalues); - P = SmallMatrix<double>(nb_eigenvalues, nb_eigenvalues); +void +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&) +{} - Array<double> eigenvector(nb_eigenvalues); - for (PetscInt i = 0; i < nb_eigenvalues; ++i) { - Vec Vr; - VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr); - EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr); - VecDestroy(&Vr); - for (size_t j = 0; j < eigenvector.size(); ++j) { - P(j, i) = eigenvector[j]; - } - } +void +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&, SmallArray<double>&) +{} - EPSDestroy(&eps); -} +void +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>&, + SmallArray<double>&, + std::vector<SmallVector<double>>&) +{} + +void +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&, + SmallArray<double>&, + std::vector<SmallVector<double>>&) +{} + +void +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&, SmallMatrix<double>&) +{} + +void +EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&, + SmallArray<double>&, + SmallMatrix<double>&) +{} #endif // PUGS_HAS_SLEPC -EigenvalueSolver::EigenvalueSolver() {} +bool +EigenvalueSolver::hasLibrary(const ESLibrary library) +{ + switch (library) { + case ESLibrary::eigen3: { +#ifdef PUGS_HAS_EIGEN3 + return true; +#else + return false; +#endif + } + case ESLibrary::slepsc: { +#ifdef PUGS_HAS_SLEPC + return true; +#else + return false; +#endif + } + // LCOV_EXCL_START + default: { + throw UnexpectedError("Eigenvalue library (" + ::name(library) + ") was not set!"); + } + // LCOV_EXCL_STOP + } +} + +void +EigenvalueSolver::checkHasLibrary(const ESLibrary library) +{ + if (not hasLibrary(library)) { + // LCOV_EXCL_START + throw NormalError(::name(library) + " is not linked to pugs. Cannot use it!"); + // LCOV_EXCL_STOP + } +} + +EigenvalueSolver::EigenvalueSolver(const EigenvalueSolverOptions& options) : m_options{options} +{ + checkHasLibrary(options.library()); +} diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp index a9b462509514ace4b736424c65f7c9abe4572f26..eebd536e4417ee953ce8cc6317a9d5a53d07bce6 100644 --- a/src/algebra/EigenvalueSolver.hpp +++ b/src/algebra/EigenvalueSolver.hpp @@ -1,10 +1,13 @@ #ifndef EIGENVALUE_SOLVER_HPP #define EIGENVALUE_SOLVER_HPP -#include <algebra/PETScUtils.hpp> +#include <algebra/EigenvalueSolverOptions.hpp> +#include <algebra/CRSMatrix.hpp> +#include <algebra/DenseMatrixWrapper.hpp> #include <algebra/SmallMatrix.hpp> #include <algebra/SmallVector.hpp> +#include <algebra/TinyMatrix.hpp> #include <utils/Exceptions.hpp> #include <utils/SmallArray.hpp> @@ -13,57 +16,97 @@ class EigenvalueSolver private: struct Internals; -#ifdef PUGS_HAS_SLEPC - void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues); + const EigenvalueSolverOptions m_options; - void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, - SmallArray<double>& eigenvalues, - std::vector<SmallVector<double>>& eigenvectors); + void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues); - void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, - SmallArray<double>& eigenvalues, - SmallMatrix<double>& P); -#endif // PUGS_HAS_SLEPC + void _slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues); + + void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, + SmallArray<double>& eigenvalues, + std::vector<SmallVector<double>>& eigenvectors); + + void _slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, + SmallArray<double>& eigenvalues, + std::vector<SmallVector<double>>& eigenvectors); + + void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, + SmallArray<double>& eigenvalues, + SmallMatrix<double>& P); + + void _slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, + SmallArray<double>& eigenvalues, + SmallMatrix<double>& P); public: + static bool hasLibrary(ESLibrary library); + static void checkHasLibrary(const ESLibrary library); + template <typename MatrixType> void - computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] SmallArray<double>& eigenvalues) + computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues) { -#ifdef PUGS_HAS_SLEPC - this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues); -#else // PUGS_HAS_SLEPC - throw NotImplementedError("SLEPc is required to solve eigenvalue problems"); -#endif // PUGS_HAS_SLEPC + Assert((eigenvalues.size() - A.numberOfRows() == 0) and A.isSquare()); + + switch (m_options.library()) { + case ESLibrary::eigen3: { + throw NotImplementedError("Eigen3"); + } + case ESLibrary::slepsc: { + this->_slepscComputeForSymmetricMatrix(A, eigenvalues); + break; + } + default: { + throw UnexpectedError(::name(m_options.library()) + " cannot compute for symmetric matrices"); + } + } } template <typename MatrixType> void - computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, - [[maybe_unused]] SmallArray<double>& eigenvalues, - [[maybe_unused]] std::vector<SmallVector<double>>& eigenvectors) + computeForSymmetricMatrix(const MatrixType& A, + SmallArray<double>& eigenvalues, + std::vector<SmallVector<double>>& eigenvectors) { -#ifdef PUGS_HAS_SLEPC - this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, eigenvectors); -#else // PUGS_HAS_SLEPC - throw NotImplementedError("SLEPc is required to solve eigenvalue problems"); -#endif // PUGS_HAS_SLEPC + Assert((eigenvalues.size() - A.numberOfRows() == 0) and (A.numberOfRows() - eigenvectors.size() == 0) and + A.isSquare()); + + switch (m_options.library()) { + case ESLibrary::eigen3: { + throw NotImplementedError("Eigen3"); + } + case ESLibrary::slepsc: { + this->_slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors); + break; + } + default: { + throw UnexpectedError(::name(m_options.library()) + " cannot compute for symmetric matrices"); + } + } } template <typename MatrixType> void - computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, - [[maybe_unused]] SmallArray<double>& eigenvalues, - [[maybe_unused]] SmallMatrix<double>& P) + computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues, SmallMatrix<double>& P) { -#ifdef PUGS_HAS_SLEPC - this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, P); -#else // PUGS_HAS_SLEPC - throw NotImplementedError("SLEPc is required to solve eigenvalue problems"); -#endif // PUGS_HAS_SLEPC + Assert((eigenvalues.size() - A.numberOfRows() == 0) and (A.numberOfRows() - P.numberOfRows() == 0) and + A.isSquare() and P.isSquare()); + + switch (m_options.library()) { + case ESLibrary::eigen3: { + throw NotImplementedError("Eigen3"); + } + case ESLibrary::slepsc: { + this->_slepscComputeForSymmetricMatrix(A, eigenvalues, P); + break; + } + default: { + throw UnexpectedError(::name(m_options.library()) + " cannot compute for symmetric matrices"); + } + } } - EigenvalueSolver(); + EigenvalueSolver(const EigenvalueSolverOptions& options = EigenvalueSolverOptions::default_options); ~EigenvalueSolver() = default; }; diff --git a/src/algebra/EigenvalueSolverOptions.cpp b/src/algebra/EigenvalueSolverOptions.cpp new file mode 100644 index 0000000000000000000000000000000000000000..2a4ef9549789da993760d72a5d7f289a76d9f2d6 --- /dev/null +++ b/src/algebra/EigenvalueSolverOptions.cpp @@ -0,0 +1,10 @@ +#include <algebra/EigenvalueSolverOptions.hpp> + +#include <rang.hpp> + +std::ostream& +operator<<(std::ostream& os, const EigenvalueSolverOptions& options) +{ + os << " library: " << rang::style::bold << name(options.library()) << rang::style::reset << '\n'; + return os; +} diff --git a/src/algebra/EigenvalueSolverOptions.hpp b/src/algebra/EigenvalueSolverOptions.hpp new file mode 100644 index 0000000000000000000000000000000000000000..d9be86be087993197949196569e341f694386a17 --- /dev/null +++ b/src/algebra/EigenvalueSolverOptions.hpp @@ -0,0 +1,108 @@ +#ifndef EIGENVALUE_SOLVER_OPTIONS_HPP +#define EIGENVALUE_SOLVER_OPTIONS_HPP + +#include <utils/Exceptions.hpp> + +#include <iostream> + +enum class ESLibrary : int8_t +{ + ES__begin = 0, + // + eigen3 = ES__begin, + slepsc, + // + ES__end +}; +; + +inline std::string +name(const ESLibrary library) +{ + switch (library) { + case ESLibrary::eigen3: { + return "Eigen3"; + } + case ESLibrary::slepsc: { + return "SLEPSc"; + } + case ESLibrary::ES__end: { + } + } + throw UnexpectedError("Eigenvalue solver library name is not defined!"); +} + +template <typename ESEnumType> +inline ESEnumType +getESEnumFromName(const std::string& enum_name) +{ + using BaseT = std::underlying_type_t<ESEnumType>; + for (BaseT enum_value = static_cast<BaseT>(ESEnumType::ES__begin); + enum_value < static_cast<BaseT>(ESEnumType::ES__end); ++enum_value) { + if (name(ESEnumType{enum_value}) == enum_name) { + return ESEnumType{enum_value}; + } + } + throw NormalError(std::string{"could not find '"} + enum_name + "' associate type!"); +} + +template <typename ESEnumType> +inline void +printESEnumListNames(std::ostream& os) +{ + using BaseT = std::underlying_type_t<ESEnumType>; + for (BaseT enum_value = static_cast<BaseT>(ESEnumType::ES__begin); + enum_value < static_cast<BaseT>(ESEnumType::ES__end); ++enum_value) { + os << " - " << name(ESEnumType{enum_value}) << '\n'; + } +} + +class EigenvalueSolverOptions +{ + private: + ESLibrary m_library = ESLibrary::slepsc; + + double m_epsilon = 1E-6; + size_t m_maximum_iteration = 200; + + bool m_verbose = false; + + public: + static EigenvalueSolverOptions default_options; + + friend std::ostream& operator<<(std::ostream& os, const EigenvalueSolverOptions& options); + + ESLibrary& + library() + { + return m_library; + } + + ESLibrary + library() const + { + return m_library; + } + + bool& + verbose() + { + return m_verbose; + }; + + bool + verbose() const + { + return m_verbose; + }; + + EigenvalueSolverOptions(const EigenvalueSolverOptions&) = default; + EigenvalueSolverOptions(EigenvalueSolverOptions&&) = default; + + EigenvalueSolverOptions() = default; + ~EigenvalueSolverOptions() = default; +}; + +inline EigenvalueSolverOptions EigenvalueSolverOptions::default_options; + +#endif // EIGENVALUE_SOLVER_OPTIONS_HPP diff --git a/src/language/modules/LinearSolverModule.cpp b/src/language/modules/LinearSolverModule.cpp index d3e2de0faeb52c213f372e3a751f075ea2d735d7..1ff8be36d3d0649caaae56e250477fda2ee7690d 100644 --- a/src/language/modules/LinearSolverModule.cpp +++ b/src/language/modules/LinearSolverModule.cpp @@ -2,7 +2,8 @@ #include <language/modules/ModuleRepository.hpp> -#include <algebra/LinearSolver.hpp> +#include <algebra/EigenvalueSolverOptions.hpp> +#include <algebra/LinearSolverOptions.hpp> #include <language/utils/BuiltinFunctionEmbedder.hpp> #include <language/utils/TypeDescriptor.hpp> @@ -90,6 +91,43 @@ LinearSolverModule::LinearSolverModule() return os.str(); } + )); + + this->_addBuiltinFunction("setESLibrary", std::function( + + [](const std::string& library_name) -> void { + EigenvalueSolverOptions::default_options.library() = + getESEnumFromName<ESLibrary>(library_name); + } + + )); + + this->_addBuiltinFunction("getESOptions", std::function( + + []() -> std::string { + std::ostringstream os; + os << rang::fgB::yellow << "Eigenvalue solver options" + << rang::style::reset << '\n'; + os << EigenvalueSolverOptions::default_options; + return os.str(); + } + + )); + + this->_addBuiltinFunction("getESAvailable", std::function( + + []() -> std::string { + std::ostringstream os; + + os << rang::fgB::yellow << "Available eigenvalue solver options" + << rang::style::reset << '\n'; + os << rang::fgB::blue << " libraries" << rang::style::reset << '\n'; + + printESEnumListNames<ESLibrary>(os); + + return os.str(); + } + )); } diff --git a/src/utils/checkpointing/Checkpoint.cpp b/src/utils/checkpointing/Checkpoint.cpp index 5976cb2fb5b3c5024bc56a1b77d384abd698e97c..ed56dcb83b4cc33fc0cdbaac4c6e704423f13336 100644 --- a/src/utils/checkpointing/Checkpoint.cpp +++ b/src/utils/checkpointing/Checkpoint.cpp @@ -4,6 +4,7 @@ #ifdef PUGS_HAS_HDF5 +#include <algebra/EigenvalueSolverOptions.hpp> #include <algebra/LinearSolverOptions.hpp> #include <dev/ParallelChecker.hpp> #include <language/ast/ASTExecutionStack.hpp> @@ -21,12 +22,12 @@ #include <utils/HighFivePugsUtils.hpp> #include <utils/RandomEngine.hpp> #include <utils/checkpointing/DualMeshTypeHFType.hpp> +#include <utils/checkpointing/EigenvalueSolverOptionsHFType.hpp> #include <utils/checkpointing/LinearSolverOptionsHFType.hpp> #include <utils/checkpointing/ParallelCheckerHFType.hpp> #include <utils/checkpointing/ResumingManager.hpp> #include <iostream> -#include <map> void checkpoint() @@ -122,6 +123,14 @@ checkpoint() linear_solver_options_default_group.createAttribute("method", default_options.method()); linear_solver_options_default_group.createAttribute("precond", default_options.precond()); } + { + HighFive::Group eigenvalue_solver_options_default_group = + checkpoint.createGroup("singleton/eigenvalue_solver_options_default"); + + const EigenvalueSolverOptions& default_options = EigenvalueSolverOptions::default_options; + + eigenvalue_solver_options_default_group.createAttribute("library", default_options.library()); + } { const auto& primal_to_dual_connectivity_info_map = DualConnectivityManager::instance().primalToDualConnectivityInfoMap(); diff --git a/src/utils/checkpointing/EigenvalueSolverOptionsHFType.hpp b/src/utils/checkpointing/EigenvalueSolverOptionsHFType.hpp new file mode 100644 index 0000000000000000000000000000000000000000..32995c1f5475b04534ac42f344d9e5e9ccef1fa2 --- /dev/null +++ b/src/utils/checkpointing/EigenvalueSolverOptionsHFType.hpp @@ -0,0 +1,15 @@ +#ifndef EIGENVALUE_SOLVER_OPTIONS_HF_TYPE_HPP +#define EIGENVALUE_SOLVER_OPTIONS_HF_TYPE_HPP + +#include <algebra/EigenvalueSolverOptions.hpp> +#include <utils/HighFivePugsUtils.hpp> +#include <utils/PugsMacros.hpp> + +HighFive::EnumType<ESLibrary> PUGS_INLINE +create_enum_ESOptions_library_type() +{ + return {{"eigen3", ESLibrary::eigen3}, {"slepsc", ESLibrary::slepsc}}; +} +HIGHFIVE_REGISTER_TYPE(ESLibrary, create_enum_ESOptions_library_type) + +#endif // EIGENVALUE_SOLVER_OPTIONS_HF_TYPE_HPP diff --git a/src/utils/checkpointing/LinearSolverOptionsHFType.hpp b/src/utils/checkpointing/LinearSolverOptionsHFType.hpp index 94f88fc60088c297479af36975d533d1c94bdd71..1c94d74a948a27f9bc8fc54915ef819171379152 100644 --- a/src/utils/checkpointing/LinearSolverOptionsHFType.hpp +++ b/src/utils/checkpointing/LinearSolverOptionsHFType.hpp @@ -8,7 +8,7 @@ HighFive::EnumType<LSLibrary> PUGS_INLINE create_enum_LSOptions_library_type() { - return {{"builtin", LSLibrary::builtin}, {"petsc", LSLibrary::petsc}}; + return {{"builtin", LSLibrary::builtin}, {"eigen3", LSLibrary::eigen3}, {"petsc", LSLibrary::petsc}}; } HIGHFIVE_REGISTER_TYPE(LSLibrary, create_enum_LSOptions_library_type) diff --git a/src/utils/checkpointing/Resume.cpp b/src/utils/checkpointing/Resume.cpp index 9d58e1d9dec4050f218d3efad71e5194ccaa5111..effd6b8080c5f5c8c3d64bbfe7d30fafbbbb1e7e 100644 --- a/src/utils/checkpointing/Resume.cpp +++ b/src/utils/checkpointing/Resume.cpp @@ -4,6 +4,7 @@ #ifdef PUGS_HAS_HDF5 +#include <algebra/EigenvalueSolverOptions.hpp> #include <algebra/LinearSolverOptions.hpp> #include <language/ast/ASTExecutionStack.hpp> #include <language/utils/ASTCheckpointsInfo.hpp> @@ -14,6 +15,7 @@ #include <utils/ExecutionStatManager.hpp> #include <utils/HighFivePugsUtils.hpp> #include <utils/RandomEngine.hpp> +#include <utils/checkpointing/EigenvalueSolverOptionsHFType.hpp> #include <utils/checkpointing/LinearSolverOptionsHFType.hpp> #include <utils/checkpointing/ParallelCheckerHFType.hpp> #include <utils/checkpointing/ResumingData.hpp> @@ -81,6 +83,14 @@ resume() default_options.method() = linear_solver_options_default_group.getAttribute("method").read<LSMethod>(); default_options.precond() = linear_solver_options_default_group.getAttribute("precond").read<LSPrecond>(); } + { + HighFive::Group eigenvalue_solver_options_default_group = + checkpoint.getGroup("singleton/eigenvalue_solver_options_default"); + + EigenvalueSolverOptions& default_options = EigenvalueSolverOptions::default_options; + + default_options.library() = eigenvalue_solver_options_default_group.getAttribute("library").read<ESLibrary>(); + } checkpointing::ResumingData::instance().readData(checkpoint, p_symbol_table); diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp index fb7bf07907bf6c9e1db9333c9713365350918b06..4bb1d339d3004694ef646c4914b343b55c591384 100644 --- a/tests/test_EigenvalueSolver.cpp +++ b/tests/test_EigenvalueSolver.cpp @@ -12,109 +12,321 @@ TEST_CASE("EigenvalueSolver", "[algebra]") { - SECTION("Sparse Matrices") + SECTION("SLEPc") { - SECTION("symmetric system") + 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}; + 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(0, 0) = 3; + S(0, 1) = 2; + S(0, 2) = 4; - S(1, 0) = 2; - S(1, 1) = 0; - S(1, 2) = 2; + S(1, 0) = 2; + S(1, 1) = 0; + S(1, 2) = 2; - S(2, 0) = 4; - S(2, 1) = 2; - S(2, 2) = 3; + S(2, 0) = 4; + S(2, 1) = 2; + S(2, 2) = 3; - CRSMatrix A{S.getCRSMatrix()}; + CRSMatrix A{S.getCRSMatrix()}; - SECTION("eigenvalues") - { - SmallArray<double> eigenvalues; + 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{}.computeForSymmetricMatrix(A, eigenvalues); - REQUIRE(eigenvalues[0] == Catch::Approx(-1)); - REQUIRE(eigenvalues[1] == Catch::Approx(-1)); - REQUIRE(eigenvalues[2] == Catch::Approx(8)); + 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{}.computeForSymmetricMatrix(A, eigenvalues), - "not implemented yet: SLEPc is required to solve eigenvalue problems"); + 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("eigenvalues and eigenvectors") + SECTION("SmallMatrix") + { + SECTION("symmetric system") { - SmallArray<double> eigenvalue_list; - std::vector<SmallVector<double>> eigenvector_list; + 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{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list); + 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)); + 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; + 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]; + 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]; } - 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)); + 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{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list), - "not implemented yet: SLEPc is required to solve eigenvalue problems"); + 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("eigenvalues and passage matrix") + SECTION("TinyMatrix") + { + SECTION("symmetric system") { - SmallArray<double> eigenvalue_list; - SmallMatrix<double> P{}; + TinyMatrix<3> A = zero; -#ifdef PUGS_HAS_SLEPC - EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P); + A(0, 0) = 3; + A(0, 1) = 2; + A(0, 2) = 4; - REQUIRE(eigenvalue_list[0] == Catch::Approx(-1)); - REQUIRE(eigenvalue_list[1] == Catch::Approx(-1)); - REQUIRE(eigenvalue_list[2] == Catch::Approx(8)); + A(1, 0) = 2; + A(1, 1) = 0; + A(1, 2) = 2; - SmallMatrix<double> D{3}; - D = zero; - for (size_t i = 0; i < 3; ++i) { - D(i, i) = eigenvalue_list[i]; + 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 } - 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)); + 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{}.computeForSymmetricMatrix(A, eigenvalue_list, P), - "not implemented yet: SLEPc is required to solve eigenvalue problems"); + REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P), + "error: SLEPSc is not linked to pugs. Cannot use it!"); #endif // PUGS_HAS_SLEPC + } } } } diff --git a/tests/test_LinearSolverOptions.cpp b/tests/test_LinearSolverOptions.cpp index 6209597147cb03fed823935d8a2a69110018d76a..7d078a7fb73e094ca09f888ab43ed12980081f2f 100644 --- a/tests/test_LinearSolverOptions.cpp +++ b/tests/test_LinearSolverOptions.cpp @@ -68,6 +68,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]") SECTION("library name") { REQUIRE(name(LSLibrary::builtin) == "builtin"); + REQUIRE(name(LSLibrary::eigen3) == "Eigen3"); REQUIRE(name(LSLibrary::petsc) == "PETSc"); REQUIRE_THROWS_WITH(name(LSLibrary::LS__end), "unexpected error: Linear system library name is not defined!"); } diff --git a/tests/test_checkpointing_HFTypes.cpp b/tests/test_checkpointing_HFTypes.cpp index 41fc6948c382e3d8f4a3f0bf100cd94405c50875..1a5f02b4a4328910fe9caca1ff3c780954f77197 100644 --- a/tests/test_checkpointing_HFTypes.cpp +++ b/tests/test_checkpointing_HFTypes.cpp @@ -3,6 +3,7 @@ #include <utils/checkpointing/DiscreteFunctionTypeHFType.hpp> #include <utils/checkpointing/DualMeshTypeHFType.hpp> +#include <utils/checkpointing/EigenvalueSolverOptionsHFType.hpp> #include <utils/checkpointing/IBoundaryConditionDescriptorHFType.hpp> #include <utils/checkpointing/IBoundaryDescriptorHFType.hpp> #include <utils/checkpointing/IInterfaceDescriptorHFType.hpp> @@ -160,6 +161,7 @@ TEST_CASE("HFTypes", "[utils/checkpointing]") SECTION("LinearSolverOptionsHFType") { file.createAttribute("builtin", LSLibrary::builtin); + file.createAttribute("eigen3", LSLibrary::eigen3); file.createAttribute("petsc", LSLibrary::petsc); file.createAttribute("cg", LSMethod::cg); @@ -177,6 +179,7 @@ TEST_CASE("HFTypes", "[utils/checkpointing]") REQUIRE(file.getAttribute("builtin").read<LSLibrary>() == LSLibrary::builtin); REQUIRE(file.getAttribute("petsc").read<LSLibrary>() == LSLibrary::petsc); + REQUIRE(file.getAttribute("eigen3").read<LSLibrary>() == LSLibrary::eigen3); REQUIRE(file.getAttribute("cg").read<LSMethod>() == LSMethod::cg); REQUIRE(file.getAttribute("bicgstab").read<LSMethod>() == LSMethod::bicgstab);