diff --git a/src/algebra/Eigen3Utils.hpp b/src/algebra/Eigen3Utils.hpp
index 14434fe940878366171fc7b264e1aa50419124ab..5224a4b0bb10dc252b0f299cc9cbbbfe9fa87712 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 7504ced5f7390e3d69c4587c0d2dbf1f1b153517..1ebb59663405883aac2405979ff60bdb7ccb3e7e 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 01c71cbea82041747ee137ca76d15a17de89c170..4d01df2f0605e2bd141c1de3d95f28776f3178a0 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
+ }
}
}
}