From ef90e8c14cce8467fb0b9a0cfdac0fabff5049fe Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 24 Jun 2021 20:16:54 +0200
Subject: [PATCH] Add support of DenseMatrix and TinyMatrix in LinearSolver

---
 src/algebra/BiCGStab.hpp     |   4 +-
 src/algebra/CG.hpp           |   4 +-
 src/algebra/LinearSolver.cpp |  74 ++--
 src/algebra/LinearSolver.hpp |  22 +-
 tests/test_LinearSolver.cpp  | 819 ++++++++++++++++++++++++++---------
 5 files changed, 671 insertions(+), 252 deletions(-)

diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
index dc62e370e..82d344679 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 7e653d8f9..4ab52012a 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 a7c6188c5..98e2a3a1c 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 3ed5660bd..387337400 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 5504ae390..1219c48e6 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
+        }
       }
     }
   }
-- 
GitLab