diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2e0b299005e12999c0dcdcffe495f2eb69bf3662..52579983f4dc7deaa42079eae59532695cdcf9ae 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -245,6 +245,26 @@ else()
   endif()
 endif()
 
+# -----------------------------------------------------
+# Check for Eigen3
+# defaults use Eigen3
+set(PUGS_ENABLE_EIGEN3 AUTO CACHE STRING
+  "Choose one of: AUTO ON OFF")
+
+if (PUGS_ENABLE_EIGEN3 MATCHES "^(AUTO|ON)$")
+  find_package (Eigen3 NO_MODULE)
+
+  if (TARGET Eigen3::Eigen)
+    set(EIGEN3_TARGET Eigen3::Eigen)
+    set(EIGEN3_FOUND TRUE)
+  else()
+    set(EIGEN3_FOUND FALSE)
+  endif (TARGET Eigen3::Eigen)
+  set(PUGS_HAS_EIGEN3 ${EIGEN3_FOUND})
+else()
+  unset(PUGS_HAS_EIGEN3)
+endif()
+
 # -----------------------------------------------------
 
 if (${MPI_FOUND})
@@ -829,6 +849,16 @@ else()
   endif()
 endif()
 
+if (EIGEN3_FOUND)
+  message(" Eigen3: ${Eigen3_VERSION}")
+else()
+  if (PUGS_ENABLE_EIGEN3 MATCHES "^(AUTO|ON)$")
+    message(" Eigen3: not found!")
+  else()
+      message(" Eigen3: explicitly deactivated!")
+  endif()
+endif()
+
 if (HDF5_FOUND)
   message(" HDF5: ${HDF5_VERSION} parallel: ${HDF5_IS_PARALLEL}")
 else()
