Skip to content
Snippets Groups Projects
Commit ebbc546e authored by Emmanuel Labourasse's avatar Emmanuel Labourasse
Browse files

add exponential off a symmetric matrix using slepc

parent 7aac11e5
No related branches found
No related tags found
No related merge requests found
#include <algebra/EigenvalueSolver.hpp>
#include <cmath>
#include <utils/pugs_config.hpp>
#ifdef PUGS_HAS_SLEPC
#include <slepc.h>
#include <slepceps.h>
struct EigenvalueSolver::Internals
{
static PetscReal
......@@ -152,6 +152,43 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
EPSDestroy(&eps);
}
void
EigenvalueSolver::computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallMatrix<double>& expA)
{
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);
SmallArray<double> eigenvalues(nb_eigenvalues);
SmallMatrix<double> P(nb_eigenvalues, nb_eigenvalues);
SmallMatrix<double> D(nb_eigenvalues, nb_eigenvalues);
D = zero;
expA = 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];
}
}
for (size_t i = 0; i < eigenvalues.size(); ++i) {
D(i, i) = std::exp(eigenvalues[i]);
}
expA = P * D * transpose(P);
EPSDestroy(&eps);
}
#endif // PUGS_HAS_SLEPC
EigenvalueSolver::EigenvalueSolver() {}
......@@ -23,6 +23,7 @@ class EigenvalueSolver
void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
SmallArray<double>& eigenvalues,
SmallMatrix<double>& P);
void computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallMatrix<double>& expA);
#endif // PUGS_HAS_SLEPC
public:
......@@ -63,6 +64,17 @@ class EigenvalueSolver
#endif // PUGS_HAS_SLEPC
}
template <typename MatrixType>
void
computeExpForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] SmallMatrix<double>& expA)
{
#ifdef PUGS_HAS_SLEPC
this->computeExpForSymmetricMatrix(PETScAijMatrixEmbedder{A}, expA);
#else // PUGS_HAS_SLEPC
throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
#endif // PUGS_HAS_SLEPC
}
EigenvalueSolver();
~EigenvalueSolver() = default;
};
......
......@@ -114,6 +114,34 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
#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, expP),
"not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif // PUGS_HAS_SLEPC
}
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment