diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp index dc62e370e6ec304cb2e4fd44d220218d9bd94283..82d3446791cd0a65065169d4f8d15f87120931a1 100644 --- a/src/algebra/BiCGStab.hpp +++ b/src/algebra/BiCGStab.hpp @@ -9,10 +9,10 @@ struct BiCGStab { - template <typename VectorType, typename MatrixType> + template <typename MatrixType, typename VectorType, typename RHSVectorType> BiCGStab(const MatrixType& A, VectorType& x, - const VectorType& b, + const RHSVectorType& b, const double epsilon, const size_t maximum_iteration, const bool verbose) diff --git a/src/algebra/CG.hpp b/src/algebra/CG.hpp index 7e653d8f95cddab4d5d1b7eb36481fd7690ef170..4ab52012a08732f21d57cf114bd873406bfd3898 100644 --- a/src/algebra/CG.hpp +++ b/src/algebra/CG.hpp @@ -8,10 +8,10 @@ struct CG { - template <typename VectorType, typename MatrixType> + template <typename MatrixType, typename VectorType, typename RHSVectorType> CG(const MatrixType& A, VectorType& x, - const VectorType& f, + const RHSVectorType& f, const double epsilon, const size_t maximum_iteration, const bool verbose = false) diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp index a7c6188c502b267cd54562b9cabb39b985b798d3..98e2a3a1c2caf8a041b6e69e71502bb6441647e7 100644 --- a/src/algebra/LinearSolver.cpp +++ b/src/algebra/LinearSolver.cpp @@ -131,10 +131,11 @@ struct LinearSolver::Internals } } + template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType> static void - builtinSolveLocalSystem(const CRSMatrix<double, size_t>& A, - Vector<double>& x, - const Vector<double>& b, + builtinSolveLocalSystem(const MatrixType& A, + SolutionVectorType& x, + const RHSVectorType& b, const LinearSolverOptions& options) { if (options.precond() != LSPrecond::none) { @@ -167,10 +168,11 @@ struct LinearSolver::Internals return 0; } + template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType> static void - petscSolveLocalSystem(const CRSMatrix<double, size_t>& A, - Vector<double>& x, - const Vector<double>& b, + petscSolveLocalSystem(const MatrixType& A, + SolutionVectorType& x, + const RHSVectorType& b, const LinearSolverOptions& options) { Assert(x.size() == b.size() and x.size() == A.numberOfColumns() and A.isSquare()); @@ -275,11 +277,9 @@ struct LinearSolver::Internals #else // PUGS_HAS_PETSC // LCOV_EXCL_START + template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType> static void - petscSolveLocalSystem(const CRSMatrix<double, size_t>&, - Vector<double>&, - const Vector<double>&, - const LinearSolverOptions&) + petscSolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&) { checkHasLibrary(LSLibrary::petsc); throw UnexpectedError("unexpected situation should not reach this point!"); @@ -287,6 +287,33 @@ 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::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 @@ -302,26 +329,15 @@ LinearSolver::checkOptions(const LinearSolverOptions& options) const } void -LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<double>& b) +LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<const double>& b) { - switch (m_options.library()) { - case LSLibrary::builtin: { - 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 - } + Internals::solveLocalSystem(A, x, b, m_options); +} + +void +LinearSolver::solveLocalSystem(const DenseMatrix<double>& A, Vector<double>& x, const Vector<const double>& b) +{ + Internals::solveLocalSystem(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 3ed5660bd9ce51c1e5ba5de3fe00c595aa96f1e7..3873374002a146c79a447a815d4ef6a57669883e 100644 --- a/src/algebra/LinearSolver.hpp +++ b/src/algebra/LinearSolver.hpp @@ -2,6 +2,7 @@ #define LINEAR_SOLVER_HPP #include <algebra/CRSMatrix.hpp> +#include <algebra/DenseMatrix.hpp> #include <algebra/LinearSolverOptions.hpp> #include <algebra/TinyMatrix.hpp> #include <algebra/TinyVector.hpp> @@ -15,21 +16,32 @@ class LinearSolver const LinearSolverOptions m_options; - void _solveLocalDense(size_t N, const double* A, double* x, const double* b); - 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); - template <size_t N> void solveLocalSystem(const TinyMatrix<N>& A, TinyVector<N>& x, const TinyVector<N>& b) { - this->_solveLocalDense(N, &A(0, 0), &x[0], &b[0]); + DenseMatrix A_dense{A}; + Vector<double> x_vector{N}; + Vector<double> b_vector{N}; + for (size_t i = 0; i < N; ++i) { + x_vector[i] = x[i]; + b_vector[i] = b[i]; + } + + this->solveLocalSystem(A_dense, x_vector, b_vector); + + for (size_t i = 0; i < N; ++i) { + x[i] = x_vector[i]; + } } + void solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<const double>& b); + void solveLocalSystem(const DenseMatrix<double>& A, Vector<double>& x, const Vector<const double>& b); + LinearSolver(); LinearSolver(const LinearSolverOptions& options); }; diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp index 5504ae390fd259eaa700a911b158bf6a6e46a625..1219c48e62d28a749383e115aef8b6e51a17e5ec 100644 --- a/tests/test_LinearSolver.cpp +++ b/tests/test_LinearSolver.cpp @@ -141,77 +141,52 @@ TEST_CASE("LinearSolver", "[algebra]") #endif // PUGS_HAS_PETSC } - SECTION("check linear solvers") + SECTION("Sparse Matrices") { - SECTION("symmetric system") + SECTION("check linear solvers") { - SparseMatrixDescriptor S{5}; - S(0, 0) = 2; - S(0, 1) = -1; + 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(1, 0) = -1; + S(1, 1) = 2; + S(1, 2) = -1; - S(2, 1) = -1; - S(2, 2) = 2; - S(2, 3) = -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(3, 2) = -1; + S(3, 3) = 2; + S(3, 4) = -1; - S(4, 3) = -1; - S(4, 4) = 2; + S(4, 3) = -1; + S(4, 4) = 2; - CRSMatrix A{S}; + 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<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; + Vector<double> b = A * x_exact; - SECTION("builtin") - { - SECTION("CG no preconditioner") + SECTION("builtin") { - 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(dot(error, error)) < 1E-10 * std::sqrt(dot(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") { + LinearSolverOptions options; + options.library() = LSLibrary::builtin; + options.method() = LSMethod::cg; options.precond() = LSPrecond::none; Vector<double> x{5}; @@ -223,24 +198,83 @@ TEST_CASE("LinearSolver", "[algebra]") Vector error = x - x_exact; REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact))); } + } - SECTION("CG Diagonal") + SECTION("PETSc") + { +#ifdef PUGS_HAS_PETSC + + SECTION("CG") { - options.precond() = LSPrecond::diagonal; + LinearSolverOptions options; + options.library() = LSLibrary::petsc; + options.method() = LSMethod::cg; + options.precond() = LSPrecond::none; + options.verbose() = true; - Vector<double> x{5}; - x = 0; + SECTION("CG no preconditioner") + { + options.precond() = LSPrecond::none; - LinearSolver solver{options}; + Vector<double> x{5}; + x = 0; - 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}; + + 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 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(dot(error, error)) < 1E-10 * std::sqrt(dot(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(dot(error, error)) < 1E-10 * std::sqrt(dot(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(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact))); + } } - SECTION("CG ICholeski") + SECTION("Choleski") { - options.precond() = LSPrecond::incomplete_choleski; + LinearSolverOptions options; + options.library() = LSLibrary::petsc; + options.method() = LSMethod::choleski; + options.precond() = LSPrecond::none; Vector<double> x{5}; x = 0; @@ -252,9 +286,63 @@ TEST_CASE("LinearSolver", "[algebra]") REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact))); } - SECTION("CG AMG") +#else // PUGS_HAS_PETSC + SECTION("PETSc not linked") { - options.precond() = LSPrecond::amg; + 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; @@ -267,121 +355,153 @@ TEST_CASE("LinearSolver", "[algebra]") } } - SECTION("Choleski") + SECTION("PETSc") { - LinearSolverOptions options; - options.library() = LSLibrary::petsc; - options.method() = LSMethod::choleski; - options.precond() = LSPrecond::none; +#ifdef PUGS_HAS_PETSC - Vector<double> x{5}; - x = 0; + SECTION("BICGStab") + { + LinearSolverOptions options; + options.library() = LSLibrary::petsc; + options.method() = LSMethod::bicgstab; + options.precond() = LSPrecond::none; + options.verbose() = true; - LinearSolver solver{options}; + SECTION("BICGStab no preconditioner") + { + options.precond() = LSPrecond::none; - 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 = 0; -#else // PUGS_HAS_PETSC - SECTION("PETSc not linked") - { - LinearSolverOptions options; - options.library() = LSLibrary::petsc; - options.method() = LSMethod::cg; - options.precond() = LSPrecond::none; + LinearSolver solver{options}; - REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!"); - } -#endif // PUGS_HAS_PETSC - } - } + 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("none symmetric system") - { - SparseMatrixDescriptor S{5}; - S(0, 0) = 2; - S(0, 1) = -1; + SECTION("BICGStab Diagonal") + { + options.precond() = LSPrecond::diagonal; + + Vector<double> x{5}; + x = 0; - S(1, 0) = -0.2; - S(1, 1) = 2; - S(1, 2) = -1; + LinearSolver solver{options}; - S(2, 1) = -1; - S(2, 2) = 4; - S(2, 3) = -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))); + } - S(3, 2) = -1; - S(3, 3) = 2; - S(3, 4) = -0.1; + SECTION("BICGStab ILU") + { + options.precond() = LSPrecond::incomplete_LU; - S(4, 3) = 1; - S(4, 4) = 3; + Vector<double> x{5}; + x = 0; - CRSMatrix A{S}; + LinearSolver solver{options}; - 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; - }(); + 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> b = A * x_exact; + SECTION("BICGStab2") + { + LinearSolverOptions options; + options.library() = LSLibrary::petsc; + options.method() = LSMethod::bicgstab2; + options.precond() = LSPrecond::none; - SECTION("builtin") - { - SECTION("BICGStab no preconditioner") - { - LinearSolverOptions options; - options.library() = LSLibrary::builtin; - options.method() = LSMethod::bicgstab; - options.precond() = LSPrecond::none; + SECTION("BICGStab2 no preconditioner") + { + options.precond() = LSPrecond::none; - Vector<double> x{5}; - x = 0; + Vector<double> x{5}; + x = 0; - LinearSolver solver{options}; + 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))); - } - } + 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("PETSc") - { -#ifdef PUGS_HAS_PETSC + SECTION("BICGStab2 Diagonal") + { + options.precond() = LSPrecond::diagonal; - SECTION("BICGStab") - { - LinearSolverOptions options; - options.library() = LSLibrary::petsc; - options.method() = LSMethod::bicgstab; - options.precond() = LSPrecond::none; - options.verbose() = true; + Vector<double> x{5}; + x = 0; - SECTION("BICGStab no preconditioner") + 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; - Vector<double> x{5}; - x = 0; + SECTION("GMRES no preconditioner") + { + options.precond() = LSPrecond::none; - LinearSolver solver{options}; + Vector<double> x{5}; + x = 0; - 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}; + + 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 = 0; + + 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 = 0; + + 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") + SECTION("LU") { - options.precond() = LSPrecond::diagonal; + LinearSolverOptions options; + options.library() = LSLibrary::petsc; + options.method() = LSMethod::lu; + options.precond() = LSPrecond::none; Vector<double> x{5}; x = 0; @@ -393,9 +513,69 @@ TEST_CASE("LinearSolver", "[algebra]") REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact))); } - SECTION("BICGStab ILU") +#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("Dense Matrices") + { + SECTION("check linear solvers") + { + SECTION("symmetric system") + { + DenseMatrix<double> A{5}; + A = zero; + + A(0, 0) = 2; + A(0, 1) = -1; + + A(1, 0) = -1; + A(1, 1) = 2; + A(1, 2) = -1; + + A(2, 1) = -1; + A(2, 2) = 2; + A(2, 3) = -1; + + A(3, 2) = -1; + A(3, 3) = 2; + A(3, 4) = -1; + + A(4, 3) = -1; + A(4, 4) = 2; + + 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") { - options.precond() = LSPrecond::incomplete_LU; + LinearSolverOptions options; + options.library() = LSLibrary::builtin; + options.method() = LSMethod::cg; + options.precond() = LSPrecond::none; Vector<double> x{5}; x = 0; @@ -408,15 +588,80 @@ TEST_CASE("LinearSolver", "[algebra]") } } - SECTION("BICGStab2") + SECTION("PETSc") { - LinearSolverOptions options; - options.library() = LSLibrary::petsc; - options.method() = LSMethod::bicgstab2; - options.precond() = LSPrecond::none; +#ifdef PUGS_HAS_PETSC - SECTION("BICGStab2 no preconditioner") + 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(dot(error, error)) < 1E-10 * std::sqrt(dot(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(dot(error, error)) < 1E-10 * std::sqrt(dot(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(dot(error, error)) < 1E-10 * std::sqrt(dot(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(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact))); + } + } + + SECTION("Choleski") + { + LinearSolverOptions options; + options.library() = LSLibrary::petsc; + options.method() = LSMethod::choleski; options.precond() = LSPrecond::none; Vector<double> x{5}; @@ -429,9 +674,63 @@ TEST_CASE("LinearSolver", "[algebra]") REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact))); } - SECTION("BICGStab2 Diagonal") +#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") + { + DenseMatrix<double> A{5}; + A = zero; + + 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; + + 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") { - options.precond() = LSPrecond::diagonal; + LinearSolverOptions options; + options.library() = LSLibrary::builtin; + options.method() = LSMethod::bicgstab; + options.precond() = LSPrecond::none; Vector<double> x{5}; x = 0; @@ -444,44 +743,153 @@ TEST_CASE("LinearSolver", "[algebra]") } } - 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; - Vector<double> x{5}; - x = 0; + SECTION("BICGStab no preconditioner") + { + options.precond() = LSPrecond::none; - LinearSolver solver{options}; + Vector<double> x{5}; + x = 0; - 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}; + + 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 = 0; + + 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 = 0; + + 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") + SECTION("BICGStab2") { - options.precond() = LSPrecond::diagonal; + LinearSolverOptions options; + options.library() = LSLibrary::petsc; + options.method() = LSMethod::bicgstab2; + options.precond() = LSPrecond::none; - Vector<double> x{5}; - x = 0; + SECTION("BICGStab2 no preconditioner") + { + options.precond() = LSPrecond::none; - LinearSolver solver{options}; + Vector<double> x{5}; + x = 0; - 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}; + + 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 = 0; + + 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 = 0; + + 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 = 0; + + 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 = 0; + + 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") + SECTION("LU") { - options.precond() = LSPrecond::incomplete_LU; + LinearSolverOptions options; + options.library() = LSLibrary::petsc; + options.method() = LSMethod::lu; + options.precond() = LSPrecond::none; Vector<double> x{5}; x = 0; @@ -492,36 +900,19 @@ TEST_CASE("LinearSolver", "[algebra]") 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 = 0; - - 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_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 + } } } }