diff --git a/src/algebra/CMakeLists.txt b/src/algebra/CMakeLists.txt index 5fd848bca482abf1285b1b36e7572deb2708e398..6310f48dcf69eb93e2c1178e05d32dbc7fe5afad 100644 --- a/src/algebra/CMakeLists.txt +++ b/src/algebra/CMakeLists.txt @@ -2,6 +2,13 @@ add_library( PugsAlgebra + EigenvalueSolver.cpp LinearSolver.cpp LinearSolverOptions.cpp PETScUtils.cpp) + +target_link_libraries( + PugsAlgebra + ${PETSC_LIBRARIES} + ${SLEPC_LIBRARIES} +) diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1d96d01f992d22a8346417dce85750d9b186749b --- /dev/null +++ b/src/algebra/EigenvalueSolver.cpp @@ -0,0 +1,152 @@ +#include <algebra/EigenvalueSolver.hpp> +#include <utils/pugs_config.hpp> + +#ifdef PUGS_HAS_SLEPC +#include <slepc.h> + +struct EigenvalueSolver::Internals +{ + static PetscReal + computeSmallestRealEigenvalueOfSymmetricMatrix(EPS& eps) + { + EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL); + EPSSolve(eps); + + PetscReal smallest_eigenvalue; + EPSGetEigenpair(eps, 0, &smallest_eigenvalue, nullptr, nullptr, nullptr); + return smallest_eigenvalue; + } + + static PetscReal + computeLargestRealEigenvalueOfSymmetricMatrix(EPS& eps) + { + EPSSetWhichEigenpairs(eps, EPS_LARGEST_REAL); + EPSSolve(eps); + + PetscReal largest_eigenvalue; + EPSGetEigenpair(eps, 0, &largest_eigenvalue, nullptr, nullptr, nullptr); + return largest_eigenvalue; + } + + static void + computeAllEigenvaluesOfSymmetricMatrixInInterval(EPS& eps, const PetscReal left_bound, const PetscReal right_bound) + { + Assert(left_bound < right_bound); + EPSSetWhichEigenpairs(eps, EPS_ALL); + EPSSetInterval(eps, left_bound - 0.01 * std::abs(left_bound), right_bound + 0.01 * std::abs(right_bound)); + + ST st; + EPSGetST(eps, &st); + STSetType(st, STSINVERT); + + KSP ksp; + STGetKSP(st, &ksp); + KSPSetType(ksp, KSPPREONLY); + + PC pc; + KSPGetPC(ksp, &pc); + PCSetType(pc, PCCHOLESKY); + EPSSetFromOptions(eps); + + EPSSolve(eps); + } + + static void + computeAllEigenvaluesOfSymmetricMatrix(EPS& eps) + { + const PetscReal left_bound = computeSmallestRealEigenvalueOfSymmetricMatrix(eps); + const PetscReal right_bound = computeLargestRealEigenvalueOfSymmetricMatrix(eps); + + computeAllEigenvaluesOfSymmetricMatrixInInterval(eps, left_bound - 0.01 * std::abs(left_bound), + right_bound + 0.01 * std::abs(right_bound)); + } +}; + +void +EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<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 = Array<double>(nb_eigenvalues); + for (PetscInt i = 0; i < nb_eigenvalues; ++i) { + EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, nullptr, nullptr); + } + + EPSDestroy(&eps); +} + +void +EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, + Array<double>& eigenvalues, + std::vector<Vector<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 = Array<double>(nb_eigenvalues); + eigenvector_list.reserve(nb_eigenvalues); + for (PetscInt i = 0; i < nb_eigenvalues; ++i) { + Vec Vr; + Vector<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); +} + +void +EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, + Array<double>& eigenvalues, + DenseMatrix<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 = Array<double>(nb_eigenvalues); + P = DenseMatrix<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 + +EigenvalueSolver::EigenvalueSolver() {} diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp new file mode 100644 index 0000000000000000000000000000000000000000..a9dd941b9bd87d3f0a79be5f36e5683ae8a1afb6 --- /dev/null +++ b/src/algebra/EigenvalueSolver.hpp @@ -0,0 +1,68 @@ +#ifndef EIGENVALUE_SOLVER_HPP +#define EIGENVALUE_SOLVER_HPP + +#include <algebra/PETScUtils.hpp> + +#include <algebra/DenseMatrix.hpp> +#include <algebra/Vector.hpp> +#include <utils/Array.hpp> +#include <utils/Exceptions.hpp> + +class EigenvalueSolver +{ + private: + struct Internals; + +#ifdef PUGS_HAS_SLEPC + void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues); + + void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, + Array<double>& eigenvalues, + std::vector<Vector<double>>& eigenvectors); + + void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues, DenseMatrix<double>& P); +#endif // PUGS_HAS_SLEPC + + public: + template <typename MatrixType> + void + computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] Array<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 + } + + template <typename MatrixType> + void + computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, + [[maybe_unused]] Array<double>& eigenvalues, + [[maybe_unused]] std::vector<Vector<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 + } + + template <typename MatrixType> + void + computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, + [[maybe_unused]] Array<double>& eigenvalues, + [[maybe_unused]] DenseMatrix<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 + } + + EigenvalueSolver(); + ~EigenvalueSolver() = default; +}; + +#endif // EIGENVALUE_SOLVER_HPP