Skip to content
Snippets Groups Projects
Commit b1ebe1e7 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Plug a bunch of Eigen3 linear solvers

parent c72934b2
No related branches found
No related tags found
1 merge request!201Feature/eigen3
Show whitespace changes
Inline Side-by-side
src/algebra/Eigen3Utils.hpp
+ 10
− 8
View file @ b1ebe1e7
...@@ -14,10 +14,11 @@ ...@@ -14,10 +14,11 @@
class Eigen3DenseMatrixEmbedder class Eigen3DenseMatrixEmbedder
{ {
public: 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: private:
MatrixType m_eigen_matrix; Eigen3MapMatrixType m_eigen_matrix;
Eigen3DenseMatrixEmbedder(const size_t nb_rows, const size_t nb_columns, const double* A) Eigen3DenseMatrixEmbedder(const size_t nb_rows, const size_t nb_columns, const double* A)
: m_eigen_matrix{A, int(nb_rows), int(nb_columns)} : m_eigen_matrix{A, int(nb_rows), int(nb_columns)}
...@@ -39,14 +40,14 @@ class Eigen3DenseMatrixEmbedder ...@@ -39,14 +40,14 @@ class Eigen3DenseMatrixEmbedder
} }
PUGS_INLINE PUGS_INLINE
MatrixType& Eigen3MapMatrixType&
matrix() matrix()
{ {
return m_eigen_matrix; return m_eigen_matrix;
} }
PUGS_INLINE PUGS_INLINE
const MatrixType& const Eigen3MapMatrixType&
matrix() const matrix() const
{ {
return m_eigen_matrix; return m_eigen_matrix;
...@@ -66,10 +67,11 @@ class Eigen3DenseMatrixEmbedder ...@@ -66,10 +67,11 @@ class Eigen3DenseMatrixEmbedder
class Eigen3SparseMatrixEmbedder class Eigen3SparseMatrixEmbedder
{ {
public: 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: private:
MatrixType m_eigen_matrix; Eigen3MapMatrixType m_eigen_matrix;
public: public:
PUGS_INLINE PUGS_INLINE
...@@ -87,14 +89,14 @@ class Eigen3SparseMatrixEmbedder ...@@ -87,14 +89,14 @@ class Eigen3SparseMatrixEmbedder
} }
PUGS_INLINE PUGS_INLINE
MatrixType& Eigen3MapMatrixType&
matrix() matrix()
{ {
return m_eigen_matrix; return m_eigen_matrix;
} }
PUGS_INLINE PUGS_INLINE
const MatrixType& const Eigen3MapMatrixType&
matrix() const matrix() const
{ {
return m_eigen_matrix; return m_eigen_matrix;
......
src/algebra/LinearSolver.cpp
+ 180
− 17
View file @ b1ebe1e7
...@@ -3,8 +3,13 @@ ...@@ -3,8 +3,13 @@
#include <algebra/BiCGStab.hpp> #include <algebra/BiCGStab.hpp>
#include <algebra/CG.hpp> #include <algebra/CG.hpp>
#include <algebra/Eigen3Utils.hpp>
#include <algebra/PETScUtils.hpp> #include <algebra/PETScUtils.hpp>
#ifdef PUGS_HAS_EIGEN3
#include <eigen3/unsupported/Eigen/IterativeSolvers>
#endif // PUGS_HAS_EIGEN3
#ifdef PUGS_HAS_PETSC #ifdef PUGS_HAS_PETSC
#include <petsc.h> #include <petsc.h>
#endif // PUGS_HAS_PETSC #endif // PUGS_HAS_PETSC
...@@ -69,13 +74,10 @@ struct LinearSolver::Internals ...@@ -69,13 +74,10 @@ struct LinearSolver::Internals
{ {
switch (method) { switch (method) {
case LSMethod::cg: case LSMethod::cg:
#warning clean-up case LSMethod::bicgstab:
// case LSMethod::bicgstab: case LSMethod::gmres:
// case LSMethod::bicgstab2: case LSMethod::lu:
// case LSMethod::gmres: case LSMethod::cholesky: {
// case LSMethod::lu:
// case LSMethod::cholesky:
{
break; break;
} }
// LCOV_EXCL_START // LCOV_EXCL_START
...@@ -124,12 +126,9 @@ struct LinearSolver::Internals ...@@ -124,12 +126,9 @@ struct LinearSolver::Internals
{ {
switch (precond) { switch (precond) {
case LSPrecond::none: case LSPrecond::none:
#warning clean-up case LSPrecond::diagonal:
// case LSPrecond::amg: case LSPrecond::incomplete_cholesky:
// case LSPrecond::diagonal: case LSPrecond::incomplete_LU: {
// case LSPrecond::incomplete_cholesky:
// case LSPrecond::incomplete_LU:
{
break; break;
} }
// LCOV_EXCL_START // LCOV_EXCL_START
...@@ -215,6 +214,163 @@ struct LinearSolver::Internals ...@@ -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 #ifdef PUGS_HAS_PETSC
static int static int
petscMonitor(KSP, int i, double residu, void*) petscMonitor(KSP, int i, double residu, void*)
...@@ -357,6 +513,13 @@ struct LinearSolver::Internals ...@@ -357,6 +513,13 @@ struct LinearSolver::Internals
break; break;
} }
// LCOV_EXCL_START // 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: { case LSLibrary::petsc: {
// not covered since if PETSc is not linked this point is // not covered since if PETSc is not linked this point is
// unreachable: LinearSolver throws an exception at construction // unreachable: LinearSolver throws an exception at construction
......
tests/test_LinearSolver.cpp
+ 496
− 258
View file @ b1ebe1e7
...@@ -167,9 +167,6 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -167,9 +167,6 @@ TEST_CASE("LinearSolver", "[algebra]")
options.method() = LSMethod::bicgstab; options.method() = LSMethod::bicgstab;
REQUIRE_NOTHROW(linear_solver.checkOptions(options)); REQUIRE_NOTHROW(linear_solver.checkOptions(options));
options.method() = LSMethod::bicgstab2;
REQUIRE_NOTHROW(linear_solver.checkOptions(options));
options.method() = LSMethod::lu; options.method() = LSMethod::lu;
REQUIRE_NOTHROW(linear_solver.checkOptions(options)); REQUIRE_NOTHROW(linear_solver.checkOptions(options));
...@@ -196,9 +193,6 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -196,9 +193,6 @@ TEST_CASE("LinearSolver", "[algebra]")
options.precond() = LSPrecond::incomplete_cholesky; options.precond() = LSPrecond::incomplete_cholesky;
REQUIRE_NOTHROW(linear_solver.checkOptions(options)); REQUIRE_NOTHROW(linear_solver.checkOptions(options));
options.precond() = LSPrecond::amg;
REQUIRE_NOTHROW(linear_solver.checkOptions(options));
} }
} }
...@@ -302,48 +296,44 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -302,48 +296,44 @@ TEST_CASE("LinearSolver", "[algebra]")
Vector error = x - x_exact; Vector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, 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}; SECTION("CG Diagonal")
// x = zero; {
options.precond() = LSPrecond::diagonal;
// LinearSolver solver{options}; Vector<double> x{5};
x = zero;
// solver.solveLocalSystem(A, x, b); LinearSolver solver{options};
// Vector error = x - x_exact;
// REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
// }
// SECTION("CG ICholesky") solver.solveLocalSystem(A, x, b);
// { Vector error = x - x_exact;
// options.precond() = LSPrecond::incomplete_cholesky; REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
// Vector<double> x{5}; SECTION("CG ICholesky")
// x = zero; {
options.precond() = LSPrecond::incomplete_cholesky;
// LinearSolver solver{options}; Vector<double> x{5};
x = zero;
// solver.solveLocalSystem(A, x, b); LinearSolver solver{options};
// Vector error = x - x_exact;
// REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
// }
// SECTION("CG AMG") solver.solveLocalSystem(A, x, b);
// { Vector error = x - x_exact;
// options.precond() = LSPrecond::amg; REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
// Vector<double> x{5}; SECTION("CG AMG")
// x = zero; {
options.precond() = LSPrecond::amg;
// LinearSolver solver{options}; Vector<double> x{5};
x = zero;
// solver.solveLocalSystem(A, x, b); REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not an Eigen3 preconditioner!");
// Vector error = x - x_exact; }
// REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
// }
} }
SECTION("Cholesky") SECTION("Cholesky")
...@@ -706,6 +696,177 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -706,6 +696,177 @@ TEST_CASE("LinearSolver", "[algebra]")
} }
#endif // PUGS_HAS_PETSC #endif // PUGS_HAS_PETSC
} }
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;
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("BICGStab 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("BICGStab 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("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.method() = LSMethod::gmres;
options.precond() = LSPrecond::none;
SECTION("GMRES 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("GMRES Diagonal")
{
options.precond() = LSPrecond::diagonal;
Vector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
Vector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
SECTION("GMRES ILU")
{
options.precond() = LSPrecond::incomplete_LU;
Vector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
Vector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
}
SECTION("LU")
{
LinearSolverOptions options;
options.library() = LSLibrary::petsc;
options.method() = LSMethod::lu;
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::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
}
} }
} }
} }
...@@ -769,6 +930,100 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -769,6 +930,100 @@ 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;
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("CG 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("CG ICholesky")
{
options.precond() = LSPrecond::incomplete_cholesky;
SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
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("Cholesky")
{
LinearSolverOptions options;
options.library() = LSLibrary::eigen3;
options.method() = LSMethod::cholesky;
options.precond() = LSPrecond::none;
SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
#else // PUGS_HAS_EIGEN3
SECTION("Eigen3 not linked")
{
LinearSolverOptions options;
options.library() = LSLibrary::eigen3;
options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;
REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
}
#endif // PUGS_HAS_EIGEN3
}
SECTION("PETSc") SECTION("PETSc")
{ {
#ifdef PUGS_HAS_PETSC #ifdef PUGS_HAS_PETSC
...@@ -870,6 +1125,8 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -870,6 +1125,8 @@ TEST_CASE("LinearSolver", "[algebra]")
} }
SECTION("none symmetric system") SECTION("none symmetric system")
{
SECTION("Dense matrix")
{ {
SmallMatrix<double> A{5}; SmallMatrix<double> A{5};
A = zero; A = zero;
...@@ -973,16 +1230,15 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -973,16 +1230,15 @@ TEST_CASE("LinearSolver", "[algebra]")
LinearSolver solver{options}; LinearSolver solver{options};
solver.solveLocalSystem(A, x, b); REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
SmallVector error = x - x_exact; "error: incomplete LU is not available for dense matrices in Eigen3");
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
} }
} }
SECTION("BICGStab2") SECTION("BICGStab2")
{ {
LinearSolverOptions options; LinearSolverOptions options;
options.library() = LSLibrary::petsc; options.library() = LSLibrary::eigen3;
options.method() = LSMethod::bicgstab2; options.method() = LSMethod::bicgstab2;
options.precond() = LSPrecond::none; options.precond() = LSPrecond::none;
...@@ -993,32 +1249,14 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -993,32 +1249,14 @@ TEST_CASE("LinearSolver", "[algebra]")
SmallVector<double> x{5}; SmallVector<double> x{5};
x = zero; x = zero;
LinearSolver solver{options}; REQUIRE_THROWS_WITH(LinearSolver{options}, "error: BICGStab2 is not an Eigen3 linear solver!");
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") SECTION("GMRES")
{ {
LinearSolverOptions options; LinearSolverOptions options;
options.library() = LSLibrary::petsc; options.library() = LSLibrary::eigen3;
options.method() = LSMethod::gmres; options.method() = LSMethod::gmres;
options.precond() = LSPrecond::none; options.precond() = LSPrecond::none;
...@@ -1059,16 +1297,15 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -1059,16 +1297,15 @@ TEST_CASE("LinearSolver", "[algebra]")
LinearSolver solver{options}; LinearSolver solver{options};
solver.solveLocalSystem(A, x, b); REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
SmallVector error = x - x_exact; "error: incomplete LU is not available for dense matrices in Eigen3");
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
} }
} }
SECTION("LU") SECTION("LU")
{ {
LinearSolverOptions options; LinearSolverOptions options;
options.library() = LSLibrary::petsc; options.library() = LSLibrary::eigen3;
options.method() = LSMethod::lu; options.method() = LSMethod::lu;
options.precond() = LSPrecond::none; options.precond() = LSPrecond::none;
...@@ -1086,7 +1323,7 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -1086,7 +1323,7 @@ TEST_CASE("LinearSolver", "[algebra]")
SECTION("Eigen3 not linked") SECTION("Eigen3 not linked")
{ {
LinearSolverOptions options; LinearSolverOptions options;
options.library() = LSLibrary::petsc; options.library() = LSLibrary::eigen3;
options.method() = LSMethod::cg; options.method() = LSMethod::cg;
options.precond() = LSPrecond::none; options.precond() = LSPrecond::none;
...@@ -1269,3 +1506,4 @@ TEST_CASE("LinearSolver", "[algebra]") ...@@ -1269,3 +1506,4 @@ TEST_CASE("LinearSolver", "[algebra]")
} }
} }
} }
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment