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 1/8] 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
From b1ebe1e7048fd1947e00507c72451f6ba52da0d3 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Sat, 25 Jan 2025 12:57:19 +0100
Subject: [PATCH 2/8] Plug a bunch of Eigen3 linear solvers
---
src/algebra/Eigen3Utils.hpp | 18 +-
src/algebra/LinearSolver.cpp | 197 ++++++++-
tests/test_LinearSolver.cpp | 754 +++++++++++++++++++++++------------
3 files changed, 686 insertions(+), 283 deletions(-)
diff --git a/src/algebra/Eigen3Utils.hpp b/src/algebra/Eigen3Utils.hpp
index 14434fe94..5224a4b0b 100644
--- a/src/algebra/Eigen3Utils.hpp
+++ b/src/algebra/Eigen3Utils.hpp
@@ -14,10 +14,11 @@
class Eigen3DenseMatrixEmbedder
{
public:
- using MatrixType = Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>;
+ using Eigen3MatrixType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+ using Eigen3MapMatrixType = Eigen::Map<const Eigen3MatrixType>;
private:
- MatrixType m_eigen_matrix;
+ 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)}
@@ -39,14 +40,14 @@ class Eigen3DenseMatrixEmbedder
}
PUGS_INLINE
- MatrixType&
+ Eigen3MapMatrixType&
matrix()
{
return m_eigen_matrix;
}
PUGS_INLINE
- const MatrixType&
+ const Eigen3MapMatrixType&
matrix() const
{
return m_eigen_matrix;
@@ -66,10 +67,11 @@ class Eigen3DenseMatrixEmbedder
class Eigen3SparseMatrixEmbedder
{
public:
- using MatrixType = Eigen::Map<const Eigen::SparseMatrix<double, Eigen::RowMajor>>;
+ using Eigen3MatrixType = Eigen::SparseMatrix<double, Eigen::RowMajor>;
+ using Eigen3MapMatrixType = Eigen::Map<const Eigen3MatrixType>;
private:
- MatrixType m_eigen_matrix;
+ Eigen3MapMatrixType m_eigen_matrix;
public:
PUGS_INLINE
@@ -87,14 +89,14 @@ class Eigen3SparseMatrixEmbedder
}
PUGS_INLINE
- MatrixType&
+ Eigen3MapMatrixType&
matrix()
{
return m_eigen_matrix;
}
PUGS_INLINE
- const MatrixType&
+ const Eigen3MapMatrixType&
matrix() const
{
return m_eigen_matrix;
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index 7504ced5f..1ebb59663 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
@@ -69,15 +74,12 @@ struct LinearSolver::Internals
{
switch (method) {
case LSMethod::cg:
-#warning clean-up
- // case LSMethod::bicgstab:
- // case LSMethod::bicgstab2:
- // case LSMethod::gmres:
- // case LSMethod::lu:
- // case LSMethod::cholesky:
- {
- break;
- }
+ 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!");
@@ -124,14 +126,11 @@ struct LinearSolver::Internals
{
switch (precond) {
case LSPrecond::none:
-#warning clean-up
- // case LSPrecond::amg:
- // case LSPrecond::diagonal:
- // case LSPrecond::incomplete_cholesky:
- // case LSPrecond::incomplete_LU:
- {
- break;
- }
+ 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!");
@@ -215,6 +214,163 @@ struct LinearSolver::Internals
}
}
+#ifdef PUGS_HAS_EIGEN3
+ template <typename PreconditionerType,
+ typename Eigen3MatrixEmbedderType,
+ typename SolutionVectorType,
+ typename RHSVectorType>
+ static void
+ _preconditionedEigen3SolveLocalSystem(const Eigen3MatrixEmbedderType& eigen3_A,
+ SolutionVectorType& x,
+ const RHSVectorType& b,
+ const LinearSolverOptions& options)
+ {
+ Eigen::Map<Eigen::VectorX<double>> eigen3_x{(x.size() > 0) ? (&x[0]) : nullptr, static_cast<int>(x.size())};
+ Eigen::Map<const Eigen::VectorX<double>> eigen3_b{(b.size() > 0) ? (&b[0]) : nullptr, 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, typename SolutionVectorType, typename RHSVectorType>
+ static void
+ _eigen3SolveLocalSystem(const Eigen3MatrixEmbedderType& eigen3_A,
+ SolutionVectorType& x,
+ const RHSVectorType& b,
+ const LinearSolverOptions& options)
+ {
+ switch (options.precond()) {
+ case LSPrecond::none: {
+ _preconditionedEigen3SolveLocalSystem<Eigen::IdentityPreconditioner, Eigen3MatrixEmbedderType, SolutionVectorType,
+ RHSVectorType>(eigen3_A, x, b, options);
+ break;
+ }
+ case LSPrecond::diagonal: {
+ _preconditionedEigen3SolveLocalSystem<Eigen::DiagonalPreconditioner<double>, Eigen3MatrixEmbedderType,
+ SolutionVectorType, RHSVectorType>(eigen3_A, x, b, options);
+ break;
+ }
+ case LSPrecond::incomplete_cholesky: {
+ if constexpr (std::is_same_v<Eigen3SparseMatrixEmbedder, Eigen3MatrixEmbedderType>) {
+ _preconditionedEigen3SolveLocalSystem<Eigen::IncompleteCholesky<double>, Eigen3MatrixEmbedderType,
+ SolutionVectorType, RHSVectorType>(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,
+ SolutionVectorType, RHSVectorType>(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
+ }
+ }
+
+ template <typename DataType, typename SolutionVectorType, typename RHSVectorType>
+ static void
+ eigen3SolveLocalSystem(const CRSMatrix<DataType>& A,
+ SolutionVectorType& x,
+ const RHSVectorType& b,
+ const LinearSolverOptions& options)
+ {
+ Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
+
+ Eigen3SparseMatrixEmbedder eigen3_A{A};
+ _eigen3SolveLocalSystem(eigen3_A, x, b, options);
+ }
+
+ template <typename DataType, typename SolutionVectorType, typename RHSVectorType>
+ static void
+ eigen3SolveLocalSystem(const SmallMatrix<DataType>& A,
+ SolutionVectorType& x,
+ const RHSVectorType& b,
+ const LinearSolverOptions& options)
+ {
+ Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
+
+ Eigen3DenseMatrixEmbedder eigen3_A{A};
+ _eigen3SolveLocalSystem(eigen3_A, x, b, options);
+ }
+
+#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*)
@@ -357,6 +513,13 @@ struct LinearSolver::Internals
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, options);
+ break;
+ }
case LSLibrary::petsc: {
// not covered since if PETSc is not linked this point is
// unreachable: LinearSolver throws an exception at construction
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 01c71cbea..4d01df2f0 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -167,9 +167,6 @@ TEST_CASE("LinearSolver", "[algebra]")
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));
@@ -196,9 +193,6 @@ TEST_CASE("LinearSolver", "[algebra]")
options.precond() = LSPrecond::incomplete_cholesky;
REQUIRE_NOTHROW(linear_solver.checkOptions(options));
-
- options.precond() = LSPrecond::amg;
- REQUIRE_NOTHROW(linear_solver.checkOptions(options));
}
}
@@ -302,48 +296,44 @@ TEST_CASE("LinearSolver", "[algebra]")
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;
+ SECTION("CG Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
- // LinearSolver solver{options};
+ Vector<double> x{5};
+ x = zero;
- // 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)));
- // }
+ LinearSolver solver{options};
- // SECTION("CG ICholesky")
- // {
- // options.precond() = LSPrecond::incomplete_cholesky;
+ 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)));
+ }
- // Vector<double> x{5};
- // x = zero;
+ SECTION("CG ICholesky")
+ {
+ options.precond() = LSPrecond::incomplete_cholesky;
- // LinearSolver solver{options};
+ Vector<double> x{5};
+ x = zero;
- // 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)));
- // }
+ LinearSolver solver{options};
- // SECTION("CG AMG")
- // {
- // options.precond() = LSPrecond::amg;
+ 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)));
+ }
- // Vector<double> x{5};
- // x = zero;
+ SECTION("CG AMG")
+ {
+ options.precond() = LSPrecond::amg;
- // LinearSolver solver{options};
+ Vector<double> x{5};
+ x = zero;
- // 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)));
- // }
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not an Eigen3 preconditioner!");
+ }
}
SECTION("Cholesky")
@@ -706,170 +696,186 @@ TEST_CASE("LinearSolver", "[algebra]")
}
#endif // PUGS_HAS_PETSC
}
- }
- }
- }
- SECTION("Dense Matrices")
- {
- SECTION("check linear solvers")
- {
- SECTION("symmetric system")
- {
- SmallMatrix<double> A{5};
- A = zero;
+ SECTION("Eigen3")
+ {
+#ifdef PUGS_HAS_EIGEN3
- A(0, 0) = 2;
- A(0, 1) = -1;
+ SECTION("BICGStab")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::bicgstab;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
- A(1, 0) = -1;
- A(1, 1) = 2;
- A(1, 2) = -1;
+ SECTION("BICGStab no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- A(2, 1) = -1;
- A(2, 2) = 2;
- A(2, 3) = -1;
+ Vector<double> x{5};
+ x = zero;
- A(3, 2) = -1;
- A(3, 3) = 2;
- A(3, 4) = -1;
+ LinearSolver solver{options};
- A(4, 3) = -1;
- A(4, 4) = 2;
+ 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)));
+ }
- 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;
- }();
+ SECTION("BICGStab Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
- SmallVector<double> b = A * x_exact;
+ Vector<double> x{5};
+ x = zero;
- SECTION("builtin")
- {
- SECTION("CG no preconditioner")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::builtin;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
+ LinearSolver solver{options};
- SmallVector<double> x{5};
- x = zero;
+ 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)));
+ }
- LinearSolver solver{options};
+ SECTION("BICGStab ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
- 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)));
- }
- }
+ Vector<double> x{5};
+ x = zero;
- SECTION("PETSc")
- {
-#ifdef PUGS_HAS_PETSC
+ LinearSolver solver{options};
- SECTION("CG")
+ 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("BICGStab2")
{
LinearSolverOptions options;
options.library() = LSLibrary::petsc;
- options.method() = LSMethod::cg;
+ options.method() = LSMethod::bicgstab2;
options.precond() = LSPrecond::none;
- options.verbose() = true;
- SECTION("CG no preconditioner")
+ SECTION("BICGStab2 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("BICGStab2 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("GMRES")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::gmres;
+ options.precond() = LSPrecond::none;
+
+ SECTION("GMRES no preconditioner")
{
- options.precond() = LSPrecond::incomplete_cholesky;
+ 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 AMG")
+ SECTION("GMRES Diagonal")
{
- options.precond() = LSPrecond::amg;
+ 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("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("Cholesky")
+ SECTION("LU")
{
LinearSolverOptions options;
options.library() = LSLibrary::petsc;
- options.method() = LSMethod::cholesky;
+ options.method() = LSMethod::lu;
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)));
}
-#else // PUGS_HAS_PETSC
- SECTION("PETSc not linked")
+#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: 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;
@@ -877,20 +883,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};
@@ -906,11 +912,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};
@@ -928,15 +934,15 @@ TEST_CASE("LinearSolver", "[algebra]")
{
#ifdef PUGS_HAS_EIGEN3
- SECTION("BICGStab")
+ SECTION("CG")
{
LinearSolverOptions options;
options.library() = LSLibrary::eigen3;
- options.method() = LSMethod::bicgstab;
+ options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;
options.verbose() = true;
- SECTION("BICGStab no preconditioner")
+ SECTION("CG no preconditioner")
{
options.precond() = LSPrecond::none;
@@ -950,7 +956,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;
@@ -964,45 +970,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;
@@ -1013,18 +1049,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;
@@ -1036,9 +1064,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;
@@ -1050,9 +1078,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;
@@ -1065,11 +1093,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};
@@ -1082,33 +1110,64 @@ TEST_CASE("LinearSolver", "[algebra]")
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
-#else // PUGS_HAS_EIGEN3
- SECTION("Eigen3 not linked")
+#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: Eigen3 is not linked to pugs. Cannot use it!");
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
}
-#endif // PUGS_HAS_EIGEN3
+#endif // PUGS_HAS_PETSC
}
+ }
- SECTION("PETSc")
+ SECTION("none symmetric system")
+ {
+ SECTION("Dense matrix")
{
-#ifdef PUGS_HAS_PETSC
+ SmallMatrix<double> A{5};
+ A = zero;
- SECTION("BICGStab")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::bicgstab;
- options.precond() = LSPrecond::none;
- options.verbose() = true;
+ 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};
@@ -1120,60 +1179,135 @@ 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("BICGStab Diagonal")
+ SECTION("Eigen3")
+ {
+#ifdef PUGS_HAS_EIGEN3
+
+ SECTION("BICGStab")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::bicgstab;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("BICGStab no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ SmallVector<double> x{5};
+ x = zero;
- 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)));
+ 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("BICGStab ILU")
+ SECTION("BICGStab2")
{
- options.precond() = LSPrecond::incomplete_LU;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::bicgstab2;
+ options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("BICGStab2 no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ SmallVector<double> x{5};
+ x = zero;
- 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(LinearSolver{options}, "error: BICGStab2 is not an Eigen3 linear solver!");
+ }
}
- }
-
- SECTION("BICGStab2")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::bicgstab2;
- options.precond() = LSPrecond::none;
- SECTION("BICGStab2 no preconditioner")
+ SECTION("GMRES")
{
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::gmres;
options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("GMRES no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ SmallVector<double> x{5};
+ x = zero;
- 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)));
+ 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("BICGStab2 Diagonal")
+ SECTION("LU")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::lu;
+ options.precond() = LSPrecond::none;
SmallVector<double> x{5};
x = zero;
@@ -1184,46 +1318,167 @@ TEST_CASE("LinearSolver", "[algebra]")
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("GMRES")
+ SECTION("PETSc")
{
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::gmres;
- options.precond() = LSPrecond::none;
+#ifdef PUGS_HAS_PETSC
- SECTION("GMRES no preconditioner")
+ SECTION("BICGStab")
{
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab;
options.precond() = LSPrecond::none;
+ options.verbose() = true;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("BICGStab no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ SmallVector<double> x{5};
+ x = zero;
- 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)));
+ 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("GMRES Diagonal")
+ SECTION("BICGStab2")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab2;
+ options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("BICGStab2 no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ SmallVector<double> x{5};
+ x = zero;
- 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)));
+ 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 ILU")
+ SECTION("GMRES")
{
- options.precond() = LSPrecond::incomplete_LU;
+ 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;
@@ -1234,36 +1489,19 @@ 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("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;
+ 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!");
- }
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
#endif // PUGS_HAS_PETSC
+ }
}
}
}
--
GitLab
From a9e3c1ffd2b64b58655f2d7220f08aaac362e81f Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Wed, 29 Jan 2025 02:05:16 +0100
Subject: [PATCH 3/8] Refactor linear solver (and improve slightly
performances)
---
src/algebra/DenseMatrixWrapper.hpp | 59 ++++++++++
src/algebra/Eigen3Utils.hpp | 43 ++++++--
src/algebra/LinearSolver.cpp | 170 ++++++++++++-----------------
src/algebra/LinearSolver.hpp | 131 ++++++++++++++++++++--
src/algebra/PETScUtils.hpp | 9 +-
src/algebra/VectorWrapper.hpp | 47 ++++++++
src/utils/PugsTraits.hpp | 22 ++++
tests/test_LinearSolver.cpp | 42 +------
8 files changed, 362 insertions(+), 161 deletions(-)
create mode 100644 src/algebra/DenseMatrixWrapper.hpp
create mode 100644 src/algebra/VectorWrapper.hpp
diff --git a/src/algebra/DenseMatrixWrapper.hpp b/src/algebra/DenseMatrixWrapper.hpp
new file mode 100644
index 000000000..dee57f27e
--- /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
index 5224a4b0b..d9d61e027 100644
--- a/src/algebra/Eigen3Utils.hpp
+++ b/src/algebra/Eigen3Utils.hpp
@@ -6,8 +6,7 @@
#ifdef PUGS_HAS_EIGEN3
#include <algebra/CRSMatrix.hpp>
-#include <algebra/SmallMatrix.hpp>
-#include <algebra/TinyMatrix.hpp>
+#include <algebra/DenseMatrixWrapper.hpp>
#include <eigen3/Eigen/Eigen>
@@ -25,6 +24,13 @@ class Eigen3DenseMatrixEmbedder
{}
public:
+ PUGS_INLINE
+ size_t
+ isSquare() const
+ {
+ return (m_eigen_matrix.rows() == m_eigen_matrix.cols());
+ }
+
PUGS_INLINE
size_t
numberOfRows() const
@@ -53,12 +59,8 @@ class Eigen3DenseMatrixEmbedder
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(const DenseMatrixWrapper<double>& A)
+ : Eigen3DenseMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), A.ptr()}
{}
~Eigen3DenseMatrixEmbedder() = default;
@@ -88,6 +90,13 @@ class Eigen3SparseMatrixEmbedder
return m_eigen_matrix.cols();
}
+ PUGS_INLINE
+ size_t
+ isSquare() const
+ {
+ return (m_eigen_matrix.rows() == m_eigen_matrix.cols());
+ }
+
PUGS_INLINE
Eigen3MapMatrixType&
matrix()
@@ -113,6 +122,24 @@ class Eigen3SparseMatrixEmbedder
~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/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index 1ebb59663..52bf8fc39 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -187,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
@@ -215,18 +215,15 @@ struct LinearSolver::Internals
}
#ifdef PUGS_HAS_EIGEN3
- template <typename PreconditionerType,
- typename Eigen3MatrixEmbedderType,
- typename SolutionVectorType,
- typename RHSVectorType>
+ template <typename PreconditionerType, typename Eigen3MatrixEmbedderType>
static void
_preconditionedEigen3SolveLocalSystem(const Eigen3MatrixEmbedderType& eigen3_A,
- SolutionVectorType& x,
- const RHSVectorType& b,
+ VectorWrapper<double>& x,
+ const VectorWrapper<const double>& b,
const LinearSolverOptions& options)
{
- Eigen::Map<Eigen::VectorX<double>> eigen3_x{(x.size() > 0) ? (&x[0]) : nullptr, static_cast<int>(x.size())};
- Eigen::Map<const Eigen::VectorX<double>> eigen3_b{(b.size() > 0) ? (&b[0]) : nullptr, static_cast<int>(b.size())};
+ 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;
@@ -287,28 +284,29 @@ struct LinearSolver::Internals
}
}
- template <typename Eigen3MatrixEmbedderType, typename SolutionVectorType, typename RHSVectorType>
+ template <typename Eigen3MatrixEmbedderType>
static void
_eigen3SolveLocalSystem(const Eigen3MatrixEmbedderType& eigen3_A,
- SolutionVectorType& x,
- const RHSVectorType& b,
+ VectorWrapper<double>& x,
+ const VectorWrapper<const double>& b,
const LinearSolverOptions& options)
{
switch (options.precond()) {
case LSPrecond::none: {
- _preconditionedEigen3SolveLocalSystem<Eigen::IdentityPreconditioner, Eigen3MatrixEmbedderType, SolutionVectorType,
- RHSVectorType>(eigen3_A, x, b, options);
+ _preconditionedEigen3SolveLocalSystem<Eigen::IdentityPreconditioner, Eigen3MatrixEmbedderType>(eigen3_A, x, b,
+ options);
break;
}
case LSPrecond::diagonal: {
- _preconditionedEigen3SolveLocalSystem<Eigen::DiagonalPreconditioner<double>, Eigen3MatrixEmbedderType,
- SolutionVectorType, RHSVectorType>(eigen3_A, x, b, options);
+ _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,
- SolutionVectorType, RHSVectorType>(eigen3_A, x, b, options);
+ _preconditionedEigen3SolveLocalSystem<Eigen::IncompleteCholesky<double>, Eigen3MatrixEmbedderType>(eigen3_A, x,
+ b, options);
} else {
throw NormalError("incomplete cholesky is not available for dense matrices in Eigen3");
}
@@ -316,8 +314,8 @@ struct LinearSolver::Internals
}
case LSPrecond::incomplete_LU: {
if constexpr (std::is_same_v<Eigen3SparseMatrixEmbedder, Eigen3MatrixEmbedderType>) {
- _preconditionedEigen3SolveLocalSystem<Eigen::IncompleteLUT<double>, Eigen3MatrixEmbedderType,
- SolutionVectorType, RHSVectorType>(eigen3_A, x, b, options);
+ _preconditionedEigen3SolveLocalSystem<Eigen::IncompleteLUT<double>, Eigen3MatrixEmbedderType>(eigen3_A, x, b,
+ options);
} else {
throw NormalError("incomplete LU is not available for dense matrices in Eigen3");
}
@@ -330,39 +328,12 @@ struct LinearSolver::Internals
// LCOV_EXCL_STOP
}
}
-
- template <typename DataType, typename SolutionVectorType, typename RHSVectorType>
- static void
- eigen3SolveLocalSystem(const CRSMatrix<DataType>& A,
- SolutionVectorType& x,
- const RHSVectorType& b,
- const LinearSolverOptions& options)
- {
- Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
-
- Eigen3SparseMatrixEmbedder eigen3_A{A};
- _eigen3SolveLocalSystem(eigen3_A, x, b, options);
- }
-
- template <typename DataType, typename SolutionVectorType, typename RHSVectorType>
- static void
- eigen3SolveLocalSystem(const SmallMatrix<DataType>& A,
- SolutionVectorType& x,
- const RHSVectorType& b,
- const LinearSolverOptions& options)
- {
- Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
-
- Eigen3DenseMatrixEmbedder eigen3_A{A};
- _eigen3SolveLocalSystem(eigen3_A, x, b, options);
- }
-
#else // PUGS_HAS_EIGEN3
// LCOV_EXCL_START
template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
static void
- eigen3SolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&)
+ _eigen3SolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&)
{
checkHasLibrary(LSLibrary::eigen3);
throw UnexpectedError("unexpected situation should not reach this point!");
@@ -379,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);
@@ -491,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!");
@@ -499,40 +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::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, options);
- 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, options);
- break;
- }
- default: {
- throw UnexpectedError(::name(options.library()) + " cannot solve local systems");
- }
- // LCOV_EXCL_STOP
- }
- }
};
bool
@@ -548,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 e90402c8d..2d03a9cf5 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/PETScUtils.hpp b/src/algebra/PETScUtils.hpp
index c4db9ee59..c6d2e8be8 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 000000000..0bdeefbfd
--- /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/utils/PugsTraits.hpp b/src/utils/PugsTraits.hpp
index 2e1d40957..3c7a4e40a 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/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 4d01df2f0..76ef6aa2b 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -752,46 +752,10 @@ TEST_CASE("LinearSolver", "[algebra]")
}
}
- SECTION("BICGStab2")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::bicgstab2;
- options.precond() = LSPrecond::none;
-
- SECTION("BICGStab2 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("BICGStab2 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")
{
LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
+ options.library() = LSLibrary::eigen3;
options.method() = LSMethod::gmres;
options.precond() = LSPrecond::none;
@@ -841,7 +805,7 @@ TEST_CASE("LinearSolver", "[algebra]")
SECTION("LU")
{
LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
+ options.library() = LSLibrary::eigen3;
options.method() = LSMethod::lu;
options.precond() = LSPrecond::none;
@@ -859,7 +823,7 @@ TEST_CASE("LinearSolver", "[algebra]")
SECTION("Eigen3 not linked")
{
LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
+ options.library() = LSLibrary::eigen3;
options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;
--
GitLab
From 1dcf50f0e6d6c8a9dc8a5de1590f5237b05c866b Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Thu, 30 Jan 2025 08:07:22 +0100
Subject: [PATCH 4/8] Prepare EigenvalueSolver to incoming use of Eigen3
---
src/algebra/CMakeLists.txt | 1 +
src/algebra/EigenvalueSolver.cpp | 250 +++++++++----
src/algebra/EigenvalueSolver.hpp | 109 ++++--
src/algebra/EigenvalueSolverOptions.cpp | 10 +
src/algebra/EigenvalueSolverOptions.hpp | 108 ++++++
src/language/modules/LinearSolverModule.cpp | 40 +-
src/utils/checkpointing/Checkpoint.cpp | 11 +-
.../EigenvalueSolverOptionsHFType.hpp | 15 +
.../LinearSolverOptionsHFType.hpp | 2 +-
src/utils/checkpointing/Resume.cpp | 10 +
tests/test_EigenvalueSolver.cpp | 344 ++++++++++++++----
tests/test_LinearSolverOptions.cpp | 1 +
tests/test_checkpointing_HFTypes.cpp | 3 +
13 files changed, 737 insertions(+), 167 deletions(-)
create mode 100644 src/algebra/EigenvalueSolverOptions.cpp
create mode 100644 src/algebra/EigenvalueSolverOptions.hpp
create mode 100644 src/utils/checkpointing/EigenvalueSolverOptionsHFType.hpp
diff --git a/src/algebra/CMakeLists.txt b/src/algebra/CMakeLists.txt
index d215b09b5..78c041c5c 100644
--- a/src/algebra/CMakeLists.txt
+++ b/src/algebra/CMakeLists.txt
@@ -3,6 +3,7 @@
add_library(
PugsAlgebra
EigenvalueSolver.cpp
+ EigenvalueSolverOptions.cpp
LinearSolver.cpp
LinearSolverOptions.cpp
PETScUtils.cpp)
diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index c71be58dd..5888d1c0b 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -2,10 +2,13 @@
#include <utils/pugs_config.hpp>
#ifdef PUGS_HAS_SLEPC
+#include <algebra/PETScUtils.hpp>
#include <slepc.h>
+#endif // PUGS_HAS_SLEPC
struct EigenvalueSolver::Internals
{
+#ifdef PUGS_HAS_SLEPC
static PetscReal
computeSmallestRealEigenvalueOfSymmetricMatrix(EPS& eps)
{
@@ -65,93 +68,210 @@ struct EigenvalueSolver::Internals
computeAllEigenvaluesOfSymmetricMatrixInInterval(eps, left_bound - 0.01 * std::abs(left_bound),
right_bound + 0.01 * std::abs(right_bound));
}
-};
-void
-EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues)
-{
- EPS eps;
+ static void
+ slepscComputeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues)
+ {
+ EPS eps;
- EPSCreate(PETSC_COMM_SELF, &eps);
- EPSSetOperators(eps, A, nullptr);
- EPSSetProblemType(eps, EPS_HEP);
+ EPSCreate(PETSC_COMM_SELF, &eps);
+ EPSSetOperators(eps, A, nullptr);
+ EPSSetProblemType(eps, EPS_HEP);
- Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+ Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
- PetscInt nb_eigenvalues;
- EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+ PetscInt nb_eigenvalues;
+ EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
- eigenvalues = SmallArray<double>(nb_eigenvalues);
- for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
- EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, nullptr, nullptr);
+ eigenvalues = SmallArray<double>(nb_eigenvalues);
+ for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
+ EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, nullptr, nullptr);
+ }
+
+ EPSDestroy(&eps);
}
- EPSDestroy(&eps);
+ static void
+ slepscComputeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvector_list)
+ {
+ EPS eps;
+
+ EPSCreate(PETSC_COMM_SELF, &eps);
+ EPSSetOperators(eps, A, nullptr);
+ EPSSetProblemType(eps, EPS_HEP);
+
+ Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+
+ PetscInt nb_eigenvalues;
+ EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+
+ eigenvalues = SmallArray<double>(nb_eigenvalues);
+ eigenvector_list.reserve(nb_eigenvalues);
+ for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
+ Vec Vr;
+ SmallVector<double> eigenvector{A.numberOfRows()};
+ VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr);
+ EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr);
+ VecDestroy(&Vr);
+ eigenvector_list.push_back(eigenvector);
+ }
+
+ EPSDestroy(&eps);
+ }
+
+ static void
+ slepscComputeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P)
+ {
+ EPS eps;
+
+ EPSCreate(PETSC_COMM_SELF, &eps);
+ EPSSetOperators(eps, A, nullptr);
+ EPSSetProblemType(eps, EPS_HEP);
+
+ Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+
+ PetscInt nb_eigenvalues;
+ EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+
+ eigenvalues = SmallArray<double>(nb_eigenvalues);
+ P = SmallMatrix<double>(nb_eigenvalues, nb_eigenvalues);
+
+ Array<double> eigenvector(nb_eigenvalues);
+ for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
+ Vec Vr;
+ VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr);
+ EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr);
+ VecDestroy(&Vr);
+ for (size_t j = 0; j < eigenvector.size(); ++j) {
+ P(j, i) = eigenvector[j];
+ }
+ }
+
+ EPSDestroy(&eps);
+ }
+#endif // PUGS_HAS_SLEPC
+};
+
+#ifdef PUGS_HAS_SLEPC
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues)
+{
+ Internals::slepscComputeForSymmetricMatrix(A, eigenvalues);
}
void
-EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
- SmallArray<double>& eigenvalues,
- std::vector<SmallVector<double>>& eigenvector_list)
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues)
{
- EPS eps;
-
- EPSCreate(PETSC_COMM_SELF, &eps);
- EPSSetOperators(eps, A, nullptr);
- EPSSetProblemType(eps, EPS_HEP);
-
- Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
-
- PetscInt nb_eigenvalues;
- EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
-
- eigenvalues = SmallArray<double>(nb_eigenvalues);
- eigenvector_list.reserve(nb_eigenvalues);
- for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
- Vec Vr;
- SmallVector<double> eigenvector{A.numberOfRows()};
- VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr);
- EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr);
- VecDestroy(&Vr);
- eigenvector_list.push_back(eigenvector);
- }
+ Internals::slepscComputeForSymmetricMatrix(A, eigenvalues);
+}
- EPSDestroy(&eps);
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors)
+{
+ Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
}
void
-EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
- SmallArray<double>& eigenvalues,
- SmallMatrix<double>& P)
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors)
{
- EPS eps;
+ Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
+}
- EPSCreate(PETSC_COMM_SELF, &eps);
- EPSSetOperators(eps, A, nullptr);
- EPSSetProblemType(eps, EPS_HEP);
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P)
+{
+ Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, P);
+}
- Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P)
+{
+ Internals::slepscComputeForSymmetricMatrix(A, eigenvalues, P);
+}
- PetscInt nb_eigenvalues;
- EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+#else // PUGS_HAS_SLEPC
- eigenvalues = SmallArray<double>(nb_eigenvalues);
- P = SmallMatrix<double>(nb_eigenvalues, nb_eigenvalues);
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&)
+{}
- Array<double> eigenvector(nb_eigenvalues);
- for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
- Vec Vr;
- VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr);
- EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr);
- VecDestroy(&Vr);
- for (size_t j = 0; j < eigenvector.size(); ++j) {
- P(j, i) = eigenvector[j];
- }
- }
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&, SmallArray<double>&)
+{}
- EPSDestroy(&eps);
-}
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>&,
+ SmallArray<double>&,
+ std::vector<SmallVector<double>>&)
+{}
+
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&,
+ SmallArray<double>&,
+ std::vector<SmallVector<double>>&)
+{}
+
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&, SmallMatrix<double>&)
+{}
+
+void
+EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&,
+ SmallArray<double>&,
+ SmallMatrix<double>&)
+{}
#endif // PUGS_HAS_SLEPC
-EigenvalueSolver::EigenvalueSolver() {}
+bool
+EigenvalueSolver::hasLibrary(const ESLibrary library)
+{
+ switch (library) {
+ case ESLibrary::eigen3: {
+#ifdef PUGS_HAS_EIGEN3
+ return true;
+#else
+ return false;
+#endif
+ }
+ case ESLibrary::slepsc: {
+#ifdef PUGS_HAS_SLEPC
+ return true;
+#else
+ return false;
+#endif
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("Eigenvalue library (" + ::name(library) + ") was not set!");
+ }
+ // LCOV_EXCL_STOP
+ }
+}
+
+void
+EigenvalueSolver::checkHasLibrary(const ESLibrary library)
+{
+ if (not hasLibrary(library)) {
+ // LCOV_EXCL_START
+ throw NormalError(::name(library) + " is not linked to pugs. Cannot use it!");
+ // LCOV_EXCL_STOP
+ }
+}
+
+EigenvalueSolver::EigenvalueSolver(const EigenvalueSolverOptions& options) : m_options{options}
+{
+ checkHasLibrary(options.library());
+}
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index a9b462509..eebd536e4 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -1,10 +1,13 @@
#ifndef EIGENVALUE_SOLVER_HPP
#define EIGENVALUE_SOLVER_HPP
-#include <algebra/PETScUtils.hpp>
+#include <algebra/EigenvalueSolverOptions.hpp>
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/DenseMatrixWrapper.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/SmallVector.hpp>
+#include <algebra/TinyMatrix.hpp>
#include <utils/Exceptions.hpp>
#include <utils/SmallArray.hpp>
@@ -13,57 +16,97 @@ class EigenvalueSolver
private:
struct Internals;
-#ifdef PUGS_HAS_SLEPC
- void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues);
+ const EigenvalueSolverOptions m_options;
- void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
- SmallArray<double>& eigenvalues,
- std::vector<SmallVector<double>>& eigenvectors);
+ void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues);
- void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
- SmallArray<double>& eigenvalues,
- SmallMatrix<double>& P);
-#endif // PUGS_HAS_SLEPC
+ void _slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues);
+
+ void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors);
+
+ void _slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors);
+
+ void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P);
+
+ void _slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P);
public:
+ static bool hasLibrary(ESLibrary library);
+ static void checkHasLibrary(const ESLibrary library);
+
template <typename MatrixType>
void
- computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] SmallArray<double>& eigenvalues)
+ computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues)
{
-#ifdef PUGS_HAS_SLEPC
- this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues);
-#else // PUGS_HAS_SLEPC
- throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
-#endif // PUGS_HAS_SLEPC
+ Assert((eigenvalues.size() - A.numberOfRows() == 0) and A.isSquare());
+
+ switch (m_options.library()) {
+ case ESLibrary::eigen3: {
+ throw NotImplementedError("Eigen3");
+ }
+ case ESLibrary::slepsc: {
+ this->_slepscComputeForSymmetricMatrix(A, eigenvalues);
+ break;
+ }
+ default: {
+ throw UnexpectedError(::name(m_options.library()) + " cannot compute for symmetric matrices");
+ }
+ }
}
template <typename MatrixType>
void
- computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A,
- [[maybe_unused]] SmallArray<double>& eigenvalues,
- [[maybe_unused]] std::vector<SmallVector<double>>& eigenvectors)
+ computeForSymmetricMatrix(const MatrixType& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors)
{
-#ifdef PUGS_HAS_SLEPC
- this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, eigenvectors);
-#else // PUGS_HAS_SLEPC
- throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
-#endif // PUGS_HAS_SLEPC
+ Assert((eigenvalues.size() - A.numberOfRows() == 0) and (A.numberOfRows() - eigenvectors.size() == 0) and
+ A.isSquare());
+
+ switch (m_options.library()) {
+ case ESLibrary::eigen3: {
+ throw NotImplementedError("Eigen3");
+ }
+ case ESLibrary::slepsc: {
+ this->_slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
+ break;
+ }
+ default: {
+ throw UnexpectedError(::name(m_options.library()) + " cannot compute for symmetric matrices");
+ }
+ }
}
template <typename MatrixType>
void
- computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A,
- [[maybe_unused]] SmallArray<double>& eigenvalues,
- [[maybe_unused]] SmallMatrix<double>& P)
+ computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues, SmallMatrix<double>& P)
{
-#ifdef PUGS_HAS_SLEPC
- this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, P);
-#else // PUGS_HAS_SLEPC
- throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
-#endif // PUGS_HAS_SLEPC
+ Assert((eigenvalues.size() - A.numberOfRows() == 0) and (A.numberOfRows() - P.numberOfRows() == 0) and
+ A.isSquare() and P.isSquare());
+
+ switch (m_options.library()) {
+ case ESLibrary::eigen3: {
+ throw NotImplementedError("Eigen3");
+ }
+ case ESLibrary::slepsc: {
+ this->_slepscComputeForSymmetricMatrix(A, eigenvalues, P);
+ break;
+ }
+ default: {
+ throw UnexpectedError(::name(m_options.library()) + " cannot compute for symmetric matrices");
+ }
+ }
}
- EigenvalueSolver();
+ EigenvalueSolver(const EigenvalueSolverOptions& options = EigenvalueSolverOptions::default_options);
~EigenvalueSolver() = default;
};
diff --git a/src/algebra/EigenvalueSolverOptions.cpp b/src/algebra/EigenvalueSolverOptions.cpp
new file mode 100644
index 000000000..2a4ef9549
--- /dev/null
+++ b/src/algebra/EigenvalueSolverOptions.cpp
@@ -0,0 +1,10 @@
+#include <algebra/EigenvalueSolverOptions.hpp>
+
+#include <rang.hpp>
+
+std::ostream&
+operator<<(std::ostream& os, const EigenvalueSolverOptions& options)
+{
+ os << " library: " << rang::style::bold << name(options.library()) << rang::style::reset << '\n';
+ return os;
+}
diff --git a/src/algebra/EigenvalueSolverOptions.hpp b/src/algebra/EigenvalueSolverOptions.hpp
new file mode 100644
index 000000000..d9be86be0
--- /dev/null
+++ b/src/algebra/EigenvalueSolverOptions.hpp
@@ -0,0 +1,108 @@
+#ifndef EIGENVALUE_SOLVER_OPTIONS_HPP
+#define EIGENVALUE_SOLVER_OPTIONS_HPP
+
+#include <utils/Exceptions.hpp>
+
+#include <iostream>
+
+enum class ESLibrary : int8_t
+{
+ ES__begin = 0,
+ //
+ eigen3 = ES__begin,
+ slepsc,
+ //
+ ES__end
+};
+;
+
+inline std::string
+name(const ESLibrary library)
+{
+ switch (library) {
+ case ESLibrary::eigen3: {
+ return "Eigen3";
+ }
+ case ESLibrary::slepsc: {
+ return "SLEPSc";
+ }
+ case ESLibrary::ES__end: {
+ }
+ }
+ throw UnexpectedError("Eigenvalue solver library name is not defined!");
+}
+
+template <typename ESEnumType>
+inline ESEnumType
+getESEnumFromName(const std::string& enum_name)
+{
+ using BaseT = std::underlying_type_t<ESEnumType>;
+ for (BaseT enum_value = static_cast<BaseT>(ESEnumType::ES__begin);
+ enum_value < static_cast<BaseT>(ESEnumType::ES__end); ++enum_value) {
+ if (name(ESEnumType{enum_value}) == enum_name) {
+ return ESEnumType{enum_value};
+ }
+ }
+ throw NormalError(std::string{"could not find '"} + enum_name + "' associate type!");
+}
+
+template <typename ESEnumType>
+inline void
+printESEnumListNames(std::ostream& os)
+{
+ using BaseT = std::underlying_type_t<ESEnumType>;
+ for (BaseT enum_value = static_cast<BaseT>(ESEnumType::ES__begin);
+ enum_value < static_cast<BaseT>(ESEnumType::ES__end); ++enum_value) {
+ os << " - " << name(ESEnumType{enum_value}) << '\n';
+ }
+}
+
+class EigenvalueSolverOptions
+{
+ private:
+ ESLibrary m_library = ESLibrary::slepsc;
+
+ double m_epsilon = 1E-6;
+ size_t m_maximum_iteration = 200;
+
+ bool m_verbose = false;
+
+ public:
+ static EigenvalueSolverOptions default_options;
+
+ friend std::ostream& operator<<(std::ostream& os, const EigenvalueSolverOptions& options);
+
+ ESLibrary&
+ library()
+ {
+ return m_library;
+ }
+
+ ESLibrary
+ library() const
+ {
+ return m_library;
+ }
+
+ bool&
+ verbose()
+ {
+ return m_verbose;
+ };
+
+ bool
+ verbose() const
+ {
+ return m_verbose;
+ };
+
+ EigenvalueSolverOptions(const EigenvalueSolverOptions&) = default;
+ EigenvalueSolverOptions(EigenvalueSolverOptions&&) = default;
+
+ EigenvalueSolverOptions() = default;
+ ~EigenvalueSolverOptions() = default;
+};
+
+inline EigenvalueSolverOptions EigenvalueSolverOptions::default_options;
+
+#endif // EIGENVALUE_SOLVER_OPTIONS_HPP
diff --git a/src/language/modules/LinearSolverModule.cpp b/src/language/modules/LinearSolverModule.cpp
index d3e2de0fa..1ff8be36d 100644
--- a/src/language/modules/LinearSolverModule.cpp
+++ b/src/language/modules/LinearSolverModule.cpp
@@ -2,7 +2,8 @@
#include <language/modules/ModuleRepository.hpp>
-#include <algebra/LinearSolver.hpp>
+#include <algebra/EigenvalueSolverOptions.hpp>
+#include <algebra/LinearSolverOptions.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/TypeDescriptor.hpp>
@@ -90,6 +91,43 @@ LinearSolverModule::LinearSolverModule()
return os.str();
}
+ ));
+
+ this->_addBuiltinFunction("setESLibrary", std::function(
+
+ [](const std::string& library_name) -> void {
+ EigenvalueSolverOptions::default_options.library() =
+ getESEnumFromName<ESLibrary>(library_name);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("getESOptions", std::function(
+
+ []() -> std::string {
+ std::ostringstream os;
+ os << rang::fgB::yellow << "Eigenvalue solver options"
+ << rang::style::reset << '\n';
+ os << EigenvalueSolverOptions::default_options;
+ return os.str();
+ }
+
+ ));
+
+ this->_addBuiltinFunction("getESAvailable", std::function(
+
+ []() -> std::string {
+ std::ostringstream os;
+
+ os << rang::fgB::yellow << "Available eigenvalue solver options"
+ << rang::style::reset << '\n';
+ os << rang::fgB::blue << " libraries" << rang::style::reset << '\n';
+
+ printESEnumListNames<ESLibrary>(os);
+
+ return os.str();
+ }
+
));
}
diff --git a/src/utils/checkpointing/Checkpoint.cpp b/src/utils/checkpointing/Checkpoint.cpp
index 5976cb2fb..ed56dcb83 100644
--- a/src/utils/checkpointing/Checkpoint.cpp
+++ b/src/utils/checkpointing/Checkpoint.cpp
@@ -4,6 +4,7 @@
#ifdef PUGS_HAS_HDF5
+#include <algebra/EigenvalueSolverOptions.hpp>
#include <algebra/LinearSolverOptions.hpp>
#include <dev/ParallelChecker.hpp>
#include <language/ast/ASTExecutionStack.hpp>
@@ -21,12 +22,12 @@
#include <utils/HighFivePugsUtils.hpp>
#include <utils/RandomEngine.hpp>
#include <utils/checkpointing/DualMeshTypeHFType.hpp>
+#include <utils/checkpointing/EigenvalueSolverOptionsHFType.hpp>
#include <utils/checkpointing/LinearSolverOptionsHFType.hpp>
#include <utils/checkpointing/ParallelCheckerHFType.hpp>
#include <utils/checkpointing/ResumingManager.hpp>
#include <iostream>
-#include <map>
void
checkpoint()
@@ -122,6 +123,14 @@ checkpoint()
linear_solver_options_default_group.createAttribute("method", default_options.method());
linear_solver_options_default_group.createAttribute("precond", default_options.precond());
}
+ {
+ HighFive::Group eigenvalue_solver_options_default_group =
+ checkpoint.createGroup("singleton/eigenvalue_solver_options_default");
+
+ const EigenvalueSolverOptions& default_options = EigenvalueSolverOptions::default_options;
+
+ eigenvalue_solver_options_default_group.createAttribute("library", default_options.library());
+ }
{
const auto& primal_to_dual_connectivity_info_map =
DualConnectivityManager::instance().primalToDualConnectivityInfoMap();
diff --git a/src/utils/checkpointing/EigenvalueSolverOptionsHFType.hpp b/src/utils/checkpointing/EigenvalueSolverOptionsHFType.hpp
new file mode 100644
index 000000000..32995c1f5
--- /dev/null
+++ b/src/utils/checkpointing/EigenvalueSolverOptionsHFType.hpp
@@ -0,0 +1,15 @@
+#ifndef EIGENVALUE_SOLVER_OPTIONS_HF_TYPE_HPP
+#define EIGENVALUE_SOLVER_OPTIONS_HF_TYPE_HPP
+
+#include <algebra/EigenvalueSolverOptions.hpp>
+#include <utils/HighFivePugsUtils.hpp>
+#include <utils/PugsMacros.hpp>
+
+HighFive::EnumType<ESLibrary> PUGS_INLINE
+create_enum_ESOptions_library_type()
+{
+ return {{"eigen3", ESLibrary::eigen3}, {"slepsc", ESLibrary::slepsc}};
+}
+HIGHFIVE_REGISTER_TYPE(ESLibrary, create_enum_ESOptions_library_type)
+
+#endif // EIGENVALUE_SOLVER_OPTIONS_HF_TYPE_HPP
diff --git a/src/utils/checkpointing/LinearSolverOptionsHFType.hpp b/src/utils/checkpointing/LinearSolverOptionsHFType.hpp
index 94f88fc60..1c94d74a9 100644
--- a/src/utils/checkpointing/LinearSolverOptionsHFType.hpp
+++ b/src/utils/checkpointing/LinearSolverOptionsHFType.hpp
@@ -8,7 +8,7 @@
HighFive::EnumType<LSLibrary> PUGS_INLINE
create_enum_LSOptions_library_type()
{
- return {{"builtin", LSLibrary::builtin}, {"petsc", LSLibrary::petsc}};
+ return {{"builtin", LSLibrary::builtin}, {"eigen3", LSLibrary::eigen3}, {"petsc", LSLibrary::petsc}};
}
HIGHFIVE_REGISTER_TYPE(LSLibrary, create_enum_LSOptions_library_type)
diff --git a/src/utils/checkpointing/Resume.cpp b/src/utils/checkpointing/Resume.cpp
index 9d58e1d9d..effd6b808 100644
--- a/src/utils/checkpointing/Resume.cpp
+++ b/src/utils/checkpointing/Resume.cpp
@@ -4,6 +4,7 @@
#ifdef PUGS_HAS_HDF5
+#include <algebra/EigenvalueSolverOptions.hpp>
#include <algebra/LinearSolverOptions.hpp>
#include <language/ast/ASTExecutionStack.hpp>
#include <language/utils/ASTCheckpointsInfo.hpp>
@@ -14,6 +15,7 @@
#include <utils/ExecutionStatManager.hpp>
#include <utils/HighFivePugsUtils.hpp>
#include <utils/RandomEngine.hpp>
+#include <utils/checkpointing/EigenvalueSolverOptionsHFType.hpp>
#include <utils/checkpointing/LinearSolverOptionsHFType.hpp>
#include <utils/checkpointing/ParallelCheckerHFType.hpp>
#include <utils/checkpointing/ResumingData.hpp>
@@ -81,6 +83,14 @@ resume()
default_options.method() = linear_solver_options_default_group.getAttribute("method").read<LSMethod>();
default_options.precond() = linear_solver_options_default_group.getAttribute("precond").read<LSPrecond>();
}
+ {
+ HighFive::Group eigenvalue_solver_options_default_group =
+ checkpoint.getGroup("singleton/eigenvalue_solver_options_default");
+
+ EigenvalueSolverOptions& default_options = EigenvalueSolverOptions::default_options;
+
+ default_options.library() = eigenvalue_solver_options_default_group.getAttribute("library").read<ESLibrary>();
+ }
checkpointing::ResumingData::instance().readData(checkpoint, p_symbol_table);
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index fb7bf0790..4bb1d339d 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -12,109 +12,321 @@
TEST_CASE("EigenvalueSolver", "[algebra]")
{
- SECTION("Sparse Matrices")
+ SECTION("SLEPc")
{
- SECTION("symmetric system")
+ EigenvalueSolverOptions options;
+ options.library() = ESLibrary::slepsc;
+
+ SECTION("Sparse Matrices")
{
- Array<int> non_zeros(3);
- non_zeros.fill(3);
- CRSMatrixDescriptor<double> S{3, 3, non_zeros};
+ SECTION("symmetric system")
+ {
+ Array<int> non_zeros(3);
+ non_zeros.fill(3);
+ CRSMatrixDescriptor<double> S{3, 3, non_zeros};
- S(0, 0) = 3;
- S(0, 1) = 2;
- S(0, 2) = 4;
+ S(0, 0) = 3;
+ S(0, 1) = 2;
+ S(0, 2) = 4;
- S(1, 0) = 2;
- S(1, 1) = 0;
- S(1, 2) = 2;
+ S(1, 0) = 2;
+ S(1, 1) = 0;
+ S(1, 2) = 2;
- S(2, 0) = 4;
- S(2, 1) = 2;
- S(2, 2) = 3;
+ S(2, 0) = 4;
+ S(2, 1) = 2;
+ S(2, 2) = 3;
- CRSMatrix A{S.getCRSMatrix()};
+ CRSMatrix A{S.getCRSMatrix()};
- SECTION("eigenvalues")
- {
- SmallArray<double> eigenvalues;
+ SECTION("eigenvalues")
+ {
+ SmallArray<double> eigenvalues;
+
+#ifdef PUGS_HAS_SLEPC
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
+ REQUIRE(eigenvalues[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[2] == Catch::Approx(8));
+#else // PUGS_HAS_SLEPC
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
+ "error: SLEPSc is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_SLEPC
+ }
+
+ SECTION("eigenvalues and eigenvectors")
+ {
+ SmallArray<double> eigenvalue_list;
+ std::vector<SmallVector<double>> eigenvector_list;
+
+#ifdef PUGS_HAS_SLEPC
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> P{3};
+ SmallMatrix<double> PT{3};
+ SmallMatrix<double> D{3};
+ D = zero;
+
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ P(i, j) = eigenvector_list[j][i];
+ PT(i, j) = eigenvector_list[i][j];
+ }
+ D(i, i) = eigenvalue_list[i];
+ }
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+#else // PUGS_HAS_SLEPC
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
+ "error: SLEPSc is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_SLEPC
+ }
+
+ SECTION("eigenvalues and transition matrix")
+ {
+ SmallArray<double> eigenvalue_list;
+ SmallMatrix<double> P{};
#ifdef PUGS_HAS_SLEPC
- EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalues);
- REQUIRE(eigenvalues[0] == Catch::Approx(-1));
- REQUIRE(eigenvalues[1] == Catch::Approx(-1));
- REQUIRE(eigenvalues[2] == Catch::Approx(8));
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> D{3};
+ D = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ D(i, i) = eigenvalue_list[i];
+ }
+ SmallMatrix PT = transpose(P);
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
#else // PUGS_HAS_SLEPC
- REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalues),
- "not implemented yet: SLEPc is required to solve eigenvalue problems");
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
+ "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif // PUGS_HAS_SLEPC
+ }
}
+ }
- SECTION("eigenvalues and eigenvectors")
+ SECTION("SmallMatrix")
+ {
+ SECTION("symmetric system")
{
- SmallArray<double> eigenvalue_list;
- std::vector<SmallVector<double>> eigenvector_list;
+ SmallMatrix<double> A{3, 3};
+ A.fill(0);
+ A(0, 0) = 3;
+ A(0, 1) = 2;
+ A(0, 2) = 4;
+
+ A(1, 0) = 2;
+ A(1, 1) = 0;
+ A(1, 2) = 2;
+
+ A(2, 0) = 4;
+ A(2, 1) = 2;
+ A(2, 2) = 3;
+
+ SECTION("eigenvalues")
+ {
+ SmallArray<double> eigenvalues;
+
+#ifdef PUGS_HAS_SLEPC
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
+ REQUIRE(eigenvalues[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[2] == Catch::Approx(8));
+#else // PUGS_HAS_SLEPC
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
+ "error: SLEPSc is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_SLEPC
+ }
+
+ SECTION("eigenvalues and eigenvectors")
+ {
+ SmallArray<double> eigenvalue_list;
+ std::vector<SmallVector<double>> eigenvector_list;
#ifdef PUGS_HAS_SLEPC
- EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
- REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
- REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
- REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
- SmallMatrix<double> P{3};
- SmallMatrix<double> PT{3};
- SmallMatrix<double> D{3};
- D = zero;
+ SmallMatrix<double> P{3};
+ SmallMatrix<double> PT{3};
+ SmallMatrix<double> D{3};
+ D = zero;
- for (size_t i = 0; i < 3; ++i) {
- for (size_t j = 0; j < 3; ++j) {
- P(i, j) = eigenvector_list[j][i];
- PT(i, j) = eigenvector_list[i][j];
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ P(i, j) = eigenvector_list[j][i];
+ PT(i, j) = eigenvector_list[i][j];
+ }
+ D(i, i) = eigenvalue_list[i];
}
- D(i, i) = eigenvalue_list[i];
- }
- SmallMatrix PDPT = P * D * PT;
- for (size_t i = 0; i < 3; ++i) {
- for (size_t j = 0; j < 3; ++j) {
- REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
+ }
}
+#else // PUGS_HAS_SLEPC
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
+ "error: SLEPSc is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_SLEPC
}
+
+ SECTION("eigenvalues and transition matrix")
+ {
+ SmallArray<double> eigenvalue_list;
+ SmallMatrix<double> P{};
+
+#ifdef PUGS_HAS_SLEPC
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> D{3};
+ D = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ D(i, i) = eigenvalue_list[i];
+ }
+ SmallMatrix PT = transpose(P);
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
#else // PUGS_HAS_SLEPC
- REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
- "not implemented yet: SLEPc is required to solve eigenvalue problems");
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
+ "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif // PUGS_HAS_SLEPC
+ }
}
+ }
- SECTION("eigenvalues and passage matrix")
+ SECTION("TinyMatrix")
+ {
+ SECTION("symmetric system")
{
- SmallArray<double> eigenvalue_list;
- SmallMatrix<double> P{};
+ TinyMatrix<3> A = zero;
-#ifdef PUGS_HAS_SLEPC
- EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P);
+ A(0, 0) = 3;
+ A(0, 1) = 2;
+ A(0, 2) = 4;
- REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
- REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
- REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+ A(1, 0) = 2;
+ A(1, 1) = 0;
+ A(1, 2) = 2;
- SmallMatrix<double> D{3};
- D = zero;
- for (size_t i = 0; i < 3; ++i) {
- D(i, i) = eigenvalue_list[i];
+ A(2, 0) = 4;
+ A(2, 1) = 2;
+ A(2, 2) = 3;
+
+ SECTION("eigenvalues")
+ {
+ SmallArray<double> eigenvalues;
+
+#ifdef PUGS_HAS_SLEPC
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
+ REQUIRE(eigenvalues[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[2] == Catch::Approx(8));
+#else // PUGS_HAS_SLEPC
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
+ "error: SLEPSc is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_SLEPC
}
- SmallMatrix PT = transpose(P);
- SmallMatrix PDPT = P * D * PT;
- for (size_t i = 0; i < 3; ++i) {
- for (size_t j = 0; j < 3; ++j) {
- REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
+ SECTION("eigenvalues and eigenvectors")
+ {
+ SmallArray<double> eigenvalue_list;
+ std::vector<SmallVector<double>> eigenvector_list;
+
+#ifdef PUGS_HAS_SLEPC
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> P{3};
+ SmallMatrix<double> PT{3};
+ SmallMatrix<double> D{3};
+ D = zero;
+
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ P(i, j) = eigenvector_list[j][i];
+ PT(i, j) = eigenvector_list[i][j];
+ }
+ D(i, i) = eigenvalue_list[i];
+ }
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
+ }
}
+#else // PUGS_HAS_SLEPC
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
+ "error: SLEPSc is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_SLEPC
}
+
+ SECTION("eigenvalues and transition matrix")
+ {
+ SmallArray<double> eigenvalue_list;
+ SmallMatrix<double> P{};
+
+#ifdef PUGS_HAS_SLEPC
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> D{3};
+ D = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ D(i, i) = eigenvalue_list[i];
+ }
+ SmallMatrix PT = transpose(P);
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
#else // PUGS_HAS_SLEPC
- REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P),
- "not implemented yet: SLEPc is required to solve eigenvalue problems");
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
+ "error: SLEPSc is not linked to pugs. Cannot use it!");
#endif // PUGS_HAS_SLEPC
+ }
}
}
}
diff --git a/tests/test_LinearSolverOptions.cpp b/tests/test_LinearSolverOptions.cpp
index 620959714..7d078a7fb 100644
--- a/tests/test_LinearSolverOptions.cpp
+++ b/tests/test_LinearSolverOptions.cpp
@@ -68,6 +68,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]")
SECTION("library name")
{
REQUIRE(name(LSLibrary::builtin) == "builtin");
+ REQUIRE(name(LSLibrary::eigen3) == "Eigen3");
REQUIRE(name(LSLibrary::petsc) == "PETSc");
REQUIRE_THROWS_WITH(name(LSLibrary::LS__end), "unexpected error: Linear system library name is not defined!");
}
diff --git a/tests/test_checkpointing_HFTypes.cpp b/tests/test_checkpointing_HFTypes.cpp
index 41fc6948c..1a5f02b4a 100644
--- a/tests/test_checkpointing_HFTypes.cpp
+++ b/tests/test_checkpointing_HFTypes.cpp
@@ -3,6 +3,7 @@
#include <utils/checkpointing/DiscreteFunctionTypeHFType.hpp>
#include <utils/checkpointing/DualMeshTypeHFType.hpp>
+#include <utils/checkpointing/EigenvalueSolverOptionsHFType.hpp>
#include <utils/checkpointing/IBoundaryConditionDescriptorHFType.hpp>
#include <utils/checkpointing/IBoundaryDescriptorHFType.hpp>
#include <utils/checkpointing/IInterfaceDescriptorHFType.hpp>
@@ -160,6 +161,7 @@ TEST_CASE("HFTypes", "[utils/checkpointing]")
SECTION("LinearSolverOptionsHFType")
{
file.createAttribute("builtin", LSLibrary::builtin);
+ file.createAttribute("eigen3", LSLibrary::eigen3);
file.createAttribute("petsc", LSLibrary::petsc);
file.createAttribute("cg", LSMethod::cg);
@@ -177,6 +179,7 @@ TEST_CASE("HFTypes", "[utils/checkpointing]")
REQUIRE(file.getAttribute("builtin").read<LSLibrary>() == LSLibrary::builtin);
REQUIRE(file.getAttribute("petsc").read<LSLibrary>() == LSLibrary::petsc);
+ REQUIRE(file.getAttribute("eigen3").read<LSLibrary>() == LSLibrary::eigen3);
REQUIRE(file.getAttribute("cg").read<LSMethod>() == LSMethod::cg);
REQUIRE(file.getAttribute("bicgstab").read<LSMethod>() == LSMethod::bicgstab);
--
GitLab
From 54aed67103b6d89da5111291e9c4f76c53ae77ac Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Thu, 30 Jan 2025 21:31:14 +0100
Subject: [PATCH 5/8] Plug Eigen3 to solve eigenvalue problems for symmetric
matrices
---
src/algebra/EigenvalueSolver.cpp | 151 ++++++++++++++
src/algebra/EigenvalueSolver.hpp | 37 +++-
tests/test_EigenvalueSolver.cpp | 332 +++++++++++++++++++++++++++++--
3 files changed, 498 insertions(+), 22 deletions(-)
diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 5888d1c0b..b17e650af 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -6,6 +6,10 @@
#include <slepc.h>
#endif // PUGS_HAS_SLEPC
+#ifdef PUGS_HAS_EIGEN3
+#include <algebra/Eigen3Utils.hpp>
+#endif // PUGS_HAS_EIGEN3
+
struct EigenvalueSolver::Internals
{
#ifdef PUGS_HAS_SLEPC
@@ -154,6 +158,73 @@ struct EigenvalueSolver::Internals
EPSDestroy(&eps);
}
#endif // PUGS_HAS_SLEPC
+
+#ifdef PUGS_HAS_EIGEN3
+
+ template <typename Eigen3MatrixEmbedderType>
+ static void
+ eigen3ComputeForSymmetricMatrix(const Eigen3MatrixEmbedderType& A, SmallArray<double>& eigenvalues)
+ {
+ using Eigen3MatrixType = typename Eigen3MatrixEmbedderType::Eigen3MatrixType;
+ Eigen::SelfAdjointEigenSolver<Eigen3MatrixType> eigen_solver;
+
+ eigen_solver.compute(A.matrix(), Eigen::EigenvaluesOnly);
+
+ eigenvalues = SmallArray<double>(A.numberOfRows());
+ for (size_t i = 0; i < eigenvalues.size(); ++i) {
+ eigenvalues[i] = eigen_solver.eigenvalues()[i];
+ }
+ }
+
+ template <typename Eigen3MatrixEmbedderType>
+ static void
+ eigen3ComputeForSymmetricMatrix(const Eigen3MatrixEmbedderType& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors)
+ {
+ using Eigen3MatrixType = typename Eigen3MatrixEmbedderType::Eigen3MatrixType;
+ Eigen::SelfAdjointEigenSolver<Eigen3MatrixType> eigen_solver;
+
+ eigen_solver.compute(A.matrix(), Eigen::ComputeEigenvectors);
+
+ eigenvalues = SmallArray<double>(A.numberOfRows());
+ for (size_t i = 0; i < eigenvalues.size(); ++i) {
+ eigenvalues[i] = eigen_solver.eigenvalues()[i];
+ }
+
+ eigenvectors.resize(eigenvalues.size());
+ for (size_t i = 0; i < eigenvalues.size(); ++i) {
+ eigenvectors[i] = SmallVector<double>{eigenvalues.size()};
+ for (size_t j = 0; j < eigenvalues.size(); ++j) {
+ eigenvectors[i][j] = eigen_solver.eigenvectors().coeff(j, i);
+ }
+ }
+ }
+
+ template <typename Eigen3MatrixEmbedderType>
+ static void
+ eigen3ComputeForSymmetricMatrix(const Eigen3MatrixEmbedderType& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P)
+ {
+ using Eigen3MatrixType = typename Eigen3MatrixEmbedderType::Eigen3MatrixType;
+ Eigen::SelfAdjointEigenSolver<Eigen3MatrixType> eigen_solver;
+
+ eigen_solver.compute(A.matrix(), Eigen::ComputeEigenvectors);
+
+ eigenvalues = SmallArray<double>(A.numberOfRows());
+ for (size_t i = 0; i < eigenvalues.size(); ++i) {
+ eigenvalues[i] = eigen_solver.eigenvalues()[i];
+ }
+
+ P = SmallMatrix<double>(eigenvalues.size(), eigenvalues.size());
+ for (size_t i = 0; i < eigenvalues.size(); ++i) {
+ for (size_t j = 0; j < eigenvalues.size(); ++j) {
+ P(i, j) = eigen_solver.eigenvectors().coeff(i, j);
+ }
+ }
+ }
+#endif // PUGS_HAS_EIGEN3
};
#ifdef PUGS_HAS_SLEPC
@@ -235,6 +306,86 @@ EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<doub
#endif // PUGS_HAS_SLEPC
+#ifdef PUGS_HAS_EIGEN3
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues)
+{
+ Eigen3SparseMatrixEmbedder embedded_A{A};
+ Internals::eigen3ComputeForSymmetricMatrix(embedded_A, eigenvalues);
+}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues)
+{
+ Internals::eigen3ComputeForSymmetricMatrix(Eigen3DenseMatrixEmbedder{A}, eigenvalues);
+}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors)
+{
+ Internals::eigen3ComputeForSymmetricMatrix(Eigen3SparseMatrixEmbedder{A}, eigenvalues, eigenvectors);
+}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors)
+{
+ Internals::eigen3ComputeForSymmetricMatrix(Eigen3DenseMatrixEmbedder{A}, eigenvalues, eigenvectors);
+}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P)
+{
+ Internals::eigen3ComputeForSymmetricMatrix(Eigen3SparseMatrixEmbedder{A}, eigenvalues, P);
+}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P)
+{
+ Internals::eigen3ComputeForSymmetricMatrix(Eigen3DenseMatrixEmbedder{A}, eigenvalues, P);
+}
+
+#else // PUGS_HAS_EIGEN3
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&, SmallArray<double>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>&,
+ SmallArray<double>&,
+ std::vector<SmallVector<double>>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&,
+ SmallArray<double>&,
+ std::vector<SmallVector<double>>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>&, SmallArray<double>&, SmallMatrix<double>&)
+{}
+
+void
+EigenvalueSolver::_eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>&,
+ SmallArray<double>&,
+ SmallMatrix<double>&)
+{}
+
+#endif // PUGS_HAS_EIGEN3
+
bool
EigenvalueSolver::hasLibrary(const ESLibrary library)
{
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index eebd536e4..6567adcdc 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -38,6 +38,26 @@ class EigenvalueSolver
SmallArray<double>& eigenvalues,
SmallMatrix<double>& P);
+ void _eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues);
+
+ void _eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues);
+
+ void _eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors);
+
+ void _eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors);
+
+ void _eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P);
+
+ void _eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P);
+
public:
static bool hasLibrary(ESLibrary library);
static void checkHasLibrary(const ESLibrary library);
@@ -46,11 +66,12 @@ class EigenvalueSolver
void
computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues)
{
- Assert((eigenvalues.size() - A.numberOfRows() == 0) and A.isSquare());
+ Assert(A.isSquare());
switch (m_options.library()) {
case ESLibrary::eigen3: {
- throw NotImplementedError("Eigen3");
+ this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues);
+ break;
}
case ESLibrary::slepsc: {
this->_slepscComputeForSymmetricMatrix(A, eigenvalues);
@@ -68,12 +89,12 @@ class EigenvalueSolver
SmallArray<double>& eigenvalues,
std::vector<SmallVector<double>>& eigenvectors)
{
- Assert((eigenvalues.size() - A.numberOfRows() == 0) and (A.numberOfRows() - eigenvectors.size() == 0) and
- A.isSquare());
+ Assert(A.isSquare());
switch (m_options.library()) {
case ESLibrary::eigen3: {
- throw NotImplementedError("Eigen3");
+ this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
+ break;
}
case ESLibrary::slepsc: {
this->_slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
@@ -89,12 +110,12 @@ class EigenvalueSolver
void
computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues, SmallMatrix<double>& P)
{
- Assert((eigenvalues.size() - A.numberOfRows() == 0) and (A.numberOfRows() - P.numberOfRows() == 0) and
- A.isSquare() and P.isSquare());
+ Assert(A.isSquare());
switch (m_options.library()) {
case ESLibrary::eigen3: {
- throw NotImplementedError("Eigen3");
+ this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues, P);
+ break;
}
case ESLibrary::slepsc: {
this->_slepscComputeForSymmetricMatrix(A, eigenvalues, P);
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 4bb1d339d..325f3e39b 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -12,14 +12,14 @@
TEST_CASE("EigenvalueSolver", "[algebra]")
{
- SECTION("SLEPc")
+ SECTION("symmetric system")
{
- EigenvalueSolverOptions options;
- options.library() = ESLibrary::slepsc;
-
- SECTION("Sparse Matrices")
+ SECTION("SLEPc")
{
- SECTION("symmetric system")
+ EigenvalueSolverOptions options;
+ options.library() = ESLibrary::slepsc;
+
+ SECTION("Sparse Matrices")
{
Array<int> non_zeros(3);
non_zeros.fill(3);
@@ -122,11 +122,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
#endif // PUGS_HAS_SLEPC
}
}
- }
- SECTION("SmallMatrix")
- {
- SECTION("symmetric system")
+ SECTION("SmallMatrix")
{
SmallMatrix<double> A{3, 3};
A.fill(0);
@@ -225,11 +222,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
#endif // PUGS_HAS_SLEPC
}
}
- }
- SECTION("TinyMatrix")
- {
- SECTION("symmetric system")
+ SECTION("TinyMatrix")
{
TinyMatrix<3> A = zero;
@@ -329,5 +323,315 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
}
}
}
+
+ SECTION("Eigen3")
+ {
+ EigenvalueSolverOptions options;
+ options.library() = ESLibrary::eigen3;
+
+ SECTION("Sparse Matrices")
+ {
+ Array<int> non_zeros(3);
+ non_zeros.fill(3);
+ CRSMatrixDescriptor<double> S{3, 3, non_zeros};
+
+ S(0, 0) = 3;
+ S(0, 1) = 2;
+ S(0, 2) = 4;
+
+ S(1, 0) = 2;
+ S(1, 1) = 0;
+ S(1, 2) = 2;
+
+ S(2, 0) = 4;
+ S(2, 1) = 2;
+ S(2, 2) = 3;
+
+ CRSMatrix A{S.getCRSMatrix()};
+
+ SECTION("eigenvalues")
+ {
+ SmallArray<double> eigenvalues;
+
+#ifdef PUGS_HAS_EIGEN3
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
+ REQUIRE(eigenvalues[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[2] == Catch::Approx(8));
+#else // PUGS_HAS_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and eigenvectors")
+ {
+ SmallArray<double> eigenvalue_list;
+ std::vector<SmallVector<double>> eigenvector_list;
+
+#ifdef PUGS_HAS_EIGEN3
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> P{3};
+ SmallMatrix<double> PT{3};
+ SmallMatrix<double> D{3};
+ D = zero;
+
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ P(i, j) = eigenvector_list[j][i];
+ PT(i, j) = eigenvector_list[i][j];
+ }
+ D(i, i) = eigenvalue_list[i];
+ }
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+#else // PUGS_HAS_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and transition matrix")
+ {
+ SmallArray<double> eigenvalue_list;
+ SmallMatrix<double> P{};
+
+#ifdef PUGS_HAS_EIGEN3
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> D{3};
+ D = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ D(i, i) = eigenvalue_list[i];
+ }
+ SmallMatrix PT = transpose(P);
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+#else // PUGS_HAS_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+ }
+
+ SECTION("SmallMatrix")
+ {
+ SmallMatrix<double> A{3, 3};
+ A.fill(0);
+ A(0, 0) = 3;
+ A(0, 1) = 2;
+ A(0, 2) = 4;
+
+ A(1, 0) = 2;
+ A(1, 1) = 0;
+ A(1, 2) = 2;
+
+ A(2, 0) = 4;
+ A(2, 1) = 2;
+ A(2, 2) = 3;
+
+ SECTION("eigenvalues")
+ {
+ SmallArray<double> eigenvalues;
+
+#ifdef PUGS_HAS_EIGEN3
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
+ REQUIRE(eigenvalues[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[2] == Catch::Approx(8));
+#else // PUGS_HAS_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and eigenvectors")
+ {
+ SmallArray<double> eigenvalue_list;
+ std::vector<SmallVector<double>> eigenvector_list;
+
+#ifdef PUGS_HAS_EIGEN3
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> P{3};
+ SmallMatrix<double> PT{3};
+ SmallMatrix<double> D{3};
+ D = zero;
+
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ P(i, j) = eigenvector_list[j][i];
+ PT(i, j) = eigenvector_list[i][j];
+ }
+ D(i, i) = eigenvalue_list[i];
+ }
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+#else // PUGS_HAS_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and transition matrix")
+ {
+ SmallArray<double> eigenvalue_list;
+ SmallMatrix<double> P{};
+
+#ifdef PUGS_HAS_EIGEN3
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> D{3};
+ D = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ D(i, i) = eigenvalue_list[i];
+ }
+ SmallMatrix PT = transpose(P);
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+#else // PUGS_HAS_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+ }
+
+ SECTION("TinyMatrix")
+ {
+ TinyMatrix<3> A = zero;
+
+ A(0, 0) = 3;
+ A(0, 1) = 2;
+ A(0, 2) = 4;
+
+ A(1, 0) = 2;
+ A(1, 1) = 0;
+ A(1, 2) = 2;
+
+ A(2, 0) = 4;
+ A(2, 1) = 2;
+ A(2, 2) = 3;
+
+ SECTION("eigenvalues")
+ {
+ SmallArray<double> eigenvalues;
+
+#ifdef PUGS_HAS_EIGEN3
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues);
+ REQUIRE(eigenvalues[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[2] == Catch::Approx(8));
+#else // PUGS_HAS_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and eigenvectors")
+ {
+ SmallArray<double> eigenvalue_list;
+ std::vector<SmallVector<double>> eigenvector_list;
+
+#ifdef PUGS_HAS_EIGEN3
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> P{3};
+ SmallMatrix<double> PT{3};
+ SmallMatrix<double> D{3};
+ D = zero;
+
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ P(i, j) = eigenvector_list[j][i];
+ PT(i, j) = eigenvector_list[i][j];
+ }
+ D(i, i) = eigenvalue_list[i];
+ }
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+#else // PUGS_HAS_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("eigenvalues and transition matrix")
+ {
+ SmallArray<double> eigenvalue_list;
+ SmallMatrix<double> P{};
+
+#ifdef PUGS_HAS_EIGEN3
+ EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ SmallMatrix<double> D{3};
+ D = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ D(i, i) = eigenvalue_list[i];
+ }
+ SmallMatrix PT = transpose(P);
+
+ SmallMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - A(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+#else // PUGS_HAS_EIGEN3
+ REQUIRE_THROWS_WITH(EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalue_list, P),
+ "error: Eigen3 is not linked to pugs. Cannot use it!");
+#endif // PUGS_HAS_EIGEN3
+ }
+ }
+ }
}
}
--
GitLab
From eca2ebc50372138653ee9af2c945378e64ce85b2 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Fri, 31 Jan 2025 08:29:52 +0100
Subject: [PATCH 6/8] Add few Eigen3 installation instructions
---
README.md | 10 ++++++++++
1 file changed, 10 insertions(+)
diff --git a/README.md b/README.md
index c859bc694..473b27b07 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` |
--
GitLab
From 02dda00424f3a631e1e0e5dcf401bfa67cef3f0b Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Fri, 31 Jan 2025 09:01:44 +0100
Subject: [PATCH 7/8] Fix compilation warning
---
src/algebra/EigenvalueSolverOptions.hpp | 1 -
1 file changed, 1 deletion(-)
diff --git a/src/algebra/EigenvalueSolverOptions.hpp b/src/algebra/EigenvalueSolverOptions.hpp
index d9be86be0..2c74ce41e 100644
--- a/src/algebra/EigenvalueSolverOptions.hpp
+++ b/src/algebra/EigenvalueSolverOptions.hpp
@@ -14,7 +14,6 @@ enum class ESLibrary : int8_t
//
ES__end
};
-;
inline std::string
name(const ESLibrary library)
--
GitLab
From 8827f99e0c6615444ba737cc2ade9c26ea7e8014 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Fri, 31 Jan 2025 09:09:13 +0100
Subject: [PATCH 8/8] Fix another bunch of warnings
---
src/algebra/EigenvalueSolverOptions.hpp | 17 -----------------
1 file changed, 17 deletions(-)
diff --git a/src/algebra/EigenvalueSolverOptions.hpp b/src/algebra/EigenvalueSolverOptions.hpp
index 2c74ce41e..1c047ef3c 100644
--- a/src/algebra/EigenvalueSolverOptions.hpp
+++ b/src/algebra/EigenvalueSolverOptions.hpp
@@ -61,11 +61,6 @@ class EigenvalueSolverOptions
private:
ESLibrary m_library = ESLibrary::slepsc;
- double m_epsilon = 1E-6;
- size_t m_maximum_iteration = 200;
-
- bool m_verbose = false;
-
public:
static EigenvalueSolverOptions default_options;
@@ -83,18 +78,6 @@ class EigenvalueSolverOptions
return m_library;
}
- bool&
- verbose()
- {
- return m_verbose;
- };
-
- bool
- verbose() const
- {
- return m_verbose;
- };
-
EigenvalueSolverOptions(const EigenvalueSolverOptions&) = default;
EigenvalueSolverOptions(EigenvalueSolverOptions&&) = default;
--
GitLab