diff --git a/README.md b/README.md
index c859bc694043fc05a83ed47dae3f515eca0bbb42..473b27b07ed129357a206fa59c4af7d74471de52 100644
--- a/README.md
+++ b/README.md
@@ -82,6 +82,15 @@ To install `SLEPc` on Debian-like systems
 apt install slepc-dev
 ```
 
+#### `Eigen3`
+
+`Eigen3` is linear system solver and an eigenvalue problem solver.
+
+To install `Eigen3` on Debian-like systems
+```shell
+apt install libeigen3-dev
+```
+
 ## Documentation
 
 ### User documentation
@@ -234,6 +243,7 @@ the way `pugs` is built.
 | Description     | Variable            | Values                       |
 |:----------------|:--------------------|:----------------------------:|
 | MPI support     | `PUGS_ENABLE_MPI`   | `AUTO`(default), `ON`, `OFF` |
+| `Eigen3` support| `PUGS_ENABLE_EIGEN3`| `AUTO`(default), `ON`, `OFF` |
 | `PETSc` support | `PUGS_ENABLE_PETSC` | `AUTO`(default), `ON`, `OFF` |
 | `SLEPc` support | `PUGS_ENABLE_SLEPC` | `AUTO`(default), `ON`, `OFF` |
 
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/DenseMatrixWrapper.hpp b/src/algebra/DenseMatrixWrapper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dee57f27eff3d4c906cb256754d86ca6eeafbdfa
--- /dev/null
+++ b/src/algebra/DenseMatrixWrapper.hpp
@@ -0,0 +1,59 @@
+#ifndef DENSE_MATRIX_WRAPPER_HPP
+#define DENSE_MATRIX_WRAPPER_HPP
+
+#include <algebra/SmallMatrix.hpp>
+#include <algebra/TinyMatrix.hpp>
+
+template <typename DataType>
+class DenseMatrixWrapper
+{
+ private:
+  const size_t m_number_of_rows;
+  const size_t m_number_of_columns;
+  const DataType* const m_matrix_ptr;
+
+ public:
+  PUGS_INLINE
+  bool
+  isSquare() const
+  {
+    return (m_number_of_rows == m_number_of_columns);
+  }
+
+  PUGS_INLINE
+  size_t
+  numberOfRows() const
+  {
+    return m_number_of_rows;
+  }
+
+  PUGS_INLINE
+  size_t
+  numberOfColumns() const
+  {
+    return m_number_of_columns;
+  }
+
+  PUGS_INLINE
+  const DataType* const&
+  ptr() const
+  {
+    return m_matrix_ptr;
+  }
+
+  DenseMatrixWrapper(const SmallMatrix<DataType>& A)
+    : m_number_of_rows{A.numberOfRows()},
+      m_number_of_columns{A.numberOfColumns()},
+      m_matrix_ptr{(m_number_of_columns * m_number_of_rows > 0) ? (&A(0, 0)) : nullptr}
+  {}
+
+  template <size_t M, size_t N>
+  DenseMatrixWrapper(const TinyMatrix<M, N, DataType>& A)
+    : m_number_of_rows{M}, m_number_of_columns{N}, m_matrix_ptr{(M * N > 0) ? (&A(0, 0)) : nullptr}
+  {}
+
+  DenseMatrixWrapper(const DenseMatrixWrapper&) = delete;
+  DenseMatrixWrapper(DenseMatrixWrapper&&)      = delete;
+};
+
+#endif   // DENSE_MATRIX_WRAPPER_HPP
diff --git a/src/algebra/Eigen3Utils.hpp b/src/algebra/Eigen3Utils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d9d61e02746c77e22002fb1c1ecfb93db6c1a9bb
--- /dev/null
+++ b/src/algebra/Eigen3Utils.hpp
@@ -0,0 +1,145 @@
+#ifndef EIGEN3_UTILS_HPP
+#define EIGEN3_UTILS_HPP
+
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_EIGEN3
+
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/DenseMatrixWrapper.hpp>
+
+#include <eigen3/Eigen/Eigen>
+
+class Eigen3DenseMatrixEmbedder
+{
+ public:
+  using Eigen3MatrixType    = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+  using Eigen3MapMatrixType = Eigen::Map<const Eigen3MatrixType>;
+
+ private:
+  Eigen3MapMatrixType m_eigen_matrix;
+
+  Eigen3DenseMatrixEmbedder(const size_t nb_rows, const size_t nb_columns, const double* A)
+    : m_eigen_matrix{A, int(nb_rows), int(nb_columns)}
+  {}
+
+ public:
+  PUGS_INLINE
+  size_t
+  isSquare() const
+  {
+    return (m_eigen_matrix.rows() == m_eigen_matrix.cols());
+  }
+
+  PUGS_INLINE
+  size_t
+  numberOfRows() const
+  {
+    return m_eigen_matrix.rows();
+  }
+
+  PUGS_INLINE
+  size_t
+  numberOfColumns() const
+  {
+    return m_eigen_matrix.cols();
+  }
+
+  PUGS_INLINE
+  Eigen3MapMatrixType&
+  matrix()
+  {
+    return m_eigen_matrix;
+  }
+
+  PUGS_INLINE
+  const Eigen3MapMatrixType&
+  matrix() const
+  {
+    return m_eigen_matrix;
+  }
+
+  Eigen3DenseMatrixEmbedder(const DenseMatrixWrapper<double>& A)
+    : Eigen3DenseMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), A.ptr()}
+  {}
+
+  ~Eigen3DenseMatrixEmbedder() = default;
+};
+
+class Eigen3SparseMatrixEmbedder
+{
+ public:
+  using Eigen3MatrixType    = Eigen::SparseMatrix<double, Eigen::RowMajor>;
+  using Eigen3MapMatrixType = Eigen::Map<const Eigen3MatrixType>;
+
+ private:
+  Eigen3MapMatrixType m_eigen_matrix;
+
+ public:
+  PUGS_INLINE
+  size_t
+  numberOfRows() const
+  {
+    return m_eigen_matrix.rows();
+  }
+
+  PUGS_INLINE
+  size_t
+  numberOfColumns() const
+  {
+    return m_eigen_matrix.cols();
+  }
+
+  PUGS_INLINE
+  size_t
+  isSquare() const
+  {
+    return (m_eigen_matrix.rows() == m_eigen_matrix.cols());
+  }
+
+  PUGS_INLINE
+  Eigen3MapMatrixType&
+  matrix()
+  {
+    return m_eigen_matrix;
+  }
+
+  PUGS_INLINE
+  const Eigen3MapMatrixType&
+  matrix() const
+  {
+    return m_eigen_matrix;
+  }
+
+  Eigen3SparseMatrixEmbedder(const CRSMatrix<double>& A)
+    : m_eigen_matrix{A.numberOfRows(),
+                     A.numberOfColumns(),
+                     static_cast<int>(A.values().size()),                                  //
+                     (A.rowMap().size() > 0) ? &(A.rowMap()[0]) : nullptr,                 //
+                     (A.columnIndices().size() > 0) ? &(A.columnIndices()[0]) : nullptr,   //
+                     (A.values().size() > 0) ? &(A.values()[0]) : nullptr}
+  {}
+  ~Eigen3SparseMatrixEmbedder() = default;
+};
+
+#else   // PUGS_HAS_EIGEN3
+
+class Eigen3DenseMatrixEmbedder
+{
+ public:
+  Eigen3DenseMatrixEmbedder(const DenseMatrixWrapper<double>&) {}
+
+  ~Eigen3DenseMatrixEmbedder() = default;
+};
+
+class Eigen3SparseMatrixEmbedder
+{
+ public:
+  Eigen3SparseMatrixEmbedder(const CRSMatrix<double>&) {}
+
+  ~Eigen3SparseMatrixEmbedder() = default;
+};
+
+#endif   // PUGS_HAS_EIGEN3
+
+#endif   // EIGEN3_UTILS_HPP
diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index c71be58ddcb1c37498a172f30a1eb39b89f43aab..b17e650af775530d3e61f4561b7676dccf1efdd4 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -2,10 +2,17 @@
 #include <utils/pugs_config.hpp>
 
 #ifdef PUGS_HAS_SLEPC
+#include <algebra/PETScUtils.hpp>
 #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
   static PetscReal
   computeSmallestRealEigenvalueOfSymmetricMatrix(EPS& eps)
   {
@@ -65,93 +72,357 @@ struct EigenvalueSolver::Internals
     computeAllEigenvaluesOfSymmetricMatrixInInterval(eps, left_bound - 0.01 * std::abs(left_bound),
                                                      right_bound + 0.01 * std::abs(right_bound));
   }
+
+  static void
+  slepscComputeForSymmetricMatrix(const PETScAijMatrixEmbedder& 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);
+    for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
+      EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, nullptr, nullptr);
+    }
+
+    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_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
 void
-EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues)
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues)
 {
-  EPS eps;
-
-  EPSCreate(PETSC_COMM_SELF, &eps);
-  EPSSetOperators(eps, A, nullptr);
-  EPSSetProblemType(eps, EPS_HEP);
+  Internals::slepscComputeForSymmetricMatrix(A, eigenvalues);
+}
 
-  Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues)
+{
+  Internals::slepscComputeForSymmetricMatrix(A, eigenvalues);
+}
 
-  PetscInt nb_eigenvalues;
-  EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+                                                   SmallArray<double>& eigenvalues,
+                                                   std::vector<SmallVector<double>>& eigenvectors)
+{
+  Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
+}
 
-  eigenvalues = SmallArray<double>(nb_eigenvalues);
-  for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
-    EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, nullptr, nullptr);
-  }
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+                                                   SmallArray<double>& eigenvalues,
+                                                   std::vector<SmallVector<double>>& eigenvectors)
+{
+  Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
+}
 
-  EPSDestroy(&eps);
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+                                                   SmallArray<double>& eigenvalues,
+                                                   SmallMatrix<double>& P)
+{
+  Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, P);
 }
 
 void
-EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
-                                            SmallArray<double>& eigenvalues,
-                                            std::vector<SmallVector<double>>& eigenvector_list)
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+                                                   SmallArray<double>& eigenvalues,
+                                                   SmallMatrix<double>& P)
 {
-  EPS eps;
+  Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, P);
+}
 
-  EPSCreate(PETSC_COMM_SELF, &eps);
-  EPSSetOperators(eps, A, nullptr);
-  EPSSetProblemType(eps, EPS_HEP);
+#else   // PUGS_HAS_SLEPC
 
-  Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&)
+{}
 
-  PetscInt nb_eigenvalues;
-  EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&, SmallArray<double>&)
+{}
 
-  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);
-  }
+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
+
+#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);
+}
 
-  EPSDestroy(&eps);
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+                                                   SmallArray<double>& eigenvalues,
+                                                   SmallMatrix<double>& P)
+{
+  Internals::eigen3ComputeForSymmetricMatrix(Eigen3SparseMatrixEmbedder{A}, eigenvalues, P);
 }
 
 void
-EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
-                                            SmallArray<double>& eigenvalues,
-                                            SmallMatrix<double>& P)
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+                                                   SmallArray<double>& eigenvalues,
+                                                   SmallMatrix<double>& P)
 {
-  EPS eps;
+  Internals::eigen3ComputeForSymmetricMatrix(Eigen3DenseMatrixEmbedder{A}, eigenvalues, P);
+}
 
-  EPSCreate(PETSC_COMM_SELF, &eps);
-  EPSSetOperators(eps, A, nullptr);
-  EPSSetProblemType(eps, EPS_HEP);
+#else   // PUGS_HAS_EIGEN3
 
-  Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&)
+{}
 
-  PetscInt nb_eigenvalues;
-  EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&, SmallArray<double>&)
+{}
 
-  eigenvalues = SmallArray<double>(nb_eigenvalues);
-  P           = SmallMatrix<double>(nb_eigenvalues, nb_eigenvalues);
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>&,
+                                                   SmallArray<double>&,
+                                                   std::vector<SmallVector<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::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&,
+                                                   SmallArray<double>&,
+                                                   std::vector<SmallVector<double>>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&, SmallMatrix<double>&)
+{}
 
-  EPSDestroy(&eps);
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&,
+                                                   SmallArray<double>&,
+                                                   SmallMatrix<double>&)
+{}
+
+#endif   // PUGS_HAS_EIGEN3
+
+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
+  }
 }
 
-#endif   // PUGS_HAS_SLEPC
+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() {}
+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..6567adcdc12caabae69baede83e74166696baaf0 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,118 @@ 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);
+
+  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);
+
   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(A.isSquare());
+
+    switch (m_options.library()) {
+    case ESLibrary::eigen3: {
+      this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues);
+      break;
+    }
+    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(A.isSquare());
+
+    switch (m_options.library()) {
+    case ESLibrary::eigen3: {
+      this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
+      break;
+    }
+    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(A.isSquare());
+
+    switch (m_options.library()) {
+    case ESLibrary::eigen3: {
+      this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues, P);
+      break;
+    }
+    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..1c047ef3cb63df23cfb7f8725edf68ffa8ad26a6
--- /dev/null
+++ b/src/algebra/EigenvalueSolverOptions.hpp
@@ -0,0 +1,90 @@
+#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;
+
+ 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;
+  }
+
+  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/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index 99bff927c744c9ba216e575dcc914a2bdba1d4da..52bf8fc3969d2001685c5548f31c30ce7af20dd0 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -3,8 +3,13 @@
 
 #include <algebra/BiCGStab.hpp>
 #include <algebra/CG.hpp>
+#include <algebra/Eigen3Utils.hpp>
 #include <algebra/PETScUtils.hpp>
 
+#ifdef PUGS_HAS_EIGEN3
+#include <eigen3/unsupported/Eigen/IterativeSolvers>
+#endif   // PUGS_HAS_EIGEN3
+
 #ifdef PUGS_HAS_PETSC
 #include <petsc.h>
 #endif   // PUGS_HAS_PETSC
@@ -18,6 +23,13 @@ struct LinearSolver::Internals
     case LSLibrary::builtin: {
       return true;
     }
+    case LSLibrary::eigen3: {
+#ifdef PUGS_HAS_EIGEN3
+      return true;
+#else
+      return false;
+#endif
+    }
     case LSLibrary::petsc: {
 #ifdef PUGS_HAS_PETSC
       return true;
@@ -57,6 +69,25 @@ struct LinearSolver::Internals
     }
   }
 
+  static void
+  checkEigen3Method(const LSMethod method)
+  {
+    switch (method) {
+    case LSMethod::cg:
+    case LSMethod::bicgstab:
+    case LSMethod::gmres:
+    case LSMethod::lu:
+    case LSMethod::cholesky: {
+      break;
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw NormalError(name(method) + " is not an Eigen3 linear solver!");
+    }
+      // LCOV_EXCL_STOP
+    }
+  }
+
   static void
   checkPETScMethod(const LSMethod method)
   {
@@ -90,6 +121,24 @@ struct LinearSolver::Internals
     }
   }
 
+  static void
+  checkEigen3Precond(const LSPrecond precond)
+  {
+    switch (precond) {
+    case LSPrecond::none:
+    case LSPrecond::diagonal:
+    case LSPrecond::incomplete_cholesky:
+    case LSPrecond::incomplete_LU: {
+      break;
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw NormalError(name(precond) + " is not an Eigen3 preconditioner!");
+    }
+      // LCOV_EXCL_STOP
+    }
+  }
+
   static void
   checkPETScPrecond(const LSPrecond precond)
   {
@@ -118,6 +167,11 @@ struct LinearSolver::Internals
       checkBuiltinPrecond(options.precond());
       break;
     }
+    case LSLibrary::eigen3: {
+      checkEigen3Method(options.method());
+      checkEigen3Precond(options.precond());
+      break;
+    }
     case LSLibrary::petsc: {
       checkPETScMethod(options.method());
       checkPETScPrecond(options.precond());
@@ -133,10 +187,10 @@ struct LinearSolver::Internals
 
   template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
   static void
-  builtinSolveLocalSystem(const MatrixType& A,
-                          SolutionVectorType& x,
-                          const RHSVectorType& b,
-                          const LinearSolverOptions& options)
+  _builtinSolveLocalSystem(const MatrixType& A,
+                           SolutionVectorType& x,
+                           const RHSVectorType& b,
+                           const LinearSolverOptions& options)
   {
     if (options.precond() != LSPrecond::none) {
       // LCOV_EXCL_START
@@ -160,6 +214,134 @@ struct LinearSolver::Internals
     }
   }
 
+#ifdef PUGS_HAS_EIGEN3
+  template <typename PreconditionerType, typename Eigen3MatrixEmbedderType>
+  static void
+  _preconditionedEigen3SolveLocalSystem(const Eigen3MatrixEmbedderType& eigen3_A,
+                                        VectorWrapper<double>& x,
+                                        const VectorWrapper<const double>& b,
+                                        const LinearSolverOptions& options)
+  {
+    Eigen::Map<Eigen::VectorX<double>> eigen3_x{x.ptr(), static_cast<int>(x.size())};
+    Eigen::Map<const Eigen::VectorX<double>> eigen3_b{b.ptr(), static_cast<int>(b.size())};
+
+    using Eigen3MatrixType = typename Eigen3MatrixEmbedderType::Eigen3MatrixType;
+
+    switch (options.method()) {
+    case LSMethod::bicgstab: {
+      Eigen::BiCGSTAB<Eigen3MatrixType, PreconditionerType> solver;
+      solver.compute(eigen3_A.matrix());
+      solver.setMaxIterations(options.maximumIteration());
+      solver.setTolerance(options.epsilon());
+      eigen3_x = solver.solve(eigen3_b);
+      break;
+    }
+    case LSMethod::cg: {
+      Eigen::ConjugateGradient<Eigen3MatrixType, Eigen::Lower | Eigen::Upper, PreconditionerType> solver;
+      solver.compute(eigen3_A.matrix());
+      solver.setMaxIterations(options.maximumIteration());
+      solver.setTolerance(options.epsilon());
+      eigen3_x = solver.solve(eigen3_b);
+      break;
+    }
+    case LSMethod::gmres: {
+      Eigen::GMRES<Eigen3MatrixType, PreconditionerType> solver;
+      solver.compute(eigen3_A.matrix());
+      solver.setMaxIterations(options.maximumIteration());
+      solver.setTolerance(options.epsilon());
+      eigen3_x = solver.solve(eigen3_b);
+      break;
+    }
+    case LSMethod::lu: {
+      if constexpr (std::is_same_v<Eigen3MatrixEmbedderType, Eigen3DenseMatrixEmbedder>) {
+        Eigen::FullPivLU<Eigen3MatrixType> solver;
+        solver.compute(eigen3_A.matrix());
+        eigen3_x = solver.solve(eigen3_b);
+      } else {
+        Eigen::SparseLU<Eigen3MatrixType> solver;
+        solver.compute(eigen3_A.matrix());
+        eigen3_x = solver.solve(eigen3_b);
+      }
+      break;
+    }
+    case LSMethod::cholesky: {
+      if constexpr (std::is_same_v<Eigen3MatrixEmbedderType, Eigen3DenseMatrixEmbedder>) {
+        Eigen::LLT<Eigen3MatrixType> solver;
+        solver.compute(eigen3_A.matrix());
+        eigen3_x = solver.solve(eigen3_b);
+      } else {
+        Eigen::SimplicialLLT<Eigen3MatrixType> solver;
+        solver.compute(eigen3_A.matrix());
+        eigen3_x = solver.solve(eigen3_b);
+      }
+      break;
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw UnexpectedError("unexpected method: " + name(options.method()));
+    }
+      // LCOV_EXCL_STOP
+    }
+  }
+
+  template <typename Eigen3MatrixEmbedderType>
+  static void
+  _eigen3SolveLocalSystem(const Eigen3MatrixEmbedderType& eigen3_A,
+                          VectorWrapper<double>& x,
+                          const VectorWrapper<const double>& b,
+                          const LinearSolverOptions& options)
+  {
+    switch (options.precond()) {
+    case LSPrecond::none: {
+      _preconditionedEigen3SolveLocalSystem<Eigen::IdentityPreconditioner, Eigen3MatrixEmbedderType>(eigen3_A, x, b,
+                                                                                                     options);
+      break;
+    }
+    case LSPrecond::diagonal: {
+      _preconditionedEigen3SolveLocalSystem<Eigen::DiagonalPreconditioner<double>, Eigen3MatrixEmbedderType>(eigen3_A,
+                                                                                                             x, b,
+                                                                                                             options);
+      break;
+    }
+    case LSPrecond::incomplete_cholesky: {
+      if constexpr (std::is_same_v<Eigen3SparseMatrixEmbedder, Eigen3MatrixEmbedderType>) {
+        _preconditionedEigen3SolveLocalSystem<Eigen::IncompleteCholesky<double>, Eigen3MatrixEmbedderType>(eigen3_A, x,
+                                                                                                           b, options);
+      } else {
+        throw NormalError("incomplete cholesky is not available for dense matrices in Eigen3");
+      }
+      break;
+    }
+    case LSPrecond::incomplete_LU: {
+      if constexpr (std::is_same_v<Eigen3SparseMatrixEmbedder, Eigen3MatrixEmbedderType>) {
+        _preconditionedEigen3SolveLocalSystem<Eigen::IncompleteLUT<double>, Eigen3MatrixEmbedderType>(eigen3_A, x, b,
+                                                                                                      options);
+      } else {
+        throw NormalError("incomplete LU is not available for dense matrices in Eigen3");
+      }
+      break;
+    }
+      // LCOV_EXCL_START
+    default: {
+      throw UnexpectedError("unexpected preconditioner: " + name(options.precond()));
+    }
+      // LCOV_EXCL_STOP
+    }
+  }
+#else   // PUGS_HAS_EIGEN3
+
+  // LCOV_EXCL_START
+  template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
+  static void
+  _eigen3SolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&)
+  {
+    checkHasLibrary(LSLibrary::eigen3);
+    throw UnexpectedError("unexpected situation should not reach this point!");
+  }
+  // LCOV_EXCL_STOP
+
+#endif   // PUGS_HAS_EIGEN3
+
 #ifdef PUGS_HAS_PETSC
   static int
   petscMonitor(KSP, int i, double residu, void*)
@@ -168,19 +350,17 @@ struct LinearSolver::Internals
     return 0;
   }
 
-  template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
+  template <typename MatrixType>
   static void
-  petscSolveLocalSystem(const MatrixType& A,
-                        SolutionVectorType& x,
-                        const RHSVectorType& b,
-                        const LinearSolverOptions& options)
+  _petscSolveLocalSystem(const MatrixType& A,
+                         VectorWrapper<double>& x,
+                         const VectorWrapper<const double>& b,
+                         const LinearSolverOptions& options)
   {
-    Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
-
     Vec petscB;
-    VecCreateMPIWithArray(PETSC_COMM_SELF, 1, b.size(), b.size(), &b[0], &petscB);
+    VecCreateMPIWithArray(PETSC_COMM_SELF, 1, b.size(), b.size(), b.ptr(), &petscB);
     Vec petscX;
-    VecCreateMPIWithArray(PETSC_COMM_SELF, 1, x.size(), x.size(), &x[0], &petscX);
+    VecCreateMPIWithArray(PETSC_COMM_SELF, 1, x.size(), x.size(), x.ptr(), &petscX);
 
     PETScAijMatrixEmbedder petscA(A);
 
@@ -280,7 +460,7 @@ struct LinearSolver::Internals
   // LCOV_EXCL_START
   template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
   static void
-  petscSolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&)
+  _petscSolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&)
   {
     checkHasLibrary(LSLibrary::petsc);
     throw UnexpectedError("unexpected situation should not reach this point!");
@@ -288,33 +468,6 @@ struct LinearSolver::Internals
   // LCOV_EXCL_STOP
 
 #endif   // PUGS_HAS_PETSC
-
-  template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
-  static void
-  solveLocalSystem(const MatrixType& A,
-                   SolutionVectorType& x,
-                   const RHSVectorType& b,
-                   const LinearSolverOptions& options)
-  {
-    switch (options.library()) {
-    case LSLibrary::builtin: {
-      builtinSolveLocalSystem(A, x, b, options);
-      break;
-    }
-      // LCOV_EXCL_START
-    case LSLibrary::petsc: {
-      // not covered since if PETSc is not linked this point is
-      // unreachable: LinearSolver throws an exception at construction
-      // in this case.
-      petscSolveLocalSystem(A, x, b, options);
-      break;
-    }
-    default: {
-      throw UnexpectedError(::name(options.library()) + " cannot solve local systems");
-    }
-      // LCOV_EXCL_STOP
-    }
-  }
 };
 
 bool
@@ -330,15 +483,52 @@ LinearSolver::checkOptions(const LinearSolverOptions& options) const
 }
 
 void
-LinearSolver::solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b)
+LinearSolver::_builtinSolveLocalSystem(const CRSMatrix<double, int>& A,
+                                       Vector<double>& x,
+                                       const Vector<const double>& b)
+{
+  Internals::_builtinSolveLocalSystem(A, x, b, m_options);
+}
+
+void
+LinearSolver::_builtinSolveLocalSystem(const SmallMatrix<double>& A,
+                                       SmallVector<double>& x,
+                                       const SmallVector<const double>& b)
+{
+  Internals::_builtinSolveLocalSystem(A, x, b, m_options);
+}
+
+void
+LinearSolver::_eigen3SolveLocalSystem(const Eigen3SparseMatrixEmbedder& A,
+                                      VectorWrapper<double> x,
+                                      const VectorWrapper<const double>& b)
+{
+  Internals::_eigen3SolveLocalSystem(A, x, b, m_options);
+}
+
+void
+LinearSolver::_eigen3SolveLocalSystem(const DenseMatrixWrapper<double>& A,
+                                      VectorWrapper<double> x,
+                                      const VectorWrapper<const double>& b)
+{
+  Eigen3DenseMatrixEmbedder A_wrapper{A};
+  Internals::_eigen3SolveLocalSystem(A_wrapper, x, b, m_options);
+}
+
+void
+LinearSolver::_petscSolveLocalSystem(const CRSMatrix<double, int>& A,
+                                     VectorWrapper<double> x,
+                                     const VectorWrapper<const double>& b)
 {
-  Internals::solveLocalSystem(A, x, b, m_options);
+  Internals::_petscSolveLocalSystem(A, x, b, m_options);
 }
 
 void
-LinearSolver::solveLocalSystem(const SmallMatrix<double>& A, SmallVector<double>& x, const SmallVector<const double>& b)
+LinearSolver::_petscSolveLocalSystem(const DenseMatrixWrapper<double>& A,
+                                     VectorWrapper<double> x,
+                                     const VectorWrapper<const double>& b)
 {
-  Internals::solveLocalSystem(A, x, b, m_options);
+  Internals::_petscSolveLocalSystem(A, x, b, m_options);
 }
 
 LinearSolver::LinearSolver(const LinearSolverOptions& options) : m_options{options}
diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
index e90402c8ddc78c9426f97300ab6b10730b4c503b..2d03a9cf5e7cced7f99b6707c042395f48054199 100644
--- a/src/algebra/LinearSolver.hpp
+++ b/src/algebra/LinearSolver.hpp
@@ -2,11 +2,14 @@
 #define LINEAR_SOLVER_HPP
 
 #include <algebra/CRSMatrix.hpp>
+#include <algebra/DenseMatrixWrapper.hpp>
+#include <algebra/Eigen3Utils.hpp>
 #include <algebra/LinearSolverOptions.hpp>
 #include <algebra/SmallMatrix.hpp>
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
 #include <algebra/Vector.hpp>
+#include <algebra/VectorWrapper.hpp>
 #include <utils/Exceptions.hpp>
 
 class LinearSolver
@@ -16,29 +19,139 @@ class LinearSolver
 
   const LinearSolverOptions m_options;
 
+  void _builtinSolveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b);
+
+  void _builtinSolveLocalSystem(const SmallMatrix<double>& A,
+                                SmallVector<double>& x,
+                                const SmallVector<const double>& b);
+
+  void _eigen3SolveLocalSystem(const Eigen3SparseMatrixEmbedder& A,
+                               VectorWrapper<double> x,
+                               const VectorWrapper<const double>& b);
+
+  void _eigen3SolveLocalSystem(const DenseMatrixWrapper<double>& A,
+                               VectorWrapper<double> x,
+                               const VectorWrapper<const double>& b);
+
+  void _petscSolveLocalSystem(const CRSMatrix<double, int>& A,
+                              VectorWrapper<double> x,
+                              const VectorWrapper<const double>& b);
+
+  void _petscSolveLocalSystem(const DenseMatrixWrapper<double>& A,
+                              VectorWrapper<double> x,
+                              const VectorWrapper<const double>& b);
+
  public:
   bool hasLibrary(LSLibrary library) const;
   void checkOptions(const LinearSolverOptions& options) const;
 
+  void
+  solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b)
+  {
+    Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
+    switch (m_options.library()) {
+    case LSLibrary::builtin: {
+      _builtinSolveLocalSystem(A, x, b);
+      break;
+    }
+      // LCOV_EXCL_START
+    case LSLibrary::eigen3: {
+      // not covered since if Eigen3 is not linked this point is
+      // unreachable: LinearSolver throws an exception at construction
+      // in this case.
+      _eigen3SolveLocalSystem(A, x, b);
+      break;
+    }
+    case LSLibrary::petsc: {
+      // not covered since if PETSc is not linked this point is
+      // unreachable: LinearSolver throws an exception at construction
+      // in this case.
+      _petscSolveLocalSystem(A, x, b);
+      break;
+    }
+    default: {
+      throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems");
+    }
+      // LCOV_EXCL_STOP
+    }
+  }
+
+  void
+  solveLocalSystem(const SmallMatrix<double>& A, SmallVector<double>& x, const SmallVector<const double>& b)
+  {
+    Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
+    switch (m_options.library()) {
+    case LSLibrary::builtin: {
+      _builtinSolveLocalSystem(A, x, b);
+      break;
+    }
+      // LCOV_EXCL_START
+    case LSLibrary::eigen3: {
+      // not covered since if Eigen3 is not linked this point is
+      // unreachable: LinearSolver throws an exception at construction
+      // in this case.
+      _eigen3SolveLocalSystem(A, x, b);
+      break;
+    }
+    case LSLibrary::petsc: {
+      // not covered since if PETSc is not linked this point is
+      // unreachable: LinearSolver throws an exception at construction
+      // in this case.
+      _petscSolveLocalSystem(A, x, b);
+      break;
+    }
+    default: {
+      throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems");
+    }
+      // LCOV_EXCL_STOP
+    }
+  }
+
   template <size_t N>
   void
   solveLocalSystem(const TinyMatrix<N>& A, TinyVector<N>& x, const TinyVector<N>& b)
   {
-    SmallMatrix A_dense{A};
+    Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
+    switch (m_options.library()) {
+    case LSLibrary::builtin: {
+      // Probably expensive but builtin is not the preferred way to
+      // solve linear systems
+      SmallMatrix A_dense{A};
+
+      SmallVector x_vector{x};
+      SmallVector<const double> b_vector{b};
 
-    SmallVector x_vector{x};
-    SmallVector b_vector{b};
+      this->solveLocalSystem(A_dense, x_vector, b_vector);
 
-    this->solveLocalSystem(A_dense, x_vector, b_vector);
+      for (size_t i = 0; i < N; ++i) {
+        x[i] = x_vector[i];
+      }
 
-    for (size_t i = 0; i < N; ++i) {
-      x[i] = x_vector[i];
+      _builtinSolveLocalSystem(A, x, b);
+      break;
+    }
+      // LCOV_EXCL_START
+    case LSLibrary::eigen3: {
+      // not covered since if Eigen3 is not linked this point is
+      // unreachable: LinearSolver throws an exception at construction
+      // in this case.
+      _eigen3SolveLocalSystem(A, x, b);
+      break;
+    }
+    case LSLibrary::petsc: {
+      // not covered since if PETSc is not linked this point is
+      // unreachable: LinearSolver throws an exception at construction
+      // in this case.
+      _petscSolveLocalSystem(A, x, b);
+      break;
+    }
+    default: {
+      throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems");
+    }
+      // LCOV_EXCL_STOP
     }
   }
 
-  void solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b);
-  void solveLocalSystem(const SmallMatrix<double>& A, SmallVector<double>& x, const SmallVector<const double>& b);
-
   LinearSolver();
   LinearSolver(const LinearSolverOptions& options);
 };
diff --git a/src/algebra/LinearSolverOptions.hpp b/src/algebra/LinearSolverOptions.hpp
index 940ce1ed63a2dbbf046022aee838c3f9438c1906..97a1748780f19d28e98f299fd7679d7ab15466f0 100644
--- a/src/algebra/LinearSolverOptions.hpp
+++ b/src/algebra/LinearSolverOptions.hpp
@@ -10,6 +10,7 @@ enum class LSLibrary : int8_t
   LS__begin = 0,
   //
   builtin = LS__begin,
+  eigen3,
   petsc,
   //
   LS__end
@@ -49,6 +50,9 @@ name(const LSLibrary library)
   case LSLibrary::builtin: {
     return "builtin";
   }
+  case LSLibrary::eigen3: {
+    return "Eigen3";
+  }
   case LSLibrary::petsc: {
     return "PETSc";
   }
diff --git a/src/algebra/PETScUtils.hpp b/src/algebra/PETScUtils.hpp
index c4db9ee5910ef8e58333dbfb6dc20c8a1987bb46..c6d2e8be88c6098e2f187cda9177c120facc6f23 100644
--- a/src/algebra/PETScUtils.hpp
+++ b/src/algebra/PETScUtils.hpp
@@ -6,6 +6,7 @@
 #ifdef PUGS_HAS_PETSC
 
 #include <algebra/CRSMatrix.hpp>
+#include <algebra/DenseMatrixWrapper.hpp>
 #include <algebra/SmallMatrix.hpp>
 #include <algebra/TinyMatrix.hpp>
 
@@ -52,12 +53,8 @@ class PETScAijMatrixEmbedder
     return m_petscMat;
   }
 
-  template <size_t N>
-  PETScAijMatrixEmbedder(const TinyMatrix<N>& A) : PETScAijMatrixEmbedder{N, N, &A(0, 0)}
-  {}
-
-  PETScAijMatrixEmbedder(const SmallMatrix<double>& A)
-    : PETScAijMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), &A(0, 0)}
+  PETScAijMatrixEmbedder(const DenseMatrixWrapper<double>& A)
+    : PETScAijMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), A.ptr()}
   {}
 
   PETScAijMatrixEmbedder(const CRSMatrix<double, int>& A);
diff --git a/src/algebra/VectorWrapper.hpp b/src/algebra/VectorWrapper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0bdeefbfd3f10407358613e7f1d2487bba22b26d
--- /dev/null
+++ b/src/algebra/VectorWrapper.hpp
@@ -0,0 +1,47 @@
+#ifndef VECTOR_WRAPPER_HPP
+#define VECTOR_WRAPPER_HPP
+
+#include <algebra/SmallVector.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+
+template <typename DataType>
+class VectorWrapper
+{
+ private:
+  const size_t m_size;
+  DataType* const m_vector_ptr;
+
+ public:
+  PUGS_INLINE
+  size_t
+  size() const
+  {
+    return m_size;
+  }
+
+  PUGS_INLINE
+  DataType* const&
+  ptr() const
+  {
+    return m_vector_ptr;
+  }
+
+  VectorWrapper(const Vector<DataType>& x) : m_size{x.size()}, m_vector_ptr{(m_size > 0) ? (&x[0]) : nullptr} {}
+  VectorWrapper(Vector<DataType>& x) : m_size{x.size()}, m_vector_ptr{(m_size > 0) ? (&x[0]) : nullptr} {}
+
+  VectorWrapper(const SmallVector<DataType>& x) : m_size{x.size()}, m_vector_ptr{(m_size > 0) ? (&x[0]) : nullptr} {}
+
+  template <size_t N>
+  VectorWrapper(const TinyVector<N, DataType>& x) : m_size{N}, m_vector_ptr{(N > 0) ? (&x[0]) : nullptr}
+  {}
+
+  template <size_t N>
+  VectorWrapper(TinyVector<N, DataType>& x) : m_size{N}, m_vector_ptr{(N > 0) ? (&x[0]) : nullptr}
+  {}
+
+  VectorWrapper(const VectorWrapper&) = delete;
+  VectorWrapper(VectorWrapper&&)      = delete;
+};
+
+#endif   // VECTOR_WRAPPER_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/BuildInfo.cpp b/src/utils/BuildInfo.cpp
index 81d7567750eaef8513b580daf11eac6d3f916642..2413bdbaeb5712e3b7eaf96fd9d09df5e8f0a6e0 100644
--- a/src/utils/BuildInfo.cpp
+++ b/src/utils/BuildInfo.cpp
@@ -17,6 +17,10 @@
 #include <slepc.h>
 #endif   // PUGS_HAS_PETSC
 
+#ifdef PUGS_HAS_EIGEN3
+#include <eigen3/Eigen/Eigen>
+#endif   // PUGS_HAS_EIGEN3
+
 #ifdef PUGS_HAS_HDF5
 #include <highfive/highfive.hpp>
 #endif   // PUGS_HAS_HDF5
@@ -82,6 +86,16 @@ BuildInfo::slepcLibrary()
 #endif   // PUGS_HAS_SLEPC
 }
 
+std::string
+BuildInfo::eigen3Library()
+{
+#ifdef PUGS_HAS_EIGEN3
+  return stringify(EIGEN_WORLD_VERSION) + "." + stringify(EIGEN_MAJOR_VERSION) + "." + stringify(EIGEN_MINOR_VERSION);
+#else    // PUGS_HAS_EIGEN3
+  return "none";
+#endif   // PUGS_HAS_EIGEN3
+}
+
 std::string
 BuildInfo::hdf5Library()
 {
diff --git a/src/utils/BuildInfo.hpp b/src/utils/BuildInfo.hpp
index df4a32d93a284db7fc8d3abbea9ac2e43ea5d36e..86f7c04f035827a936a00c9b1f42e304e33e7b26 100644
--- a/src/utils/BuildInfo.hpp
+++ b/src/utils/BuildInfo.hpp
@@ -11,6 +11,7 @@ struct BuildInfo
   static std::string mpiLibrary();
   static std::string petscLibrary();
   static std::string slepcLibrary();
+  static std::string eigen3Library();
   static std::string hdf5Library();
   static std::string slurmLibrary();
 };
diff --git a/src/utils/PugsTraits.hpp b/src/utils/PugsTraits.hpp
index 2e1d409578c9295a1ad032466842c2f336ba272f..3c7a4e40a217d957dea388223d99574e6b1ecf74 100644
--- a/src/utils/PugsTraits.hpp
+++ b/src/utils/PugsTraits.hpp
@@ -15,6 +15,12 @@ class TinyVector;
 template <size_t M, size_t N, typename T>
 class TinyMatrix;
 
+template <typename DataType>
+class SmallMatrix;
+
+template <typename DataType, typename IndexType>
+class CRSMatrix;
+
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 class ItemValue;
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
@@ -114,6 +120,22 @@ inline constexpr bool is_tiny_matrix_v = false;
 template <size_t M, size_t N, typename T>
 inline constexpr bool is_tiny_matrix_v<TinyMatrix<M, N, T>> = true;
 
+// Traits is_small_matrix
+
+template <typename T>
+inline constexpr bool is_small_matrix_v = false;
+
+template <typename DataType>
+inline constexpr bool is_small_matrix_v<SmallMatrix<DataType>> = true;
+
+// Traits is_crs_matrix
+
+template <typename T>
+inline constexpr bool is_crs_matrix_v = false;
+
+template <typename DataType, typename IndexType>
+inline constexpr bool is_crs_matrix_v<CRSMatrix<DataType, IndexType>> = true;
+
 // Trais is ItemValue
 
 template <typename T>
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index 1ac38b32696417d887e91205715b9fc80d9f81af..ac64a50cfe0b026fe0185d79e30d4f6d648317ac 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -64,6 +64,7 @@ pugsBuildInfo()
   os << "MPI:      " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
   os << "PETSc:    " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n';
   os << "SLEPc:    " << rang::style::bold << BuildInfo::slepcLibrary() << rang::style::reset << '\n';
+  os << "Eigen3:   " << rang::style::bold << BuildInfo::eigen3Library() << rang::style::reset << '\n';
   os << "HDF5:     " << rang::style::bold << BuildInfo::hdf5Library() << rang::style::reset << '\n';
   os << "SLURM:    " << rang::style::bold << BuildInfo::slurmLibrary() << rang::style::reset << '\n';
   os << "-------------------------------------------------------";
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/src/utils/pugs_config.hpp.in b/src/utils/pugs_config.hpp.in
index 7402a526d5f34505dbda8524226c971e4701d4b6..a07d8c314da95f8db595886e3df5b29cdc45907a 100644
--- a/src/utils/pugs_config.hpp.in
+++ b/src/utils/pugs_config.hpp.in
@@ -5,6 +5,7 @@
 #cmakedefine PUGS_HAS_MPI
 #cmakedefine PUGS_HAS_PETSC
 #cmakedefine PUGS_HAS_SLEPC
+#cmakedefine PUGS_HAS_EIGEN3
 #cmakedefine PUGS_HAS_HDF5
 #cmakedefine PUGS_HAS_SLURM
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 23a4ac3094424515b6d488a577781cfe6bbcce5e..1dd13e6e3ce1f23c31e046db7597a1ed2ffb947d 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -104,6 +104,7 @@ add_executable (unit_tests
   test_DualMeshManager.cpp
   test_DualMeshType.cpp
   test_EdgeIntegrator.cpp
+  test_Eigen3Utils.cpp
   test_EigenvalueSolver.cpp
   test_EmbeddedData.cpp
   test_EmbeddedDiscreteFunctionUtils.cpp
diff --git a/tests/test_Eigen3Utils.cpp b/tests/test_Eigen3Utils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0c51aad3c03ee90695f837c4a9cbff0566d566f3
--- /dev/null
+++ b/tests/test_Eigen3Utils.cpp
@@ -0,0 +1,107 @@
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_EIGEN3
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <algebra/Eigen3Utils.hpp>
+
+#include <algebra/CRSMatrixDescriptor.hpp>
+
+#include <eigen3/Eigen/Core>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("Eigen3Utils", "[algebra]")
+{
+  SECTION("Eigen3DenseMatrixEmbedder")
+  {
+    SECTION("from TinyMatrix")
+    {
+      TinyMatrix<3> A{1, 2, 3, 4, 5, 6, 7, 8, 9};
+      Eigen3DenseMatrixEmbedder eigen3_A{A};
+
+      for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          REQUIRE(eigen3_A.matrix().coeff(i, j) == 1 + i * 3 + j);
+        }
+      }
+    }
+
+    SECTION("from SmallMatrix")
+    {
+      SmallMatrix<double> A(3, 3);
+      for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          A(i, j) = 1 + i * 3 + j;
+        }
+      }
+
+      Eigen3DenseMatrixEmbedder eigen3_A{A};
+
+      REQUIRE(eigen3_A.numberOfRows() == 3);
+      REQUIRE(eigen3_A.numberOfColumns() == 3);
+
+      for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          REQUIRE(eigen3_A.matrix().coeff(i, j) == 1 + i * 3 + j);
+        }
+      }
+    }
+
+    SECTION("from SmallMatrix [non-square]")
+    {
+      SmallMatrix<double> A(4, 3);
+      for (size_t i = 0; i < 4; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          A(i, j) = 1 + i * 3 + j;
+        }
+      }
+
+      Eigen3DenseMatrixEmbedder eigen3_A{A};
+
+      REQUIRE(eigen3_A.numberOfRows() == 4);
+      REQUIRE(eigen3_A.numberOfColumns() == 3);
+
+      for (size_t i = 0; i < 4; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          REQUIRE(eigen3_A.matrix().coeff(i, j) == 1 + i * 3 + j);
+        }
+      }
+    }
+  }
+
+  SECTION("from CRSMatrix")
+  {
+    Array<int> non_zeros(4);
+    non_zeros[0] = 1;
+    non_zeros[1] = 2;
+    non_zeros[2] = 3;
+    non_zeros[3] = 2;
+
+    CRSMatrixDescriptor<double, int> A(4, 3, non_zeros);
+
+    A(0, 0) = 1;
+    A(1, 0) = 2;
+    A(1, 2) = 3;
+    A(2, 0) = 4;
+    A(2, 1) = 5;
+    A(2, 2) = 6;
+    A(3, 1) = 7;
+    A(3, 2) = 8;
+
+    Eigen3SparseMatrixEmbedder eigen3_A{A.getCRSMatrix()};
+
+    REQUIRE(eigen3_A.matrix().coeff(0, 0) == 1);
+    REQUIRE(eigen3_A.matrix().coeff(1, 0) == 2);
+    REQUIRE(eigen3_A.matrix().coeff(1, 2) == 3);
+    REQUIRE(eigen3_A.matrix().coeff(2, 0) == 4);
+    REQUIRE(eigen3_A.matrix().coeff(2, 1) == 5);
+    REQUIRE(eigen3_A.matrix().coeff(2, 2) == 6);
+    REQUIRE(eigen3_A.matrix().coeff(3, 1) == 7);
+    REQUIRE(eigen3_A.matrix().coeff(3, 2) == 8);
+  }
+}
+
+#endif   // PUGS_HAS_EIGEN3
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index fb7bf07907bf6c9e1db9333c9713365350918b06..325f3e39b0a6f34e028f99372956f67a226e39cc 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -12,109 +12,625 @@
 
 TEST_CASE("EigenvalueSolver", "[algebra]")
 {
-  SECTION("Sparse Matrices")
+  SECTION("symmetric system")
   {
-    SECTION("symmetric system")
+    SECTION("SLEPc")
     {
-      Array<int> non_zeros(3);
-      non_zeros.fill(3);
-      CRSMatrixDescriptor<double> S{3, 3, non_zeros};
+      EigenvalueSolverOptions options;
+      options.library() = ESLibrary::slepsc;
 
-      S(0, 0) = 3;
-      S(0, 1) = 2;
-      S(0, 2) = 4;
+      SECTION("Sparse Matrices")
+      {
+        Array<int> non_zeros(3);
+        non_zeros.fill(3);
+        CRSMatrixDescriptor<double> S{3, 3, non_zeros};
 
-      S(1, 0) = 2;
-      S(1, 1) = 0;
-      S(1, 2) = 2;
+        S(0, 0) = 3;
+        S(0, 1) = 2;
+        S(0, 2) = 4;
 
-      S(2, 0) = 4;
-      S(2, 1) = 2;
-      S(2, 2) = 3;
+        S(1, 0) = 2;
+        S(1, 1) = 0;
+        S(1, 2) = 2;
 
-      CRSMatrix A{S.getCRSMatrix()};
+        S(2, 0) = 4;
+        S(2, 1) = 2;
+        S(2, 2) = 3;
 
-      SECTION("eigenvalues")
-      {
-        SmallArray<double> eigenvalues;
+        CRSMatrix A{S.getCRSMatrix()};
+
+        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, 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, 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{}.computeForSymmetricMatrix(A, eigenvalues),
-                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
+          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) - S(i, j) == Catch::Approx(0).margin(1E-13));
+            }
+          }
+#else    //  PUGS_HAS_SLEPC
+          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")
       {
-        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{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
+          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
+        }
 
-        REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
-        REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
-        REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+        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);
 
-        SmallMatrix<double> P{3};
-        SmallMatrix<double> PT{3};
-        SmallMatrix<double> D{3};
-        D = zero;
+          REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+          REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+          REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
 
-        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];
+          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];
           }
-          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")
       {
-        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
+        }
+      }
+    }
+
+    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
+        }
       }
     }
   }
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 11195faab3dad2f94e82dba78b9afc368d0edcd0..76ef6aa2bb75c3ea94fcf392f59ddb74350dc50a 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -21,6 +21,12 @@ TEST_CASE("LinearSolver", "[algebra]")
 #else    // PUGS_HAS_PETSC
     REQUIRE(linear_solver.hasLibrary(LSLibrary::petsc) == false);
 #endif   // PUGS_HAS_PETSC
+
+#ifdef PUGS_HAS_EIGEN3
+    REQUIRE(linear_solver.hasLibrary(LSLibrary::eigen3) == true);
+#else    // PUGS_HAS_PETSC
+    REQUIRE(linear_solver.hasLibrary(LSLibrary::eigen3) == false);
+#endif   // PUGS_HAS_PETSC
   }
 
   SECTION("check linear solver building")
@@ -108,7 +114,7 @@ TEST_CASE("LinearSolver", "[algebra]")
         REQUIRE_NOTHROW(linear_solver.checkOptions(options));
       }
 
-      SECTION("builtin precond")
+      SECTION("PETSc precond")
       {
         options.library() = LSLibrary::petsc;
         options.method()  = LSMethod::cg;
@@ -140,6 +146,66 @@ TEST_CASE("LinearSolver", "[algebra]")
       REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
     }
 #endif   // PUGS_HAS_PETSC
+
+    SECTION("Eigen3")
+    {
+      LinearSolverOptions always_valid;
+      always_valid.library() = LSLibrary::builtin;
+      always_valid.method()  = LSMethod::cg;
+      always_valid.precond() = LSPrecond::none;
+
+      LinearSolver linear_solver{always_valid};
+
+      SECTION("Eigen3 methods")
+      {
+        options.library() = LSLibrary::eigen3;
+        options.precond() = LSPrecond::none;
+
+        options.method() = LSMethod::cg;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.method() = LSMethod::bicgstab;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.method() = LSMethod::lu;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.method() = LSMethod::cholesky;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.method() = LSMethod::gmres;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+      }
+
+      SECTION("Eigen3 precond")
+      {
+        options.library() = LSLibrary::eigen3;
+        options.method()  = LSMethod::cg;
+
+        options.precond() = LSPrecond::none;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.precond() = LSPrecond::diagonal;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.precond() = LSPrecond::incomplete_LU;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+        options.precond() = LSPrecond::incomplete_cholesky;
+        REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+      }
+    }
+
+#ifndef PUGS_HAS_EIGEN3
+    SECTION("not linked Eigen3")
+    {
+      options.library() = LSLibrary::eigen3;
+      options.method()  = LSMethod::cg;
+      options.precond() = LSPrecond::none;
+
+      REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+    }
+#endif   // PUGS_HAS_EIGEN3
   }
 
   SECTION("Sparse Matrices")
@@ -205,6 +271,101 @@ TEST_CASE("LinearSolver", "[algebra]")
           }
         }
 
+        SECTION("Eigen3")
+        {
+#ifdef PUGS_HAS_EIGEN3
+
+          SECTION("CG")
+          {
+            LinearSolverOptions options;
+            options.library() = LSLibrary::eigen3;
+            options.method()  = LSMethod::cg;
+            options.precond() = LSPrecond::none;
+            options.verbose() = true;
+
+            SECTION("CG no preconditioner")
+            {
+              options.precond() = LSPrecond::none;
+
+              Vector<double> x{5};
+              x = zero;
+
+              LinearSolver solver{options};
+
+              solver.solveLocalSystem(A, x, b);
+              Vector error = x - x_exact;
+              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+            }
+
+            SECTION("CG Diagonal")
+            {
+              options.precond() = LSPrecond::diagonal;
+
+              Vector<double> x{5};
+              x = zero;
+
+              LinearSolver solver{options};
+
+              solver.solveLocalSystem(A, x, b);
+              Vector error = x - x_exact;
+              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+            }
+
+            SECTION("CG ICholesky")
+            {
+              options.precond() = LSPrecond::incomplete_cholesky;
+
+              Vector<double> x{5};
+              x = zero;
+
+              LinearSolver solver{options};
+
+              solver.solveLocalSystem(A, x, b);
+              Vector error = x - x_exact;
+              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+            }
+
+            SECTION("CG AMG")
+            {
+              options.precond() = LSPrecond::amg;
+
+              Vector<double> x{5};
+              x = zero;
+
+              REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not an Eigen3 preconditioner!");
+            }
+          }
+
+          SECTION("Cholesky")
+          {
+            LinearSolverOptions options;
+            options.library() = LSLibrary::eigen3;
+            options.method()  = LSMethod::cholesky;
+            options.precond() = LSPrecond::none;
+
+            Vector<double> x{5};
+            x = zero;
+
+            LinearSolver solver{options};
+
+            solver.solveLocalSystem(A, x, b);
+            Vector error = x - x_exact;
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+          }
+
+#else    // PUGS_HAS_EIGEN3
+          SECTION("Eigen3 not linked")
+          {
+            LinearSolverOptions options;
+            options.library() = LSLibrary::eigen3;
+            options.method()  = LSMethod::cg;
+            options.precond() = LSPrecond::none;
+
+            REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+          }
+#endif   // PUGS_HAS_EIGEN3
+        }
+
         SECTION("PETSc")
         {
 #ifdef PUGS_HAS_PETSC
@@ -535,170 +696,150 @@ TEST_CASE("LinearSolver", "[algebra]")
           }
 #endif   // PUGS_HAS_PETSC
         }
-      }
-    }
-  }
-
-  SECTION("Dense Matrices")
-  {
-    SECTION("check linear solvers")
-    {
-      SECTION("symmetric system")
-      {
-        SmallMatrix<double> A{5};
-        A = zero;
-
-        A(0, 0) = 2;
-        A(0, 1) = -1;
-
-        A(1, 0) = -1;
-        A(1, 1) = 2;
-        A(1, 2) = -1;
-
-        A(2, 1) = -1;
-        A(2, 2) = 2;
-        A(2, 3) = -1;
-
-        A(3, 2) = -1;
-        A(3, 3) = 2;
-        A(3, 4) = -1;
-
-        A(4, 3) = -1;
-        A(4, 4) = 2;
-
-        SmallVector<const double> x_exact = [] {
-          SmallVector<double> y{5};
-          y[0] = 1;
-          y[1] = 3;
-          y[2] = 2;
-          y[3] = 4;
-          y[4] = 5;
-          return y;
-        }();
-
-        SmallVector<double> b = A * x_exact;
-
-        SECTION("builtin")
-        {
-          SECTION("CG no preconditioner")
-          {
-            LinearSolverOptions options;
-            options.library() = LSLibrary::builtin;
-            options.method()  = LSMethod::cg;
-            options.precond() = LSPrecond::none;
-
-            SmallVector<double> x{5};
-            x = zero;
-
-            LinearSolver solver{options};
-
-            solver.solveLocalSystem(A, x, b);
-            SmallVector error = x - x_exact;
-            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
-          }
-        }
 
-        SECTION("PETSc")
+        SECTION("Eigen3")
         {
-#ifdef PUGS_HAS_PETSC
+#ifdef PUGS_HAS_EIGEN3
 
-          SECTION("CG")
+          SECTION("BICGStab")
           {
             LinearSolverOptions options;
-            options.library() = LSLibrary::petsc;
-            options.method()  = LSMethod::cg;
+            options.library() = LSLibrary::eigen3;
+            options.method()  = LSMethod::bicgstab;
             options.precond() = LSPrecond::none;
             options.verbose() = true;
 
-            SECTION("CG no preconditioner")
+            SECTION("BICGStab no preconditioner")
             {
               options.precond() = LSPrecond::none;
 
-              SmallVector<double> x{5};
+              Vector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
+              Vector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
-            SECTION("CG Diagonal")
+            SECTION("BICGStab Diagonal")
             {
               options.precond() = LSPrecond::diagonal;
 
-              SmallVector<double> x{5};
+              Vector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
+              Vector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
-            SECTION("CG ICholesky")
+            SECTION("BICGStab ILU")
             {
-              options.precond() = LSPrecond::incomplete_cholesky;
+              options.precond() = LSPrecond::incomplete_LU;
 
-              SmallVector<double> x{5};
+              Vector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
+              Vector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
+          }
 
-            SECTION("CG AMG")
+          SECTION("GMRES")
+          {
+            LinearSolverOptions options;
+            options.library() = LSLibrary::eigen3;
+            options.method()  = LSMethod::gmres;
+            options.precond() = LSPrecond::none;
+
+            SECTION("GMRES no preconditioner")
             {
-              options.precond() = LSPrecond::amg;
+              options.precond() = LSPrecond::none;
 
-              SmallVector<double> x{5};
+              Vector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
+              Vector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
-          }
-
-          SECTION("Cholesky")
-          {
-            LinearSolverOptions options;
-            options.library() = LSLibrary::petsc;
-            options.method()  = LSMethod::cholesky;
-            options.precond() = LSPrecond::none;
 
-            SmallVector<double> x{5};
-            x = zero;
+            SECTION("GMRES Diagonal")
+            {
+              options.precond() = LSPrecond::diagonal;
+
+              Vector<double> x{5};
+              x = zero;
+
+              LinearSolver solver{options};
+
+              solver.solveLocalSystem(A, x, b);
+              Vector error = x - x_exact;
+              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+            }
+
+            SECTION("GMRES ILU")
+            {
+              options.precond() = LSPrecond::incomplete_LU;
+
+              Vector<double> x{5};
+              x = zero;
+
+              LinearSolver solver{options};
+
+              solver.solveLocalSystem(A, x, b);
+              Vector error = x - x_exact;
+              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+            }
+          }
+
+          SECTION("LU")
+          {
+            LinearSolverOptions options;
+            options.library() = LSLibrary::eigen3;
+            options.method()  = LSMethod::lu;
+            options.precond() = LSPrecond::none;
+
+            Vector<double> x{5};
+            x = zero;
 
             LinearSolver solver{options};
 
             solver.solveLocalSystem(A, x, b);
-            SmallVector error = x - x_exact;
+            Vector error = x - x_exact;
             REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
-#else    // PUGS_HAS_PETSC
-          SECTION("PETSc not linked")
+#else    // PUGS_HAS_EIGEN3
+          SECTION("Eigen3 not linked")
           {
             LinearSolverOptions options;
-            options.library() = LSLibrary::petsc;
+            options.library() = LSLibrary::eigen3;
             options.method()  = LSMethod::cg;
             options.precond() = LSPrecond::none;
 
-            REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+            REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
           }
-#endif   // PUGS_HAS_PETSC
+#endif   // PUGS_HAS_EIGEN3
         }
       }
+    }
+  }
 
-      SECTION("none symmetric system")
+  SECTION("Dense Matrices")
+  {
+    SECTION("check linear solvers")
+    {
+      SECTION("symmetric system")
       {
         SmallMatrix<double> A{5};
         A = zero;
@@ -706,20 +847,20 @@ TEST_CASE("LinearSolver", "[algebra]")
         A(0, 0) = 2;
         A(0, 1) = -1;
 
-        A(1, 0) = -0.2;
+        A(1, 0) = -1;
         A(1, 1) = 2;
         A(1, 2) = -1;
 
         A(2, 1) = -1;
-        A(2, 2) = 4;
-        A(2, 3) = -2;
+        A(2, 2) = 2;
+        A(2, 3) = -1;
 
         A(3, 2) = -1;
         A(3, 3) = 2;
-        A(3, 4) = -0.1;
+        A(3, 4) = -1;
 
-        A(4, 3) = 1;
-        A(4, 4) = 3;
+        A(4, 3) = -1;
+        A(4, 4) = 2;
 
         SmallVector<const double> x_exact = [] {
           SmallVector<double> y{5};
@@ -735,11 +876,11 @@ TEST_CASE("LinearSolver", "[algebra]")
 
         SECTION("builtin")
         {
-          SECTION("BICGStab no preconditioner")
+          SECTION("CG no preconditioner")
           {
             LinearSolverOptions options;
             options.library() = LSLibrary::builtin;
-            options.method()  = LSMethod::bicgstab;
+            options.method()  = LSMethod::cg;
             options.precond() = LSPrecond::none;
 
             SmallVector<double> x{5};
@@ -753,19 +894,19 @@ TEST_CASE("LinearSolver", "[algebra]")
           }
         }
 
-        SECTION("PETSc")
+        SECTION("Eigen3")
         {
-#ifdef PUGS_HAS_PETSC
+#ifdef PUGS_HAS_EIGEN3
 
-          SECTION("BICGStab")
+          SECTION("CG")
           {
             LinearSolverOptions options;
-            options.library() = LSLibrary::petsc;
-            options.method()  = LSMethod::bicgstab;
+            options.library() = LSLibrary::eigen3;
+            options.method()  = LSMethod::cg;
             options.precond() = LSPrecond::none;
             options.verbose() = true;
 
-            SECTION("BICGStab no preconditioner")
+            SECTION("CG no preconditioner")
             {
               options.precond() = LSPrecond::none;
 
@@ -779,7 +920,7 @@ TEST_CASE("LinearSolver", "[algebra]")
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
-            SECTION("BICGStab Diagonal")
+            SECTION("CG Diagonal")
             {
               options.precond() = LSPrecond::diagonal;
 
@@ -793,45 +934,75 @@ TEST_CASE("LinearSolver", "[algebra]")
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
-            SECTION("BICGStab ILU")
+            SECTION("CG ICholesky")
             {
-              options.precond() = LSPrecond::incomplete_LU;
+              options.precond() = LSPrecond::incomplete_cholesky;
 
               SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
-              solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
-              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+                                  "error: incomplete cholesky is not available for dense matrices in Eigen3");
+            }
+
+            SECTION("CG AMG")
+            {
+              options.precond() = LSPrecond::amg;
+
+              SmallVector<double> x{5};
+              x = zero;
+
+              REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not an Eigen3 preconditioner!");
             }
           }
 
-          SECTION("BICGStab2")
+          SECTION("Cholesky")
           {
             LinearSolverOptions options;
-            options.library() = LSLibrary::petsc;
-            options.method()  = LSMethod::bicgstab2;
+            options.library() = LSLibrary::eigen3;
+            options.method()  = LSMethod::cholesky;
             options.precond() = LSPrecond::none;
 
-            SECTION("BICGStab2 no preconditioner")
-            {
-              options.precond() = LSPrecond::none;
+            SmallVector<double> x{5};
+            x = zero;
 
-              SmallVector<double> x{5};
-              x = zero;
+            LinearSolver solver{options};
 
-              LinearSolver solver{options};
+            solver.solveLocalSystem(A, x, b);
+            SmallVector error = x - x_exact;
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+          }
 
-              solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
-              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
-            }
+#else    // PUGS_HAS_EIGEN3
+          SECTION("Eigen3 not linked")
+          {
+            LinearSolverOptions options;
+            options.library() = LSLibrary::eigen3;
+            options.method()  = LSMethod::cg;
+            options.precond() = LSPrecond::none;
 
-            SECTION("BICGStab2 Diagonal")
+            REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+          }
+#endif   // PUGS_HAS_EIGEN3
+        }
+
+        SECTION("PETSc")
+        {
+#ifdef PUGS_HAS_PETSC
+
+          SECTION("CG")
+          {
+            LinearSolverOptions options;
+            options.library() = LSLibrary::petsc;
+            options.method()  = LSMethod::cg;
+            options.precond() = LSPrecond::none;
+            options.verbose() = true;
+
+            SECTION("CG no preconditioner")
             {
-              options.precond() = LSPrecond::diagonal;
+              options.precond() = LSPrecond::none;
 
               SmallVector<double> x{5};
               x = zero;
@@ -842,18 +1013,10 @@ TEST_CASE("LinearSolver", "[algebra]")
               SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
-          }
 
-          SECTION("GMRES")
-          {
-            LinearSolverOptions options;
-            options.library() = LSLibrary::petsc;
-            options.method()  = LSMethod::gmres;
-            options.precond() = LSPrecond::none;
-
-            SECTION("GMRES no preconditioner")
+            SECTION("CG Diagonal")
             {
-              options.precond() = LSPrecond::none;
+              options.precond() = LSPrecond::diagonal;
 
               SmallVector<double> x{5};
               x = zero;
@@ -865,9 +1028,9 @@ TEST_CASE("LinearSolver", "[algebra]")
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
-            SECTION("GMRES Diagonal")
+            SECTION("CG ICholesky")
             {
-              options.precond() = LSPrecond::diagonal;
+              options.precond() = LSPrecond::incomplete_cholesky;
 
               SmallVector<double> x{5};
               x = zero;
@@ -879,9 +1042,9 @@ TEST_CASE("LinearSolver", "[algebra]")
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
-            SECTION("GMRES ILU")
+            SECTION("CG AMG")
             {
-              options.precond() = LSPrecond::incomplete_LU;
+              options.precond() = LSPrecond::amg;
 
               SmallVector<double> x{5};
               x = zero;
@@ -894,11 +1057,11 @@ TEST_CASE("LinearSolver", "[algebra]")
             }
           }
 
-          SECTION("LU")
+          SECTION("Cholesky")
           {
             LinearSolverOptions options;
             options.library() = LSLibrary::petsc;
-            options.method()  = LSMethod::lu;
+            options.method()  = LSMethod::cholesky;
             options.precond() = LSPrecond::none;
 
             SmallVector<double> x{5};
@@ -924,6 +1087,387 @@ TEST_CASE("LinearSolver", "[algebra]")
 #endif   // PUGS_HAS_PETSC
         }
       }
+
+      SECTION("none symmetric system")
+      {
+        SECTION("Dense matrix")
+        {
+          SmallMatrix<double> A{5};
+          A = zero;
+
+          A(0, 0) = 2;
+          A(0, 1) = -1;
+
+          A(1, 0) = -0.2;
+          A(1, 1) = 2;
+          A(1, 2) = -1;
+
+          A(2, 1) = -1;
+          A(2, 2) = 4;
+          A(2, 3) = -2;
+
+          A(3, 2) = -1;
+          A(3, 3) = 2;
+          A(3, 4) = -0.1;
+
+          A(4, 3) = 1;
+          A(4, 4) = 3;
+
+          SmallVector<const double> x_exact = [] {
+            SmallVector<double> y{5};
+            y[0] = 1;
+            y[1] = 3;
+            y[2] = 2;
+            y[3] = 4;
+            y[4] = 5;
+            return y;
+          }();
+
+          SmallVector<double> b = A * x_exact;
+
+          SECTION("builtin")
+          {
+            SECTION("BICGStab no preconditioner")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::builtin;
+              options.method()  = LSMethod::bicgstab;
+              options.precond() = LSPrecond::none;
+
+              SmallVector<double> x{5};
+              x = zero;
+
+              LinearSolver solver{options};
+
+              solver.solveLocalSystem(A, x, b);
+              SmallVector error = x - x_exact;
+              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+            }
+          }
+
+          SECTION("Eigen3")
+          {
+#ifdef PUGS_HAS_EIGEN3
+
+            SECTION("BICGStab")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::bicgstab;
+              options.precond() = LSPrecond::none;
+              options.verbose() = true;
+
+              SECTION("BICGStab no preconditioner")
+              {
+                options.precond() = LSPrecond::none;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+              SECTION("BICGStab Diagonal")
+              {
+                options.precond() = LSPrecond::diagonal;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+              SECTION("BICGStab ILU")
+              {
+                options.precond() = LSPrecond::incomplete_LU;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+                                    "error: incomplete LU is not available for dense matrices in Eigen3");
+              }
+            }
+
+            SECTION("BICGStab2")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::bicgstab2;
+              options.precond() = LSPrecond::none;
+
+              SECTION("BICGStab2 no preconditioner")
+              {
+                options.precond() = LSPrecond::none;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                REQUIRE_THROWS_WITH(LinearSolver{options}, "error: BICGStab2 is not an Eigen3 linear solver!");
+              }
+            }
+
+            SECTION("GMRES")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::gmres;
+              options.precond() = LSPrecond::none;
+
+              SECTION("GMRES no preconditioner")
+              {
+                options.precond() = LSPrecond::none;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+              SECTION("GMRES Diagonal")
+              {
+                options.precond() = LSPrecond::diagonal;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+              SECTION("GMRES ILU")
+              {
+                options.precond() = LSPrecond::incomplete_LU;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+                                    "error: incomplete LU is not available for dense matrices in Eigen3");
+              }
+            }
+
+            SECTION("LU")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::lu;
+              options.precond() = LSPrecond::none;
+
+              SmallVector<double> x{5};
+              x = zero;
+
+              LinearSolver solver{options};
+
+              solver.solveLocalSystem(A, x, b);
+              SmallVector error = x - x_exact;
+              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+            }
+
+#else    // PUGS_HAS_EIGEN3
+            SECTION("Eigen3 not linked")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::cg;
+              options.precond() = LSPrecond::none;
+
+              REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+            }
+#endif   // PUGS_HAS_EIGEN3
+          }
+
+          SECTION("PETSc")
+          {
+#ifdef PUGS_HAS_PETSC
+
+            SECTION("BICGStab")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::bicgstab;
+              options.precond() = LSPrecond::none;
+              options.verbose() = true;
+
+              SECTION("BICGStab no preconditioner")
+              {
+                options.precond() = LSPrecond::none;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+              SECTION("BICGStab Diagonal")
+              {
+                options.precond() = LSPrecond::diagonal;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+              SECTION("BICGStab ILU")
+              {
+                options.precond() = LSPrecond::incomplete_LU;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+            }
+
+            SECTION("BICGStab2")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::bicgstab2;
+              options.precond() = LSPrecond::none;
+
+              SECTION("BICGStab2 no preconditioner")
+              {
+                options.precond() = LSPrecond::none;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+              SECTION("BICGStab2 Diagonal")
+              {
+                options.precond() = LSPrecond::diagonal;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+            }
+
+            SECTION("GMRES")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::gmres;
+              options.precond() = LSPrecond::none;
+
+              SECTION("GMRES no preconditioner")
+              {
+                options.precond() = LSPrecond::none;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+              SECTION("GMRES Diagonal")
+              {
+                options.precond() = LSPrecond::diagonal;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+              SECTION("GMRES ILU")
+              {
+                options.precond() = LSPrecond::incomplete_LU;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+            }
+
+            SECTION("LU")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::lu;
+              options.precond() = LSPrecond::none;
+
+              SmallVector<double> x{5};
+              x = zero;
+
+              LinearSolver solver{options};
+
+              solver.solveLocalSystem(A, x, b);
+              SmallVector error = x - x_exact;
+              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+            }
+
+#else    // PUGS_HAS_PETSC
+            SECTION("PETSc not linked")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::cg;
+              options.precond() = LSPrecond::none;
+
+              REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+            }
+#endif   // PUGS_HAS_PETSC
+          }
+        }
+      }
     }
   }
 }
diff --git a/tests/test_LinearSolverOptions.cpp b/tests/test_LinearSolverOptions.cpp
index 90519ddf9797f0d6fa745cbc34e2f80443eed33d..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!");
   }
@@ -135,6 +136,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]")
 
     const std::string library_list = R"(
   - builtin
+  - Eigen3
   - PETSc
 )";
 
diff --git a/tests/test_PugsUtils.cpp b/tests/test_PugsUtils.cpp
index fb86486527e3f1d91ce5c2f1e7d2653a1dea3bde..5ba8817fcaa34ddf31f0955ed9fd118f841bd87d 100644
--- a/tests/test_PugsUtils.cpp
+++ b/tests/test_PugsUtils.cpp
@@ -51,6 +51,7 @@ TEST_CASE("PugsUtils", "[utils]")
       os << "MPI:      " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
       os << "PETSc:    " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n';
       os << "SLEPc:    " << rang::style::bold << BuildInfo::slepcLibrary() << rang::style::reset << '\n';
+      os << "Eigen3:   " << rang::style::bold << BuildInfo::eigen3Library() << rang::style::reset << '\n';
       os << "HDF5:     " << rang::style::bold << BuildInfo::hdf5Library() << rang::style::reset << '\n';
       os << "SLURM:    " << rang::style::bold << BuildInfo::slurmLibrary() << rang::style::reset << '\n';
       os << "-------------------------------------------------------";
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);