diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 5888d1c0b3ecba96e908984f6d6d5f022e8be548..b17e650af775530d3e61f4561b7676dccf1efdd4 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -6,6 +6,10 @@
#include <slepc.h>
#endif // PUGS_HAS_SLEPC
+#ifdef PUGS_HAS_EIGEN3
+#include <algebra/Eigen3Utils.hpp>
+#endif // PUGS_HAS_EIGEN3
+
struct EigenvalueSolver::Internals
{
#ifdef PUGS_HAS_SLEPC
@@ -154,6 +158,73 @@ struct EigenvalueSolver::Internals
EPSDestroy(&eps);
}
#endif // PUGS_HAS_SLEPC
+
+#ifdef PUGS_HAS_EIGEN3
+
+ template <typename Eigen3MatrixEmbedderType>
+ static void
+ eigen3ComputeForSymmetricMatrix(const Eigen3MatrixEmbedderType& A, SmallArray<double>& eigenvalues)
+ {
+ using Eigen3MatrixType = typename Eigen3MatrixEmbedderType::Eigen3MatrixType;
+ Eigen::SelfAdjointEigenSolver<Eigen3MatrixType> eigen_solver;
+
+ eigen_solver.compute(A.matrix(), Eigen::EigenvaluesOnly);
+
+ eigenvalues = SmallArray<double>(A.numberOfRows());
+ for (size_t i = 0; i < eigenvalues.size(); ++i) {
+ eigenvalues[i] = eigen_solver.eigenvalues()[i];
+ }
+ }
+
+ template <typename Eigen3MatrixEmbedderType>
+ static void
+ eigen3ComputeForSymmetricMatrix(const Eigen3MatrixEmbedderType& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors)
+ {
+ using Eigen3MatrixType = typename Eigen3MatrixEmbedderType::Eigen3MatrixType;
+ Eigen::SelfAdjointEigenSolver<Eigen3MatrixType> eigen_solver;
+
+ eigen_solver.compute(A.matrix(), Eigen::ComputeEigenvectors);
+
+ eigenvalues = SmallArray<double>(A.numberOfRows());
+ for (size_t i = 0; i < eigenvalues.size(); ++i) {
+ eigenvalues[i] = eigen_solver.eigenvalues()[i];
+ }
+
+ eigenvectors.resize(eigenvalues.size());
+ for (size_t i = 0; i < eigenvalues.size(); ++i) {
+ eigenvectors[i] = SmallVector<double>{eigenvalues.size()};
+ for (size_t j = 0; j < eigenvalues.size(); ++j) {
+ eigenvectors[i][j] = eigen_solver.eigenvectors().coeff(j, i);
+ }
+ }
+ }
+
+ template <typename Eigen3MatrixEmbedderType>
+ static void
+ eigen3ComputeForSymmetricMatrix(const Eigen3MatrixEmbedderType& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P)
+ {
+ using Eigen3MatrixType = typename Eigen3MatrixEmbedderType::Eigen3MatrixType;
+ Eigen::SelfAdjointEigenSolver<Eigen3MatrixType> eigen_solver;
+
+ eigen_solver.compute(A.matrix(), Eigen::ComputeEigenvectors);
+
+ eigenvalues = SmallArray<double>(A.numberOfRows());
+ for (size_t i = 0; i < eigenvalues.size(); ++i) {
+ eigenvalues[i] = eigen_solver.eigenvalues()[i];
+ }
+
+ P = SmallMatrix<double>(eigenvalues.size(), eigenvalues.size());
+ for (size_t i = 0; i < eigenvalues.size(); ++i) {
+ for (size_t j = 0; j < eigenvalues.size(); ++j) {
+ P(i, j) = eigen_solver.eigenvectors().coeff(i, j);
+ }
+ }
+ }
+#endif // PUGS_HAS_EIGEN3
};
#ifdef PUGS_HAS_SLEPC
@@ -235,6 +306,86 @@ EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<doub
#endif // PUGS_HAS_SLEPC
+#ifdef PUGS_HAS_EIGEN3
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues)
+{
+ Eigen3SparseMatrixEmbedder embedded_A{A};
+ Internals::eigen3ComputeForSymmetricMatrix(embedded_A, eigenvalues);
+}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues)
+{
+ Internals::eigen3ComputeForSymmetricMatrix(Eigen3DenseMatrixEmbedder{A}, eigenvalues);
+}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors)
+{
+ Internals::eigen3ComputeForSymmetricMatrix(Eigen3SparseMatrixEmbedder{A}, eigenvalues, eigenvectors);
+}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors)
+{
+ Internals::eigen3ComputeForSymmetricMatrix(Eigen3DenseMatrixEmbedder{A}, eigenvalues, eigenvectors);
+}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P)
+{
+ Internals::eigen3ComputeForSymmetricMatrix(Eigen3SparseMatrixEmbedder{A}, eigenvalues, P);
+}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P)
+{
+ Internals::eigen3ComputeForSymmetricMatrix(Eigen3DenseMatrixEmbedder{A}, eigenvalues, P);
+}
+
+#else // PUGS_HAS_EIGEN3
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&, SmallArray<double>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>&,
+ SmallArray<double>&,
+ std::vector<SmallVector<double>>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&,
+ SmallArray<double>&,
+ std::vector<SmallVector<double>>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&, SmallMatrix<double>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&,
+ SmallArray<double>&,
+ SmallMatrix<double>&)
+{}
+
+#endif // PUGS_HAS_EIGEN3
+
bool
EigenvalueSolver::hasLibrary(const ESLibrary library)
{
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index eebd536e4417ee953ce8cc6317a9d5a53d07bce6..6567adcdc12caabae69baede83e74166696baaf0 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -38,6 +38,26 @@ class EigenvalueSolver
SmallArray<double>& eigenvalues,
SmallMatrix<double>& P);
+ void _eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues);
+
+ void _eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues);
+
+ void _eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors);
+
+ void _eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors);
+
+ void _eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P);
+
+ void _eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P);
+
public:
static bool hasLibrary(ESLibrary library);
static void checkHasLibrary(const ESLibrary library);
@@ -46,11 +66,12 @@ class EigenvalueSolver
void
computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues)
{
- Assert((eigenvalues.size() - A.numberOfRows() == 0) and A.isSquare());
+ Assert(A.isSquare());
switch (m_options.library()) {
case ESLibrary::eigen3: {
- throw NotImplementedError("Eigen3");
+ this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues);
+ break;
}
case ESLibrary::slepsc: {
this->_slepscComputeForSymmetricMatrix(A, eigenvalues);
@@ -68,12 +89,12 @@ class EigenvalueSolver
SmallArray<double>& eigenvalues,
std::vector<SmallVector<double>>& eigenvectors)
{
- Assert((eigenvalues.size() - A.numberOfRows() == 0) and (A.numberOfRows() - eigenvectors.size() == 0) and
- A.isSquare());
+ Assert(A.isSquare());
switch (m_options.library()) {
case ESLibrary::eigen3: {
- throw NotImplementedError("Eigen3");
+ this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
+ break;
}
case ESLibrary::slepsc: {
this->_slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
@@ -89,12 +110,12 @@ class EigenvalueSolver
void
computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues, SmallMatrix<double>& P)
{
- Assert((eigenvalues.size() - A.numberOfRows() == 0) and (A.numberOfRows() - P.numberOfRows() == 0) and
- A.isSquare() and P.isSquare());
+ Assert(A.isSquare());
switch (m_options.library()) {
case ESLibrary::eigen3: {
- throw NotImplementedError("Eigen3");
+ this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues, P);
+ break;
}
case ESLibrary::slepsc: {
this->_slepscComputeForSymmetricMatrix(A, eigenvalues, P);
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 4bb1d339d3004694ef646c4914b343b55c591384..325f3e39b0a6f34e028f99372956f67a226e39cc 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -12,14 +12,14 @@
TEST_CASE("EigenvalueSolver", "[algebra]")
{
- SECTION("SLEPc")
+ SECTION("symmetric system")
{
- EigenvalueSolverOptions options;
- options.library() = ESLibrary::slepsc;
-
- 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);
@@ -122,11 +122,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
#endif // PUGS_HAS_SLEPC
}
}
- }
- SECTION("SmallMatrix")
- {
- SECTION("symmetric system")
+ SECTION("SmallMatrix")
{
SmallMatrix<double> A{3, 3};
A.fill(0);
@@ -225,11 +222,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
#endif // PUGS_HAS_SLEPC
}
}
- }
- SECTION("TinyMatrix")
- {
- SECTION("symmetric system")
+ SECTION("TinyMatrix")
{
TinyMatrix<3> A = zero;
@@ -329,5 +323,315 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
}
}
}
+
+ SECTION("Eigen3")
+ {
+ EigenvalueSolverOptions options;
+ options.library() = ESLibrary::eigen3;
+
+ SECTION("Sparse Matrices")
+ {
+ 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(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.getCRSMatrix()};
+
+ SECTION("eigenvalues")
+ {
+ SmallArray<double> eigenvalues;
+
+#ifdef PUGS_HAS_EIGEN3
+ 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_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and eigenvectors")
+ {
+ SmallArray<double> eigenvalue_list;
+ std::vector<SmallVector<double>> eigenvector_list;
+
+#ifdef PUGS_HAS_EIGEN3
+ 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_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and transition matrix")
+ {
+ SmallArray<double> eigenvalue_list;
+ SmallMatrix<double> P{};
+
+#ifdef PUGS_HAS_EIGEN3
+ 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_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+ }
+
+ SECTION("SmallMatrix")
+ {
+ 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_EIGEN3
+ 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_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and eigenvectors")
+ {
+ SmallArray<double> eigenvalue_list;
+ std::vector<SmallVector<double>> eigenvector_list;
+
+#ifdef PUGS_HAS_EIGEN3
+ 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_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and transition matrix")
+ {
+ SmallArray<double> eigenvalue_list;
+ SmallMatrix<double> P{};
+
+#ifdef PUGS_HAS_EIGEN3
+ 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_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+ }
+
+ SECTION("TinyMatrix")
+ {
+ TinyMatrix<3> A = zero;
+
+ 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_EIGEN3
+ 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_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and eigenvectors")
+ {
+ SmallArray<double> eigenvalue_list;
+ std::vector<SmallVector<double>> eigenvector_list;
+
+#ifdef PUGS_HAS_EIGEN3
+ 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_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and transition matrix")
+ {
+ SmallArray<double> eigenvalue_list;
+ SmallMatrix<double> P{};
+
+#ifdef PUGS_HAS_EIGEN3
+ 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_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+ }
+ }
}
}