From c72934b2bf57e85a4edabf650c9b439eac26b7c8 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Thu, 23 Jan 2025 09:00:27 +0100
Subject: [PATCH] Begin Eigen3 interfacing [ci-skip]
---
CMakeLists.txt | 30 +++
src/algebra/Eigen3Utils.hpp | 116 ++++++++++
src/algebra/LinearSolver.cpp | 55 +++++
src/algebra/LinearSolverOptions.hpp | 4 +
src/utils/BuildInfo.cpp | 14 ++
src/utils/BuildInfo.hpp | 1 +
src/utils/PugsUtils.cpp | 1 +
src/utils/pugs_config.hpp.in | 1 +
tests/CMakeLists.txt | 1 +
tests/test_Eigen3Utils.cpp | 107 +++++++++
tests/test_LinearSolver.cpp | 344 +++++++++++++++++++++++++++-
tests/test_LinearSolverOptions.cpp | 1 +
tests/test_PugsUtils.cpp | 1 +
13 files changed, 675 insertions(+), 1 deletion(-)
create mode 100644 src/algebra/Eigen3Utils.hpp
create mode 100644 tests/test_Eigen3Utils.cpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2e0b29900..52579983f 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/src/algebra/Eigen3Utils.hpp b/src/algebra/Eigen3Utils.hpp
new file mode 100644
index 000000000..14434fe94
--- /dev/null
+++ b/src/algebra/Eigen3Utils.hpp
@@ -0,0 +1,116 @@
+#ifndef EIGEN3_UTILS_HPP
+#define EIGEN3_UTILS_HPP
+
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_EIGEN3
+
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/SmallMatrix.hpp>
+#include <algebra/TinyMatrix.hpp>
+
+#include <eigen3/Eigen/Eigen>
+
+class Eigen3DenseMatrixEmbedder
+{
+ public:
+ using MatrixType = Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>;
+
+ private:
+ MatrixType 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
+ numberOfRows() const
+ {
+ return m_eigen_matrix.rows();
+ }
+
+ PUGS_INLINE
+ size_t
+ numberOfColumns() const
+ {
+ return m_eigen_matrix.cols();
+ }
+
+ PUGS_INLINE
+ MatrixType&
+ matrix()
+ {
+ return m_eigen_matrix;
+ }
+
+ PUGS_INLINE
+ const MatrixType&
+ matrix() const
+ {
+ return m_eigen_matrix;
+ }
+
+ template <size_t N>
+ Eigen3DenseMatrixEmbedder(const TinyMatrix<N>& A) : Eigen3DenseMatrixEmbedder{N, N, &A(0, 0)}
+ {}
+
+ Eigen3DenseMatrixEmbedder(const SmallMatrix<double>& A)
+ : Eigen3DenseMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), &A(0, 0)}
+ {}
+
+ ~Eigen3DenseMatrixEmbedder() = default;
+};
+
+class Eigen3SparseMatrixEmbedder
+{
+ public:
+ using MatrixType = Eigen::Map<const Eigen::SparseMatrix<double, Eigen::RowMajor>>;
+
+ private:
+ MatrixType 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
+ MatrixType&
+ matrix()
+ {
+ return m_eigen_matrix;
+ }
+
+ PUGS_INLINE
+ const MatrixType&
+ 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;
+};
+
+#endif // PUGS_HAS_EIGEN3
+
+#endif // EIGEN3_UTILS_HPP
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index 99bff927c..7504ced5f 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -18,6 +18,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 +64,28 @@ struct LinearSolver::Internals
}
}
+ static void
+ checkEigen3Method(const LSMethod method)
+ {
+ switch (method) {
+ case LSMethod::cg:
+#warning clean-up
+ // case LSMethod::bicgstab:
+ // case LSMethod::bicgstab2:
+ // 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 +119,27 @@ struct LinearSolver::Internals
}
}
+ static void
+ checkEigen3Precond(const LSPrecond precond)
+ {
+ switch (precond) {
+ case LSPrecond::none:
+#warning clean-up
+ // case LSPrecond::amg:
+ // 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 +168,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());
diff --git a/src/algebra/LinearSolverOptions.hpp b/src/algebra/LinearSolverOptions.hpp
index 940ce1ed6..97a174878 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/utils/BuildInfo.cpp b/src/utils/BuildInfo.cpp
index 81d756775..2413bdbae 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 df4a32d93..86f7c04f0 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/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index 1ac38b326..ac64a50cf 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/pugs_config.hpp.in b/src/utils/pugs_config.hpp.in
index 7402a526d..a07d8c314 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 23a4ac309..1dd13e6e3 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 000000000..0c51aad3c
--- /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_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 11195faab..01c71cbea 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,72 @@ 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::bicgstab2;
+ 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));
+
+ options.precond() = LSPrecond::amg;
+ 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 +277,105 @@ 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)));
+ }
+#warning test all options
+ // 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;
+
+ // 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("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
@@ -753,6 +924,177 @@ TEST_CASE("LinearSolver", "[algebra]")
}
}
+ 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};
+
+ 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_EIGEN3
+ SECTION("Eigen3 not linked")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ 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
diff --git a/tests/test_LinearSolverOptions.cpp b/tests/test_LinearSolverOptions.cpp
index 90519ddf9..620959714 100644
--- a/tests/test_LinearSolverOptions.cpp
+++ b/tests/test_LinearSolverOptions.cpp
@@ -135,6 +135,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 fb8648652..5ba8817fc 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 << "-------------------------------------------------------";
--
GitLab