diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index e11e60ae0c8b72ef111e1bc64c69ffe931d98621..1965a9966cafdf181b366feb1e41909077edd8ac 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -24,9 +24,11 @@ struct LinearSolver::Internals
return false;
#endif
}
+ // LCOV_EXCL_START
default: {
throw UnexpectedError("Linear system library (" + ::name(library) + ") was not set!");
}
+ // LCOV_EXCL_STOP
}
}
@@ -34,7 +36,9 @@ struct LinearSolver::Internals
checkHasLibrary(const LSLibrary library)
{
if (not hasLibrary(library)) {
+ // LCOV_EXCL_START
throw NormalError(::name(library) + " is not linked to pugs. Cannot use it!");
+ // LCOV_EXCL_STOP
}
}
@@ -64,9 +68,11 @@ struct LinearSolver::Internals
case LSMethod::choleski: {
break;
}
+ // LCOV_EXCL_START
default: {
throw NormalError(name(method) + " is not a builtin linear solver!");
}
+ // LCOV_EXCL_STOP
}
}
@@ -94,9 +100,11 @@ struct LinearSolver::Internals
case LSPrecond::incomplete_LU: {
break;
}
+ // LCOV_EXCL_START
default: {
throw NormalError(name(precond) + " is not a PETSc preconditioner!");
}
+ // LCOV_EXCL_STOP
}
}
@@ -114,9 +122,11 @@ struct LinearSolver::Internals
checkPETScPrecond(options.precond());
break;
}
+ // LCOV_EXCL_START
default: {
throw UnexpectedError("undefined options compatibility for this library (" + ::name(options.library()) + ")!");
}
+ // LCOV_EXCL_STOP
}
}
@@ -127,7 +137,9 @@ struct LinearSolver::Internals
const LinearSolverOptions& options)
{
if (options.precond() != LSPrecond::none) {
+ // LCOV_EXCL_START
throw UnexpectedError("builtin linear solver do not allow any preconditioner!");
+ // LCOV_EXCL_STOP
}
switch (options.method()) {
case LSMethod::cg: {
@@ -138,14 +150,15 @@ struct LinearSolver::Internals
BiCGStab{A, x, b, options.epsilon(), options.maximumIteration(), options.verbose()};
break;
}
+ // LCOV_EXCL_START
default: {
throw NotImplementedError("undefined builtin method: " + name(options.method()));
}
+ // LCOV_EXCL_STOP
}
}
#ifdef PUGS_HAS_PETSC
-
static int
petscMonitor(KSP, int i, double residu, void*)
{
@@ -238,9 +251,11 @@ struct LinearSolver::Internals
direct_solver = true;
break;
}
+ // LCOV_EXCL_START
default: {
throw UnexpectedError("unexpected method: " + name(options.method()));
}
+ // LCOV_EXCL_STOP
}
if (not direct_solver) {
@@ -265,9 +280,11 @@ struct LinearSolver::Internals
PCSetType(pc, PCNONE);
break;
}
+ // LCOV_EXCL_START
default: {
throw UnexpectedError("unexpected preconditioner: " + name(options.precond()));
}
+ // LCOV_EXCL_STOP
}
}
if (options.verbose()) {
@@ -278,7 +295,9 @@ struct LinearSolver::Internals
MatDestroy(&petscMat);
}
-#else
+#else // PUGS_HAS_PETSC
+
+ // LCOV_EXCL_START
static void
petscSolveLocalSystem(const CRSMatrix<double, size_t>&,
Vector<double>&,
@@ -288,6 +307,7 @@ struct LinearSolver::Internals
checkHasLibrary(LSLibrary::petsc);
throw UnexpectedError("unexpected situation should not reach this point!");
}
+ // LCOV_EXCL_STOP
#endif // PUGS_HAS_PETSC
};
@@ -298,6 +318,12 @@ LinearSolver::hasLibrary(LSLibrary library) const
return Internals::hasLibrary(library);
}
+void
+LinearSolver::checkOptions(const LinearSolverOptions& options) const
+{
+ Internals::checkOptions(options);
+}
+
void
LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<double>& b)
{
@@ -306,13 +332,18 @@ LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double
Internals::builtinSolveLocalSystem(A, x, b, m_options);
break;
}
+ // LCOV_EXCL_START
case LSLibrary::petsc: {
+ // not covered since if PETSc is not linked this point is
+ // unreachable: LinearSolver throws an exception at construction
+ // in this case.
Internals::petscSolveLocalSystem(A, x, b, m_options);
break;
}
default: {
throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems for sparse matrices");
}
+ // LCOV_EXCL_STOP
}
}
diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
index b72b09b0f6a3255f914c96ee3df7717923931294..3ed5660bd9ce51c1e5ba5de3fe00c595aa96f1e7 100644
--- a/src/algebra/LinearSolver.hpp
+++ b/src/algebra/LinearSolver.hpp
@@ -19,6 +19,7 @@ class LinearSolver
public:
bool hasLibrary(LSLibrary library) const;
+ void checkOptions(const LinearSolverOptions& options) const;
void solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<double>& b);
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 61a2b661986a3b0ab2d8484bdd38779e197b28c0..aaf4bfbad904a6466b4442b0d99f35a41840ddb3 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -73,6 +73,7 @@ add_executable (unit_tests
test_IncDecExpressionProcessor.cpp
test_INodeProcessor.cpp
test_ItemType.cpp
+ test_LinearSolver.cpp
test_LinearSolverOptions.cpp
test_ListAffectationProcessor.cpp
test_MathModule.cpp
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f56bb59ed1d6b7df2e81d8e4076f0ee23bb27ec6
--- /dev/null
+++ b/tests/test_LinearSolver.cpp
@@ -0,0 +1,527 @@
+#include <catch2/catch.hpp>
+
+#include <utils/pugs_config.hpp>
+
+#include <algebra/LinearSolver.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("LinearSolver", "[algebra]")
+{
+ SECTION("check has library")
+ {
+ LinearSolver linear_solver;
+
+ REQUIRE(linear_solver.hasLibrary(LSLibrary::builtin) == true);
+
+#ifdef PUGS_HAS_PETSC
+ REQUIRE(linear_solver.hasLibrary(LSLibrary::petsc) == true);
+#else // PUGS_HAS_PETSC
+ REQUIRE(linear_solver.hasLibrary(LSLibrary::petsc) == false);
+#endif // PUGS_HAS_PETSC
+ }
+
+ SECTION("check linear solver building")
+ {
+ LinearSolverOptions options;
+
+ SECTION("builtin")
+ {
+ SECTION("builtin methods")
+ {
+ options.library() = LSLibrary::builtin;
+ options.precond() = LSPrecond::none;
+
+ options.method() = LSMethod::cg;
+ REQUIRE_NOTHROW(LinearSolver{options});
+
+ options.method() = LSMethod::bicgstab;
+ REQUIRE_NOTHROW(LinearSolver{options});
+
+ options.method() = LSMethod::bicgstab2;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: BICGStab2 is not a builtin linear solver!");
+
+ options.method() = LSMethod::lu;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: LU is not a builtin linear solver!");
+
+ options.method() = LSMethod::choleski;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Choleski is not a builtin linear solver!");
+
+ options.method() = LSMethod::gmres;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: GMRES is not a builtin linear solver!");
+ }
+
+ SECTION("builtin precond")
+ {
+ options.library() = LSLibrary::builtin;
+ options.method() = LSMethod::cg;
+
+ options.precond() = LSPrecond::none;
+ REQUIRE_NOTHROW(LinearSolver{options});
+
+ options.precond() = LSPrecond::diagonal;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: diagonal is not a builtin preconditioner!");
+
+ options.precond() = LSPrecond::incomplete_LU;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: ILU is not a builtin preconditioner!");
+
+ options.precond() = LSPrecond::incomplete_choleski;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: ICholeski is not a builtin preconditioner!");
+
+ options.precond() = LSPrecond::amg;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not a builtin preconditioner!");
+ }
+ }
+
+ SECTION("PETSc")
+ {
+ LinearSolverOptions always_valid;
+ always_valid.library() = LSLibrary::builtin;
+ always_valid.method() = LSMethod::cg;
+ always_valid.precond() = LSPrecond::none;
+
+ LinearSolver linear_solver{always_valid};
+
+ SECTION("PETSc methods")
+ {
+ options.library() = LSLibrary::petsc;
+ 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::choleski;
+ REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+ options.method() = LSMethod::gmres;
+ REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+ }
+
+ SECTION("builtin precond")
+ {
+ options.library() = LSLibrary::petsc;
+ 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_choleski;
+ REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+
+ options.precond() = LSPrecond::amg;
+ REQUIRE_NOTHROW(linear_solver.checkOptions(options));
+ }
+ }
+
+#ifndef PUGS_HAS_PETSC
+ SECTION("not linked PETSc")
+ {
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
+#endif // PUGS_HAS_PETSC
+ }
+
+ SECTION("check linear solvers")
+ {
+ SECTION("symmetric system")
+ {
+ SparseMatrixDescriptor S{5};
+ S(0, 0) = 2;
+ S(0, 1) = -1;
+
+ S(1, 0) = -1;
+ S(1, 1) = 2;
+ S(1, 2) = -1;
+
+ S(2, 1) = -1;
+ S(2, 2) = 2;
+ S(2, 3) = -1;
+
+ S(3, 2) = -1;
+ S(3, 3) = 2;
+ S(3, 4) = -1;
+
+ S(4, 3) = -1;
+ S(4, 4) = 2;
+
+ CRSMatrix A{S};
+
+ Vector<const double> x_exact = [] {
+ Vector<double> y{5};
+ y[0] = 1;
+ y[1] = 3;
+ y[2] = 2;
+ y[3] = 4;
+ y[4] = 5;
+ return y;
+ }();
+
+ Vector<double> b = A * x_exact;
+
+ SECTION("builtin")
+ {
+ SECTION("CG no preconditioner")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::builtin;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+ }
+
+ 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::none;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+
+ SECTION("CG Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+
+ SECTION("CG ICholeski")
+ {
+ options.precond() = LSPrecond::incomplete_choleski;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+
+ SECTION("CG AMG")
+ {
+ options.precond() = LSPrecond::amg;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+ }
+
+ SECTION("Choleski")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::choleski;
+ options.precond() = LSPrecond::none;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+
+#else // PUGS_HAS_PETSC
+ SECTION("PETSc not linked")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
+#endif // PUGS_HAS_PETSC
+ }
+ }
+
+ SECTION("none symmetric system")
+ {
+ SparseMatrixDescriptor S{5};
+ S(0, 0) = 2;
+ S(0, 1) = -1;
+
+ S(1, 0) = -0.2;
+ S(1, 1) = 2;
+ S(1, 2) = -1;
+
+ S(2, 1) = -1;
+ S(2, 2) = 4;
+ S(2, 3) = -2;
+
+ S(3, 2) = -1;
+ S(3, 3) = 2;
+ S(3, 4) = -0.1;
+
+ S(4, 3) = 1;
+ S(4, 4) = 3;
+
+ CRSMatrix A{S};
+
+ Vector<const double> x_exact = [] {
+ Vector<double> y{5};
+ y[0] = 1;
+ y[1] = 3;
+ y[2] = 2;
+ y[3] = 4;
+ y[4] = 5;
+ return y;
+ }();
+
+ Vector<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;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+ }
+
+ SECTION("PETSc")
+ {
+#ifdef PUGS_HAS_PETSC
+
+ SECTION("BICGStab")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
+
+ SECTION("BICGStab no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((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;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab2 Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((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;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+
+ SECTION("GMRES Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+
+ SECTION("GMRES ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+ }
+
+ SECTION("LU")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::lu;
+ options.precond() = LSPrecond::none;
+
+ Vector<double> x{5};
+ x = 0;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ Vector error = x - x_exact;
+ REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+ }
+
+#else // PUGS_HAS_PETSC
+ SECTION("PETSc not linked")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
+#endif // PUGS_HAS_PETSC
+ }
+ }
+ }
+}