From db22fdf658607f8a45c50995de1d5efcc51dabac Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Sun, 16 Mar 2025 19:43:11 +0100
Subject: [PATCH] Add missing tests

---
 tests/test_LinearSolver.cpp | 1254 +++++++++++++++++++++++++----------
 1 file changed, 922 insertions(+), 332 deletions(-)

diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 76ef6aa2b..286fe11ed 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -837,77 +837,52 @@ TEST_CASE("LinearSolver", "[algebra]")
 
   SECTION("Dense Matrices")
   {
-    SECTION("check linear solvers")
+    SECTION("SmallMatrix")
     {
-      SECTION("symmetric system")
+      SECTION("check linear solvers")
       {
-        SmallMatrix<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;
+        SECTION("symmetric system")
+        {
+          SmallMatrix<double> A{5};
+          A = zero;
 
-        SmallVector<const double> x_exact = [] {
-          SmallVector<double> y{5};
-          y[0] = 1;
-          y[1] = 3;
-          y[2] = 2;
-          y[3] = 4;
-          y[4] = 5;
-          return y;
-        }();
+          A(0, 0) = 2;
+          A(0, 1) = -1;
 
-        SmallVector<double> b = A * x_exact;
+          A(1, 0) = -1;
+          A(1, 1) = 2;
+          A(1, 2) = -1;
 
-        SECTION("builtin")
-        {
-          SECTION("CG no preconditioner")
-          {
-            LinearSolverOptions options;
-            options.library() = LSLibrary::builtin;
-            options.method()  = LSMethod::cg;
-            options.precond() = LSPrecond::none;
+          A(2, 1) = -1;
+          A(2, 2) = 2;
+          A(2, 3) = -1;
 
-            SmallVector<double> x{5};
-            x = zero;
+          A(3, 2) = -1;
+          A(3, 3) = 2;
+          A(3, 4) = -1;
 
-            LinearSolver solver{options};
+          A(4, 3) = -1;
+          A(4, 4) = 2;
 
-            solver.solveLocalSystem(A, x, b);
-            SmallVector error = x - x_exact;
-            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
-          }
-        }
+          SmallVector<const double> x_exact = [] {
+            SmallVector<double> y{5};
+            y[0] = 1;
+            y[1] = 3;
+            y[2] = 2;
+            y[3] = 4;
+            y[4] = 5;
+            return y;
+          }();
 
-        SECTION("Eigen3")
-        {
-#ifdef PUGS_HAS_EIGEN3
+          SmallVector<double> b = A * x_exact;
 
-          SECTION("CG")
+          SECTION("builtin")
           {
-            LinearSolverOptions options;
-            options.library() = LSLibrary::eigen3;
-            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;
 
               SmallVector<double> x{5};
@@ -919,89 +894,176 @@ TEST_CASE("LinearSolver", "[algebra]")
               SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
+          }
 
-            SECTION("CG Diagonal")
+          SECTION("Eigen3")
+          {
+#ifdef PUGS_HAS_EIGEN3
+
+            SECTION("CG")
             {
-              options.precond() = LSPrecond::diagonal;
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::cg;
+              options.precond() = LSPrecond::none;
+              options.verbose() = true;
 
-              SmallVector<double> x{5};
-              x = zero;
+              SECTION("CG no preconditioner")
+              {
+                options.precond() = LSPrecond::none;
 
-              LinearSolver solver{options};
+                SmallVector<double> x{5};
+                x = zero;
 
-              solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
-              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+              SECTION("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("CG ICholesky")
+            SECTION("Cholesky")
             {
-              options.precond() = LSPrecond::incomplete_cholesky;
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::cholesky;
+              options.precond() = LSPrecond::none;
 
               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");
+              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 AMG")
+#else    // PUGS_HAS_EIGEN3
+            SECTION("Eigen3 not linked")
             {
-              options.precond() = LSPrecond::amg;
-
-              SmallVector<double> x{5};
-              x = zero;
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::cg;
+              options.precond() = LSPrecond::none;
 
-              REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not an Eigen3 preconditioner!");
+              REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
             }
+#endif   // PUGS_HAS_EIGEN3
           }
 
-          SECTION("Cholesky")
+          SECTION("PETSc")
           {
-            LinearSolverOptions options;
-            options.library() = LSLibrary::eigen3;
-            options.method()  = LSMethod::cholesky;
-            options.precond() = LSPrecond::none;
+#ifdef PUGS_HAS_PETSC
 
-            SmallVector<double> x{5};
-            x = zero;
+            SECTION("CG")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::cg;
+              options.precond() = LSPrecond::none;
+              options.verbose() = true;
 
-            LinearSolver solver{options};
+              SECTION("CG no preconditioner")
+              {
+                options.precond() = LSPrecond::none;
 
-            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)));
-          }
+                SmallVector<double> x{5};
+                x = zero;
 
-#else    // PUGS_HAS_EIGEN3
-          SECTION("Eigen3 not linked")
-          {
-            LinearSolverOptions options;
-            options.library() = LSLibrary::eigen3;
-            options.method()  = LSMethod::cg;
-            options.precond() = LSPrecond::none;
+                LinearSolver solver{options};
 
-            REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
-          }
-#endif   // PUGS_HAS_EIGEN3
-        }
+                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("PETSc")
-        {
-#ifdef PUGS_HAS_PETSC
+              SECTION("CG Diagonal")
+              {
+                options.precond() = LSPrecond::diagonal;
 
-          SECTION("CG")
-          {
-            LinearSolverOptions options;
-            options.library() = LSLibrary::petsc;
-            options.method()  = LSMethod::cg;
-            options.precond() = LSPrecond::none;
-            options.verbose() = true;
+                SmallVector<double> x{5};
+                x = zero;
 
-            SECTION("CG no preconditioner")
+                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};
+
+                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 AMG")
+              {
+                options.precond() = LSPrecond::amg;
+
+                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("Cholesky")
             {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::cholesky;
               options.precond() = LSPrecond::none;
 
               SmallVector<double> x{5};
@@ -1014,133 +1076,448 @@ TEST_CASE("LinearSolver", "[algebra]")
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
-            SECTION("CG Diagonal")
+#else    // PUGS_HAS_PETSC
+            SECTION("PETSc not linked")
             {
-              options.precond() = LSPrecond::diagonal;
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::cg;
+              options.precond() = LSPrecond::none;
 
-              SmallVector<double> x{5};
-              x = zero;
+              REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+            }
+#endif   // PUGS_HAS_PETSC
+          }
+        }
 
-              LinearSolver solver{options};
+        SECTION("none symmetric system")
+        {
+          SECTION("Dense matrix")
+          {
+            SmallMatrix<double> A{5};
+            A = zero;
 
-              solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
-              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
-            }
+            A(0, 0) = 2;
+            A(0, 1) = -1;
 
-            SECTION("CG ICholesky")
+            A(1, 0) = -0.2;
+            A(1, 1) = 2;
+            A(1, 2) = -1;
+
+            A(2, 1) = -1;
+            A(2, 2) = 4;
+            A(2, 3) = -2;
+
+            A(3, 2) = -1;
+            A(3, 3) = 2;
+            A(3, 4) = -0.1;
+
+            A(4, 3) = 1;
+            A(4, 4) = 3;
+
+            SmallVector<const double> x_exact = [] {
+              SmallVector<double> y{5};
+              y[0] = 1;
+              y[1] = 3;
+              y[2] = 2;
+              y[3] = 4;
+              y[4] = 5;
+              return y;
+            }();
+
+            SmallVector<double> b = A * x_exact;
+
+            SECTION("builtin")
             {
-              options.precond() = LSPrecond::incomplete_cholesky;
+              SECTION("BICGStab no preconditioner")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::builtin;
+                options.method()  = LSMethod::bicgstab;
+                options.precond() = LSPrecond::none;
 
-              SmallVector<double> x{5};
-              x = zero;
+                SmallVector<double> x{5};
+                x = zero;
 
-              LinearSolver solver{options};
+                LinearSolver solver{options};
 
-              solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
-              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
             }
 
-            SECTION("CG AMG")
+            SECTION("Eigen3")
             {
-              options.precond() = LSPrecond::amg;
+#ifdef PUGS_HAS_EIGEN3
 
-              SmallVector<double> x{5};
-              x = zero;
+              SECTION("BICGStab")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                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);
-              SmallVector error = x - x_exact;
-              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  SmallVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
+
+                SECTION("BICGStab Diagonal")
+                {
+                  options.precond() = LSPrecond::diagonal;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  SmallVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
+
+                SECTION("BICGStab ILU")
+                {
+                  options.precond() = LSPrecond::incomplete_LU;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+                                      "error: incomplete LU is not available for dense matrices in Eigen3");
+                }
+              }
+
+              SECTION("BICGStab2")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::bicgstab2;
+                options.precond() = LSPrecond::none;
+
+                SECTION("BICGStab2 no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  REQUIRE_THROWS_WITH(LinearSolver{options}, "error: BICGStab2 is not an Eigen3 linear solver!");
+                }
+              }
+
+              SECTION("GMRES")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::gmres;
+                options.precond() = LSPrecond::none;
+
+                SECTION("GMRES no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  SmallVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
+
+                SECTION("GMRES Diagonal")
+                {
+                  options.precond() = LSPrecond::diagonal;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  SmallVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
+
+                SECTION("GMRES ILU")
+                {
+                  options.precond() = LSPrecond::incomplete_LU;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+                                      "error: incomplete LU is not available for dense matrices in Eigen3");
+                }
+              }
+
+              SECTION("LU")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::lu;
+                options.precond() = LSPrecond::none;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
+
+#else    // PUGS_HAS_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("Cholesky")
-          {
-            LinearSolverOptions options;
-            options.library() = LSLibrary::petsc;
-            options.method()  = LSMethod::cholesky;
-            options.precond() = LSPrecond::none;
+            SECTION("PETSc")
+            {
+#ifdef PUGS_HAS_PETSC
 
-            SmallVector<double> x{5};
-            x = zero;
+              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);
-            SmallVector error = x - x_exact;
-            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
-          }
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  SmallVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
+
+                SECTION("BICGStab Diagonal")
+                {
+                  options.precond() = LSPrecond::diagonal;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  SmallVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
+
+                SECTION("BICGStab ILU")
+                {
+                  options.precond() = LSPrecond::incomplete_LU;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  SmallVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
+              }
+
+              SECTION("BICGStab2")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::bicgstab2;
+                options.precond() = LSPrecond::none;
+
+                SECTION("BICGStab2 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("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")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::gmres;
+                options.precond() = LSPrecond::none;
+
+                SECTION("GMRES no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  SmallVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
+
+                SECTION("GMRES Diagonal")
+                {
+                  options.precond() = LSPrecond::diagonal;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  SmallVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
+
+                SECTION("GMRES ILU")
+                {
+                  options.precond() = LSPrecond::incomplete_LU;
+
+                  SmallVector<double> x{5};
+                  x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  SmallVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
+              }
+
+              SECTION("LU")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::lu;
+                options.precond() = LSPrecond::none;
+
+                SmallVector<double> x{5};
+                x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                SmallVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+              }
 
 #else    // PUGS_HAS_PETSC
-          SECTION("PETSc not linked")
-          {
-            LinearSolverOptions options;
-            options.library() = LSLibrary::petsc;
-            options.method()  = LSMethod::cg;
-            options.precond() = LSPrecond::none;
+              SECTION("PETSc not linked")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::cg;
+                options.precond() = LSPrecond::none;
 
-            REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
-          }
+                REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+              }
 #endif   // PUGS_HAS_PETSC
+            }
+          }
         }
       }
+    }
 
-      SECTION("none symmetric system")
+    SECTION("TinyMatrix")
+    {
+      SECTION("check linear solvers")
       {
-        SECTION("Dense matrix")
+        SECTION("symmetric system")
         {
-          SmallMatrix<double> A{5};
-          A = zero;
+          TinyMatrix<5> A = zero;
 
           A(0, 0) = 2;
           A(0, 1) = -1;
 
-          A(1, 0) = -0.2;
+          A(1, 0) = -1;
           A(1, 1) = 2;
           A(1, 2) = -1;
 
           A(2, 1) = -1;
-          A(2, 2) = 4;
-          A(2, 3) = -2;
+          A(2, 2) = 2;
+          A(2, 3) = -1;
 
           A(3, 2) = -1;
           A(3, 3) = 2;
-          A(3, 4) = -0.1;
+          A(3, 4) = -1;
 
-          A(4, 3) = 1;
-          A(4, 4) = 3;
+          A(4, 3) = -1;
+          A(4, 4) = 2;
 
-          SmallVector<const double> x_exact = [] {
-            SmallVector<double> y{5};
-            y[0] = 1;
-            y[1] = 3;
-            y[2] = 2;
-            y[3] = 4;
-            y[4] = 5;
-            return y;
-          }();
+          TinyVector<5> x_exact{1, 3, 2, 4, 5};
 
-          SmallVector<double> b = A * x_exact;
+          TinyVector b = A * x_exact;
 
           SECTION("builtin")
           {
-            SECTION("BICGStab no preconditioner")
+            SECTION("CG no preconditioner")
             {
               LinearSolverOptions options;
               options.library() = LSLibrary::builtin;
-              options.method()  = LSMethod::bicgstab;
+              options.method()  = LSMethod::cg;
               options.precond() = LSPrecond::none;
 
-              SmallVector<double> x{5};
-              x = zero;
+              TinyVector<5> x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
+              TinyVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
           }
@@ -1149,322 +1526,535 @@ TEST_CASE("LinearSolver", "[algebra]")
           {
 #ifdef PUGS_HAS_EIGEN3
 
-            SECTION("BICGStab")
+            SECTION("CG")
             {
               LinearSolverOptions options;
               options.library() = LSLibrary::eigen3;
-              options.method()  = LSMethod::bicgstab;
+              options.method()  = LSMethod::cg;
               options.precond() = LSPrecond::none;
               options.verbose() = true;
 
-              SECTION("BICGStab no preconditioner")
+              SECTION("CG no preconditioner")
               {
                 options.precond() = LSPrecond::none;
 
-                SmallVector<double> x{5};
-                x = zero;
+                TinyVector<5> x = zero;
 
                 LinearSolver solver{options};
 
                 solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
+                TinyVector error = x - x_exact;
                 REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
               }
 
-              SECTION("BICGStab Diagonal")
+              SECTION("CG Diagonal")
               {
                 options.precond() = LSPrecond::diagonal;
 
-                SmallVector<double> x{5};
-                x = zero;
+                TinyVector<5> x = zero;
 
                 LinearSolver solver{options};
 
                 solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
+                TinyVector error = x - x_exact;
                 REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
               }
 
-              SECTION("BICGStab ILU")
+              SECTION("CG ICholesky")
               {
-                options.precond() = LSPrecond::incomplete_LU;
+                options.precond() = LSPrecond::incomplete_cholesky;
 
-                SmallVector<double> x{5};
-                x = zero;
+                TinyVector<5> x = zero;
 
                 LinearSolver solver{options};
 
                 REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
-                                    "error: incomplete LU is not available for dense matrices in Eigen3");
+                                    "error: incomplete cholesky is not available for dense matrices in Eigen3");
+              }
+
+              SECTION("CG AMG")
+              {
+                options.precond() = LSPrecond::amg;
+
+                REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not an Eigen3 preconditioner!");
               }
             }
 
-            SECTION("BICGStab2")
+            SECTION("Cholesky")
             {
               LinearSolverOptions options;
               options.library() = LSLibrary::eigen3;
-              options.method()  = LSMethod::bicgstab2;
+              options.method()  = LSMethod::cholesky;
               options.precond() = LSPrecond::none;
 
-              SECTION("BICGStab2 no preconditioner")
-              {
-                options.precond() = LSPrecond::none;
+              TinyVector<5> x = zero;
 
-                SmallVector<double> x{5};
-                x = zero;
+              LinearSolver solver{options};
 
-                REQUIRE_THROWS_WITH(LinearSolver{options}, "error: BICGStab2 is not an Eigen3 linear solver!");
-              }
+              solver.solveLocalSystem(A, x, b);
+              TinyVector error = x - x_exact;
+              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
-            SECTION("GMRES")
+#else    // PUGS_HAS_EIGEN3
+            SECTION("Eigen3 not linked")
             {
               LinearSolverOptions options;
               options.library() = LSLibrary::eigen3;
-              options.method()  = LSMethod::gmres;
+              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")
+          {
+#ifdef PUGS_HAS_PETSC
+
+            SECTION("CG")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::cg;
               options.precond() = LSPrecond::none;
+              options.verbose() = true;
 
-              SECTION("GMRES no preconditioner")
+              SECTION("CG no preconditioner")
               {
                 options.precond() = LSPrecond::none;
 
-                SmallVector<double> x{5};
-                x = zero;
+                TinyVector<5> x = zero;
 
                 LinearSolver solver{options};
 
                 solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
+                TinyVector error = x - x_exact;
                 REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
               }
 
-              SECTION("GMRES Diagonal")
+              SECTION("CG Diagonal")
               {
                 options.precond() = LSPrecond::diagonal;
 
-                SmallVector<double> x{5};
-                x = zero;
+                TinyVector<5> x = zero;
 
                 LinearSolver solver{options};
 
                 solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
+                TinyVector error = x - x_exact;
                 REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
               }
 
-              SECTION("GMRES ILU")
+              SECTION("CG ICholesky")
               {
-                options.precond() = LSPrecond::incomplete_LU;
+                options.precond() = LSPrecond::incomplete_cholesky;
 
-                SmallVector<double> x{5};
-                x = zero;
+                TinyVector<5> x = zero;
 
                 LinearSolver solver{options};
 
-                REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
-                                    "error: incomplete LU is not available for dense matrices in Eigen3");
+                solver.solveLocalSystem(A, x, b);
+                TinyVector 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;
+
+                TinyVector<5> x = zero;
+
+                LinearSolver solver{options};
+
+                solver.solveLocalSystem(A, x, b);
+                TinyVector error = x - x_exact;
+                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
               }
             }
 
-            SECTION("LU")
+            SECTION("Cholesky")
             {
               LinearSolverOptions options;
-              options.library() = LSLibrary::eigen3;
-              options.method()  = LSMethod::lu;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::cholesky;
               options.precond() = LSPrecond::none;
 
-              SmallVector<double> x{5};
-              x = zero;
+              TinyVector<5> x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
+              TinyVector 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")
+#else    // PUGS_HAS_PETSC
+            SECTION("PETSc not linked")
             {
               LinearSolverOptions options;
-              options.library() = LSLibrary::eigen3;
+              options.library() = LSLibrary::petsc;
               options.method()  = LSMethod::cg;
               options.precond() = LSPrecond::none;
 
-              REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+              REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
             }
-#endif   // PUGS_HAS_EIGEN3
+#endif   // PUGS_HAS_PETSC
           }
+        }
 
-          SECTION("PETSc")
+        SECTION("none symmetric system")
+        {
+          SECTION("Dense matrix")
           {
-#ifdef PUGS_HAS_PETSC
+            TinyMatrix<5> A = zero;
 
-            SECTION("BICGStab")
-            {
-              LinearSolverOptions options;
-              options.library() = LSLibrary::petsc;
-              options.method()  = LSMethod::bicgstab;
-              options.precond() = LSPrecond::none;
-              options.verbose() = true;
+            A(0, 0) = 2;
+            A(0, 1) = -1;
+
+            A(1, 0) = -0.2;
+            A(1, 1) = 2;
+            A(1, 2) = -1;
+
+            A(2, 1) = -1;
+            A(2, 2) = 4;
+            A(2, 3) = -2;
+
+            A(3, 2) = -1;
+            A(3, 3) = 2;
+            A(3, 4) = -0.1;
+
+            A(4, 3) = 1;
+            A(4, 4) = 3;
+
+            const TinyVector<5> x_exact{1, 3, 2, 4, 5
+
+            };
+
+            const TinyVector b = A * x_exact;
 
+            SECTION("builtin")
+            {
               SECTION("BICGStab no preconditioner")
               {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::builtin;
+                options.method()  = LSMethod::bicgstab;
                 options.precond() = LSPrecond::none;
 
-                SmallVector<double> x{5};
-                x = zero;
+                TinyVector<5> x = zero;
 
                 LinearSolver solver{options};
 
                 solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
+                TinyVector error = x - x_exact;
                 REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
               }
+            }
+
+            SECTION("Eigen3")
+            {
+#ifdef PUGS_HAS_EIGEN3
 
-              SECTION("BICGStab Diagonal")
+              SECTION("BICGStab")
               {
-                options.precond() = LSPrecond::diagonal;
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::bicgstab;
+                options.precond() = LSPrecond::none;
+                options.verbose() = true;
 
-                SmallVector<double> x{5};
-                x = zero;
+                SECTION("BICGStab no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
 
-                LinearSolver solver{options};
+                  TinyVector<5> x = zero;
 
-                solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
-                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
-              }
+                  LinearSolver solver{options};
 
-              SECTION("BICGStab ILU")
-              {
-                options.precond() = LSPrecond::incomplete_LU;
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
 
-                SmallVector<double> x{5};
-                x = zero;
+                SECTION("BICGStab Diagonal")
+                {
+                  options.precond() = LSPrecond::diagonal;
 
-                LinearSolver solver{options};
+                  TinyVector<5> x = zero;
 
-                solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
-                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector 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;
+
+                  TinyVector<5> x = zero;
+
+                  LinearSolver solver{options};
+
+                  REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+                                      "error: incomplete LU is not available for dense matrices in Eigen3");
+                }
               }
-            }
 
-            SECTION("BICGStab2")
-            {
-              LinearSolverOptions options;
-              options.library() = LSLibrary::petsc;
-              options.method()  = LSMethod::bicgstab2;
-              options.precond() = LSPrecond::none;
+              SECTION("BICGStab2")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::bicgstab2;
+                options.precond() = LSPrecond::none;
 
-              SECTION("BICGStab2 no preconditioner")
+                SECTION("BICGStab2 no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
+
+                  REQUIRE_THROWS_WITH(LinearSolver{options}, "error: BICGStab2 is not an Eigen3 linear solver!");
+                }
+              }
+
+              SECTION("GMRES")
               {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::gmres;
                 options.precond() = LSPrecond::none;
 
-                SmallVector<double> x{5};
-                x = zero;
+                SECTION("GMRES no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
 
-                LinearSolver solver{options};
+                  TinyVector<5> x = zero;
 
-                solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
-                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector 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;
+
+                  TinyVector<5> x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector 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;
+
+                  TinyVector<5> x = zero;
+
+                  LinearSolver solver{options};
+
+                  REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+                                      "error: incomplete LU is not available for dense matrices in Eigen3");
+                }
               }
 
-              SECTION("BICGStab2 Diagonal")
+              SECTION("LU")
               {
-                options.precond() = LSPrecond::diagonal;
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::lu;
+                options.precond() = LSPrecond::none;
 
-                SmallVector<double> x{5};
-                x = zero;
+                TinyVector<5> x = zero;
 
                 LinearSolver solver{options};
 
                 solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
+                TinyVector error = x - x_exact;
                 REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
               }
+
+#else    // PUGS_HAS_EIGEN3
+              SECTION("Eigen3 not linked")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::cg;
+                options.precond() = LSPrecond::none;
+
+                REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+              }
+#endif   // PUGS_HAS_EIGEN3
             }
 
-            SECTION("GMRES")
+            SECTION("PETSc")
             {
-              LinearSolverOptions options;
-              options.library() = LSLibrary::petsc;
-              options.method()  = LSMethod::gmres;
-              options.precond() = LSPrecond::none;
+#ifdef PUGS_HAS_PETSC
 
-              SECTION("GMRES no preconditioner")
+              SECTION("BICGStab")
               {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::bicgstab;
                 options.precond() = LSPrecond::none;
+                options.verbose() = true;
 
-                SmallVector<double> x{5};
-                x = zero;
+                SECTION("BICGStab no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
 
-                LinearSolver solver{options};
+                  TinyVector<5> x = zero;
 
-                solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
-                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector 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;
+
+                  TinyVector<5> x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector 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;
+
+                  TinyVector<5> x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
               }
 
-              SECTION("GMRES Diagonal")
+              SECTION("BICGStab2")
               {
-                options.precond() = LSPrecond::diagonal;
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::bicgstab2;
+                options.precond() = LSPrecond::none;
 
-                SmallVector<double> x{5};
-                x = zero;
+                SECTION("BICGStab2 no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
 
-                LinearSolver solver{options};
+                  TinyVector<5> x = zero;
 
-                solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
-                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector 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;
+
+                  TinyVector<5> x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector error = x - x_exact;
+                  REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                }
               }
 
-              SECTION("GMRES ILU")
+              SECTION("GMRES")
               {
-                options.precond() = LSPrecond::incomplete_LU;
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::gmres;
+                options.precond() = LSPrecond::none;
 
-                SmallVector<double> x{5};
-                x = zero;
+                SECTION("GMRES no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
 
-                LinearSolver solver{options};
+                  TinyVector<5> x = zero;
 
-                solver.solveLocalSystem(A, x, b);
-                SmallVector error = x - x_exact;
-                REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector 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;
+
+                  TinyVector<5> x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector 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;
+
+                  TinyVector<5> x = zero;
+
+                  LinearSolver solver{options};
+
+                  solver.solveLocalSystem(A, x, b);
+                  TinyVector 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;
+              SECTION("LU")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::lu;
+                options.precond() = LSPrecond::none;
 
-              SmallVector<double> x{5};
-              x = zero;
+                TinyVector<5> x = zero;
 
-              LinearSolver solver{options};
+                LinearSolver solver{options};
 
-              solver.solveLocalSystem(A, x, b);
-              SmallVector error = x - x_exact;
-              REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
-            }
+                solver.solveLocalSystem(A, x, b);
+                TinyVector 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