Add support of DenseMatrix and TinyMatrix in LinearSolver

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)

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)

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}

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);
 };

tests/test_LinearSolver.cpp

@@ -141,6 +141,8 @@ TEST_CASE("LinearSolver", "[algebra]")
 #endif   // PUGS_HAS_PETSC
   }
 
+  SECTION("Sparse Matrices")
+  {
     SECTION("check linear solvers")
     {
       SECTION("symmetric system")
@@ -526,3 +528,392 @@ TEST_CASE("LinearSolver", "[algebra]")
       }
     }
   }
+
+  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")
+          {
+            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")
+            {
+              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};
+            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;
+
+            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")
+          {
+            LinearSolverOptions options;
+            options.library() = LSLibrary::builtin;
+            options.method()  = LSMethod::bicgstab;
+            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("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;
+
+              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")
+            {
+              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("BICGStab2")
+          {
+            LinearSolverOptions options;
+            options.library() = LSLibrary::petsc;
+            options.method()  = LSMethod::bicgstab2;
+            options.precond() = LSPrecond::none;
+
+            SECTION("BICGStab2 no preconditioner")
+            {
+              options.precond() = LSPrecond::none;
+
+              Vector<double> x{5};
+              x = 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("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("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;
+
+            REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+          }
+#endif   // PUGS_HAS_PETSC
+        }
+      }
+    }
+  }
+}