Merge branch 'develop' into feature/reconstruction

tests/test_LinearSolver.cpp

@@ -836,6 +836,8 @@ TEST_CASE("LinearSolver", "[algebra]")
   }
 
   SECTION("Dense Matrices")
+  {
+    SECTION("SmallMatrix")
     {
       SECTION("check linear solvers")
       {
@@ -1470,4 +1472,592 @@ TEST_CASE("LinearSolver", "[algebra]")
         }
       }
     }
+
+    SECTION("TinyMatrix")
+    {
+      SECTION("check linear solvers")
+      {
+        SECTION("symmetric system")
+        {
+          TinyMatrix<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;
+
+          TinyVector<5> x_exact{1, 3, 2, 4, 5};
+
+          TinyVector b = A * x_exact;
+
+          SECTION("builtin")
+          {
+            SECTION("CG no preconditioner")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::builtin;
+              options.method()  = LSMethod::cg;
+              options.precond() = LSPrecond::none;
+
+              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("Eigen3")
+          {
+#ifdef PUGS_HAS_EIGEN3
+
+            SECTION("CG")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::cg;
+              options.precond() = LSPrecond::none;
+              options.verbose() = true;
+
+              SECTION("CG no preconditioner")
+              {
+                options.precond() = LSPrecond::none;
+
+                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("CG 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("CG ICholesky")
+              {
+                options.precond() = LSPrecond::incomplete_cholesky;
+
+                TinyVector<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;
+
+                REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not an Eigen3 preconditioner!");
+              }
+            }
+
+            SECTION("Cholesky")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::cholesky;
+              options.precond() = LSPrecond::none;
+
+              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)));
+            }
+
+#else    // PUGS_HAS_EIGEN3
+            SECTION("Eigen3 not linked")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::eigen3;
+              options.method()  = LSMethod::cg;
+              options.precond() = LSPrecond::none;
+
+              REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+            }
+#endif   // PUGS_HAS_EIGEN3
+          }
+
+          SECTION("PETSc")
+          {
+#ifdef PUGS_HAS_PETSC
+
+            SECTION("CG")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::cg;
+              options.precond() = LSPrecond::none;
+              options.verbose() = true;
+
+              SECTION("CG no preconditioner")
+              {
+                options.precond() = LSPrecond::none;
+
+                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("CG 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("CG ICholesky")
+              {
+                options.precond() = LSPrecond::incomplete_cholesky;
+
+                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("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("Cholesky")
+            {
+              LinearSolverOptions options;
+              options.library() = LSLibrary::petsc;
+              options.method()  = LSMethod::cholesky;
+              options.precond() = LSPrecond::none;
+
+              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)));
+            }
+
+#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")
+        {
+          SECTION("Dense matrix")
+          {
+            TinyMatrix<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;
+
+            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;
+
+                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("Eigen3")
+            {
+#ifdef PUGS_HAS_EIGEN3
+
+              SECTION("BICGStab")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::bicgstab;
+                options.precond() = LSPrecond::none;
+                options.verbose() = true;
+
+                SECTION("BICGStab no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
+
+                  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 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};
+
+                  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;
+
+                  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;
+
+                  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")
+                {
+                  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("LU")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::lu;
+                options.precond() = LSPrecond::none;
+
+                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)));
+              }
+
+#else    // PUGS_HAS_EIGEN3
+              SECTION("Eigen3 not linked")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::eigen3;
+                options.method()  = LSMethod::cg;
+                options.precond() = LSPrecond::none;
+
+                REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+              }
+#endif   // PUGS_HAS_EIGEN3
+            }
+
+            SECTION("PETSc")
+            {
+#ifdef PUGS_HAS_PETSC
+
+              SECTION("BICGStab")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::bicgstab;
+                options.precond() = LSPrecond::none;
+                options.verbose() = true;
+
+                SECTION("BICGStab no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
+
+                  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 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("BICGStab2")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::bicgstab2;
+                options.precond() = LSPrecond::none;
+
+                SECTION("BICGStab2 no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
+
+                  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("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")
+              {
+                LinearSolverOptions options;
+                options.library() = LSLibrary::petsc;
+                options.method()  = LSMethod::gmres;
+                options.precond() = LSPrecond::none;
+
+                SECTION("GMRES no preconditioner")
+                {
+                  options.precond() = LSPrecond::none;
+
+                  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")
+                {
+                  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;
+
+                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)));
+              }
+
+#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
+            }
+          }
+        }
+      }
+    }
+  }
 }