From 1dcf50f0e6d6c8a9dc8a5de1590f5237b05c866b Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Thu, 30 Jan 2025 08:07:22 +0100
Subject: [PATCH] Prepare EigenvalueSolver to incoming use of Eigen3

---
 src/algebra/CMakeLists.txt                    |   1 +
 src/algebra/EigenvalueSolver.cpp              | 250 +++++++++----
 src/algebra/EigenvalueSolver.hpp              | 109 ++++--
 src/algebra/EigenvalueSolverOptions.cpp       |  10 +
 src/algebra/EigenvalueSolverOptions.hpp       | 108 ++++++
 src/language/modules/LinearSolverModule.cpp   |  40 +-
 src/utils/checkpointing/Checkpoint.cpp        |  11 +-
 .../EigenvalueSolverOptionsHFType.hpp         |  15 +
 .../LinearSolverOptionsHFType.hpp             |   2 +-
 src/utils/checkpointing/Resume.cpp            |  10 +
 tests/test_EigenvalueSolver.cpp               | 344 ++++++++++++++----
 tests/test_LinearSolverOptions.cpp            |   1 +
 tests/test_checkpointing_HFTypes.cpp          |   3 +
 13 files changed, 737 insertions(+), 167 deletions(-)
 create mode 100644 src/algebra/EigenvalueSolverOptions.cpp
 create mode 100644 src/algebra/EigenvalueSolverOptions.hpp
 create mode 100644 src/utils/checkpointing/EigenvalueSolverOptionsHFType.hpp

diff --git a/src/algebra/CMakeLists.txt b/src/algebra/CMakeLists.txt
index d215b09b5..78c041c5c 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 c71be58dd..5888d1c0b 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 a9b462509..eebd536e4 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 000000000..2a4ef9549
--- /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 000000000..d9be86be0
--- /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 d3e2de0fa..1ff8be36d 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 5976cb2fb..ed56dcb83 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 000000000..32995c1f5
--- /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 94f88fc60..1c94d74a9 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 9d58e1d9d..effd6b808 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 fb7bf0790..4bb1d339d 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 620959714..7d078a7fb 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 41fc6948c..1a5f02b4a 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);
-- 
GitLab