Begin linear solver interface

PETSc is almost plugged

CMakeLists.txt

@@ -2,6 +2,7 @@ cmake_minimum_required (VERSION 3.10)
 
 # CMake utils
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
+list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/packages/cmake-modules")
 
 # Forbids in-source builds
 include(CheckNotInSources)
@@ -113,13 +114,12 @@ if (NOT PUGS_ENABLE_MPI MATCHES "^(AUTO|ON|OFF)$")
   message(FATAL_ERROR "PUGS_ENABLE_MPI='${PUGS_ENABLE_MPI}'. Must be set to one of AUTO, ON or OFF")
 endif()
 
-# checks for MPI
+# Check for MPI
 if (PUGS_ENABLE_MPI MATCHES "^(AUTO|ON)$")
   set(MPI_DETERMINE_LIBRARY_VERSION TRUE)
   find_package(MPI)
 endif()
 
-
 #------------------------------------------------------
 # Search for ParMETIS
 
@@ -137,6 +137,16 @@ if(${MPI_FOUND})
   endif()
 endif()
 
+#------------------------------------------------------
+# Check for PETSc
+find_package(PETSc)
+set(PUGS_HAS_PETSC ${PETSC_FOUND})
+
+if (${PETSC_FOUND})
+  include_directories(SYSTEM ${PETSC_INCLUDES})
+endif()
+
+
 # -----------------------------------------------------
 
 if (${MPI_FOUND})
@@ -414,15 +424,18 @@ add_executable(
 target_link_libraries(
   pugs
   PugsMesh
+  PugsAlgebra
   PugsUtils
   PugsLanguage
   PugsLanguageAST
   PugsLanguageModules
   PugsLanguageAlgorithms
   PugsMesh
+  PugsAlgebra
   PugsUtils
   PugsLanguageUtils
   kokkos
+  ${PETSC_LIBRARIES}
   ${PARMETIS_LIBRARIES}
   ${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
   ${KOKKOS_CXX_FLAGS}
@@ -450,37 +463,43 @@ endif()
 
 message("")
 message("====== pugs build options ======")
-message(STATUS "version: ${PUGS_VERSION}")
-message(STATUS "build type: ${CMAKE_BUILD_TYPE}")
-message(STATUS "compiler: ${CMAKE_CXX_COMPILER_ID} ${CMAKE_CXX_COMPILER_VERSION}")
-message(STATUS "kokkos devices: ${PUGS_BUILD_KOKKOS_DEVICES}")
+message(" version: ${PUGS_VERSION}")
+message(" build type: ${CMAKE_BUILD_TYPE}")
+message(" compiler: ${CMAKE_CXX_COMPILER_ID} ${CMAKE_CXX_COMPILER_VERSION}")
+message(" kokkos devices: ${PUGS_BUILD_KOKKOS_DEVICES}")
 
 if (MPI_FOUND)
-  message(STATUS "MPI: ${MPI_CXX_LIBRARY_VERSION_STRING}")
+  message(" MPI: ${MPI_CXX_LIBRARY_VERSION_STRING}")
 else()
   if(NOT PARMETIS_LIBRARIES)
-    message(STATUS "MPI: deactivated: ParMETIS cannot be found!")
+    message(" MPI: deactivated: ParMETIS cannot be found!")
   else()
-    message(STATUS "MPI: not found!")
+    message(" MPI: not found!")
+  endif()
 endif()
+
+if (PETSC_FOUND)
+  message(" PETSc: ${PETSC_VERSION}")
+else()
+  message(" PETSc: not found!")
 endif()
 
 if(CLANG_FORMAT)
-  message(STATUS "clang-format: ${CLANG_FORMAT}")
+  message(" clang-format: ${CLANG_FORMAT}")
 else()
-  message(STATUS "clang-format: not found!")
+  message(" clang-format: not found!")
 endif()
 
 if(CLAZY_STANDALONE)
-  message(STATUS "clazy-standalone: ${CLAZY_STANDALONE}")
+  message(" clazy-standalone: ${CLAZY_STANDALONE}")
 else()
-  message(STATUS "clazy-standalone: no found!")
+  message(" clazy-standalone: no found!")
 endif()
 
 if (DOXYGEN_FOUND)
-  message(STATUS "doxygen: ${DOXYGEN_EXECUTABLE}")
+  message(" doxygen: ${DOXYGEN_EXECUTABLE}")
 else()
-  message(STATUS "doxygen: no found!")
+  message(" doxygen: no found!")
 endif()
 
 message("================================")

src/CMakeLists.txt

@@ -10,7 +10,7 @@ add_subdirectory(utils)
 add_subdirectory(language)
 
 # Pugs algebra
-#add_subdirectory(algebra)
+add_subdirectory(algebra)
 
 # Pugs mesh
 add_subdirectory(mesh)

src/algebra/CMakeLists.txt

+# ------------------- Source files --------------------
+
+add_library(
+  PugsAlgebra
+  LinearSystemSolver.cpp
+  PETScWrapper.cpp)

src/algebra/LinearSystemSolver.cpp

+#include <algebra/LinearSystemSolver.hpp>
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_PETSC
+#include <petsc.h>
+#endif   // PUGS_HAS_PETSC
+
+struct LinearSystemSolver::Internals
+{
+  static bool
+  hasLibrary(LSLibrary library)
+  {
+    switch (library) {
+    case LSLibrary::builtin: {
+      return true;
+    }
+    case LSLibrary::petsc: {
+#ifdef PUGS_HAS_PETSC
+      return true;
+#else
+      return false;
+#endif
+    }
+    default: {
+      throw UnexpectedError("Linear system library (" + name(library) + ") was not set!");
+    }
+    }
+  }
+
+  static void
+  checkHasLibrary(LSLibrary library)
+  {
+    if (not hasLibrary(library)) {
+      throw NormalError(name(library) + " is not linked to pugs. Cannot use it!");
+    }
+  }
+
+  static std::string
+  name(LSLibrary library)
+  {
+    switch (library) {
+    case LSLibrary::builtin: {
+      return "builtin";
+    }
+    case LSLibrary::petsc: {
+      return "PETSc";
+    }
+    default: {
+      throw UnexpectedError("Linear system library name is not defined!");
+    }
+    }
+  }
+
+  static std::string
+  name(LSMethod method)
+  {
+    switch (method) {
+    case LSMethod::cg: {
+      return "conjugate gradient";
+    }
+    case LSMethod::bicgstab: {
+      return "bicg-stab";
+    }
+    case LSMethod::gmres: {
+      return "gmres";
+    }
+    default: {
+      throw UnexpectedError("Linear system method name is not defined!");
+    }
+    }
+  }
+
+  static std::string
+  name(LSPrecond precond)
+  {
+    switch (precond) {
+    case LSPrecond::none: {
+      return "none";
+    }
+    case LSPrecond::diagonal: {
+      return "Jacobi (diagonal)";
+    }
+    case LSPrecond::incomplete_choleski: {
+      return "incomplete Choleski";
+    }
+    case LSPrecond::incomplete_LU: {
+      return "incomplete LU";
+    }
+    default: {
+      throw UnexpectedError("Linear system preconditioner name is not defined!");
+    }
+    }
+  }
+
+  static void
+  checkBuiltinMethod(LSMethod method)
+  {
+    switch (method) {
+    case LSMethod::cg:
+    case LSMethod::bicgstab: {
+      break;
+    }
+    default: {
+      throw NormalError(name(method) + " is not a builtin linear solver!");
+    }
+    }
+  }
+
+  static void
+  checkPETScMethod(LSMethod method)
+  {
+    switch (method) {
+    case LSMethod::cg:
+    case LSMethod::bicgstab:
+    case LSMethod::gmres: {
+      break;
+    }
+    default: {
+      throw NormalError(name(method) + " is not a builtin linear solver!");
+    }
+    }
+  }
+
+  static void
+  checkBuiltinPrecond(LSPrecond precond)
+  {
+    switch (precond) {
+    case LSPrecond::none: {
+      break;
+    }
+    default: {
+      throw NormalError(name(precond) + " is not a builtin preconditioner!");
+    }
+    }
+  }
+
+  static void
+  checkPETScPrecond(LSPrecond precond)
+  {
+    switch (precond) {
+    case LSPrecond::none:
+    case LSPrecond::diagonal:
+    case LSPrecond::incomplete_choleski:
+    case LSPrecond::incomplete_LU: {
+      break;
+    }
+    default: {
+      throw NormalError(name(precond) + " is not a PETSc preconditioner!");
+    }
+    }
+  }
+
+  static void
+  checkOptions(LSLibrary library, LSMethod method, LSPrecond precond)
+  {
+    switch (library) {
+    case LSLibrary::builtin: {
+      checkBuiltinMethod(method);
+      checkBuiltinPrecond(precond);
+      break;
+    }
+    case LSLibrary::petsc: {
+      checkPETScMethod(method);
+      checkPETScPrecond(precond);
+      break;
+    }
+    default: {
+      throw UnexpectedError("Undefined options compatibility for this library (" + name(library) + ")!");
+    }
+    }
+  }
+
+#ifdef PUGS_HAS_PETSC
+
+  static int
+  petscMonitor(KSP, int it, double norm, void*)
+  {
+    std::cout << "iteration: " << it << " residual norm = " << norm << '\n';
+    return 0;
+  }
+
+  static void
+  petscSolveLocalSystem(const SparseMatrixDescriptor<double, size_t>& A,
+                        Vector<double>& x,
+                        const Vector<double>& b,
+                        LSMethod method,
+                        LSPrecond precond,
+                        double epsilon,
+                        size_t maximum_iteration)
+  {
+    Assert(x.size() == b.size() and x.size() == A.numberOfRows());
+
+    Vec petscB;
+    VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, b.size(), b.size(), &b[0], &petscB);
+    Vec petscX;
+    VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, x.size(), x.size(), &x[0], &petscX);
+
+    Mat petscMat;
+    MatCreate(PETSC_COMM_WORLD, &petscMat);
+    MatSetSizes(petscMat, PETSC_DECIDE, PETSC_DECIDE, x.size(), x.size());
+
+    MatSetType(petscMat, MATAIJ);
+    MatSetUp(petscMat);
+
+    for (size_t i = 0; i < A.numberOfRows(); ++i) {
+      const auto& row = A.row(i);
+      PetscInt I      = i;
+      for (auto [j, value] : row.columnIdToValueMap()) {
+        PetscInt J = j;
+        MatSetValues(petscMat, 1, &I, 1, &J, &value, INSERT_VALUES);
+      }
+    }
+
+    MatAssemblyBegin(petscMat, MAT_FINAL_ASSEMBLY);
+    MatAssemblyEnd(petscMat, MAT_FINAL_ASSEMBLY);
+
+    PetscPrintf(PETSC_COMM_WORLD, "Matrix _before_ MatScaleUserImpl() operation\n");
+    MatView(petscMat, PETSC_VIEWER_STDOUT_SELF);
+
+    KSP ksp;
+
+    KSPCreate(PETSC_COMM_WORLD, &ksp);
+    KSPSetTolerances(ksp, epsilon, 1E-100, 1E5, maximum_iteration);
+
+    KSPSetOperators(ksp, petscMat, petscMat);
+
+    switch (method) {
+    case LSMethod::bicgstab: {
+      KSPSetType(ksp, KSPBCGS);
+      break;
+    }
+    case LSMethod::cg: {
+      KSPSetType(ksp, KSPCG);
+      break;
+    }
+    case LSMethod::gmres: {
+      KSPSetType(ksp, KSPGMRES);
+      break;
+    }
+    default: {
+      throw UnexpectedError("unexpected method: " + name(method));
+    }
+    }
+
+    PC pc;
+    KSPGetPC(ksp, &pc);
+
+    switch (precond) {
+    case LSPrecond::diagonal: {
+      PCSetType(pc, PCJACOBI);
+      break;
+    }
+    case LSPrecond::incomplete_LU: {
+      PCSetType(pc, PCILU);
+      break;
+    }
+    case LSPrecond::incomplete_choleski: {
+      PCSetType(pc, PCICC);
+      break;
+    }
+    case LSPrecond::none: {
+      PCSetType(pc, PCNONE);
+      break;
+    }
+    default: {
+      throw UnexpectedError("unexpected preconditioner: " + name(precond));
+    }
+    }
+
+    KSPMonitorSet(ksp, petscMonitor, 0, 0);
+
+    KSPSolve(ksp, petscB, petscX);
+
+    VecView(petscX, PETSC_VIEWER_STDOUT_SELF);
+
+    MatDestroy(&petscMat);
+  }
+
+#else
+
+  static void
+  petscSolveLocalSystem(const SparseMatrixDescriptor<double, size_t>&,
+                        Vector<double>&,
+                        const Vector<double>&,
+                        LSMethod,
+                        LSPrecond,
+                        double,
+                        size_t)
+  {
+    checkHasLibrary(LSLibrary::petsc);
+  }
+
+#endif   // PUGS_HAS_PETSC
+};
+
+bool
+LinearSystemSolver::hasLibrary(LSLibrary library) const
+{
+  return Internals::hasLibrary(library);
+}
+
+LinearSystemSolver::LinearSystemSolver(LSLibrary library,
+                                       LSMethod method,
+                                       LSPrecond precond,
+                                       double epsilon,
+                                       size_t maximum_iteration)
+  : m_library{library}, m_method{method}, m_precond{precond}, m_epsilon{epsilon}, m_maximum_iteration{maximum_iteration}
+{
+  Internals::checkHasLibrary(library);
+  Internals::checkOptions(library, method, precond);
+}
+
+void
+LinearSystemSolver::solveLocalSystem(const SparseMatrixDescriptor<double, size_t>& A,
+                                     Vector<double>& x,
+                                     const Vector<double>& b)
+{
+  switch (m_library) {
+  case LSLibrary::builtin: {
+    throw NotImplementedError("builtin solvers are not plugged yet!");
+    break;
+  }
+  case LSLibrary::petsc: {
+    Internals::petscSolveLocalSystem(A, x, b, m_method, m_precond, m_epsilon, m_maximum_iteration);
+    break;
+  }
+  default: {
+    throw UnexpectedError(Internals::name(m_library) + " cannot solve local systems for sparse matrices");
+  }
+  }
+}

src/algebra/LinearSystemSolver.hpp

+#ifndef LINEAR_SYSTEM_SOLVER_HPP
+#define LINEAR_SYSTEM_SOLVER_HPP
+
+#include <algebra/SparseMatrixDescriptor.hpp>
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+
+#include <utils/Exceptions.hpp>
+
+enum class LSLibrary
+{
+  builtin,
+  petsc
+};
+
+enum class LSMethod
+{
+  cg,
+  bicgstab,
+  gmres
+};
+
+enum class LSPrecond
+{
+  none,
+  diagonal,
+  incomplete_choleski,
+  incomplete_LU
+};
+
+class LinearSystemSolver
+{
+ private:
+  struct Internals;
+
+  const LSLibrary m_library;
+  const LSMethod m_method;
+  const LSPrecond m_precond;
+
+  const double m_epsilon;
+  const size_t m_maximum_iteration;
+
+  bool m_is_verbose = false;
+
+  void _solveLocalDense(size_t N, const double* A, double* x, const double* b);
+
+ public:
+  bool
+  isVerbose() const
+  {
+    return m_is_verbose;
+  }
+
+  void
+  setVerbose(bool verbose)
+  {
+    m_is_verbose = verbose;
+  }
+
+  bool hasLibrary(LSLibrary library) const;
+
+  void solveLocalSystem(const SparseMatrixDescriptor<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]);
+  }
+
+  LinearSystemSolver(LSLibrary library,
+                     LSMethod method,
+                     LSPrecond precond,
+                     double epsilon           = 1e-6,
+                     size_t maximum_iteration = 100);
+};
+
+#endif   // LINEAR_SYSTEM_SOLVER_HPP

src/algebra/PETScWrapper.cpp

+#include <algebra/PETScWrapper.hpp>
+
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_PETSC
+#include <petsc.h>
+#endif   // PUGS_HAS_PETSC
+
+namespace PETScWrapper
+{
+void
+initialize([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[])
+{
+#ifdef PUGS_HAS_PETSC
+  PetscInitialize(&argc, &argv, 0, 0);
+#endif   // PUGS_HAS_PETSC
+}
+
+void
+finalize()
+{
+#ifdef PUGS_HAS_PETSC
+  PetscFinalize();
+#endif   // PUGS_HAS_PETSC
+}
+}   // namespace PETScWrapper

src/algebra/PETScWrapper.hpp

+#ifndef PETSC_WRAPPER_HPP
+#define PETSC_WRAPPER_HPP
+
+namespace PETScWrapper
+{
+void initialize(int& argc, char* argv[]);
+void finalize();
+}   // namespace PETScWrapper
+
+#endif   // PETSC_WRAPPER_HPP

src/algebra/SparseMatrixDescriptor.hpp

@@ -22,6 +22,12 @@ class SparseMatrixDescriptor
     std::map<IndexType, DataType> m_id_value_map;
 
    public:
+    const auto&
+    columnIdToValueMap() const
+    {
+      return m_id_value_map;
+    }
+
     IndexType
     numberOfValues() const noexcept
     {

src/utils/BuildInfo.cpp

@@ -8,6 +8,10 @@
 #include <mpi.h>
 #endif   //  PUGS_HAS_MPI
 
+#ifdef PUGS_HAS_PETSC
+#include <petsc.h>
+#endif   // PUGS_HAS_PETSC
+
 std::string
 BuildInfo::type()
 {
@@ -42,3 +46,14 @@ BuildInfo::mpiLibrary()
   return "none";
 #endif   // PUGS_HAS_MPI
 }
+
+std::string
+BuildInfo::petscLibrary()
+{
+#ifdef PUGS_HAS_PETSC
+  return std::to_string(PETSC_VERSION_MAJOR) + "." + std::to_string(PETSC_VERSION_MINOR) + "." +
+         std::to_string(PETSC_VERSION_SUBMINOR);
+#else
+  return "none";
+#endif   // PUGS_HAS_PETSC
+}

src/utils/BuildInfo.hpp

@@ -9,6 +9,7 @@ struct BuildInfo
   static std::string compiler();
   static std::string kokkosDevices();
   static std::string mpiLibrary();
+  static std::string petscLibrary();
 };
 
 #endif   // BUILD_INFO_HPP

src/utils/PugsUtils.cpp

@@ -9,6 +9,8 @@
 #include <utils/FPEManager.hpp>
 #include <utils/SignalManager.hpp>
 
+#include <algebra/PETScWrapper.hpp>
+
 #include <rang.hpp>
 
 #include <Kokkos_Core.hpp>
@@ -41,7 +43,8 @@ initialize(int& argc, char* argv[])
   std::cout << "type:     " << rang::style::bold << BuildInfo::type() << rang::style::reset << '\n';
   std::cout << "compiler: " << rang::style::bold << BuildInfo::compiler() << rang::style::reset << '\n';
   std::cout << "kokkos:   " << rang::style::bold << BuildInfo::kokkosDevices() << rang::style::reset << '\n';
-  std::cout << "mpi:      " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
+  std::cout << "MPI:      " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
+  std::cout << "PETSc:    " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n';
   std::cout << "-------------------------------------------------------\n";
 
   std::string filename;
@@ -81,6 +84,8 @@ initialize(int& argc, char* argv[])
     SignalManager::init(enable_signals);
   }
 
+  PETScWrapper::initialize(argc, argv);
+
   Kokkos::initialize(argc, argv);
   std::cout << "-------------------- " << rang::fg::green << "exec info" << rang::fg::reset
             << " ------------------------" << '\n';
@@ -97,5 +102,6 @@ void
 finalize()
 {
   Kokkos::finalize();
+  PETScWrapper::finalize();
   parallel::Messenger::destroy();
 }

src/utils/pugs_config.hpp.in

@@ -3,6 +3,7 @@
 
 #cmakedefine PUGS_HAS_FENV_H
 #cmakedefine PUGS_HAS_MPI
+#cmakedefine PUGS_HAS_PETSC
 
 #cmakedefine SYSTEM_IS_LINUX
 #cmakedefine SYSTEM_IS_DARWIN

tests/CMakeLists.txt

@@ -93,6 +93,7 @@ target_link_libraries (unit_tests
   PugsLanguageAlgorithms
   PugsLanguageUtils
   PugsMesh
+  PugsAlgebra
   PugsUtils
   kokkos
   ${PARMETIS_LIBRARIES}
@@ -102,6 +103,7 @@ target_link_libraries (unit_tests
 
 target_link_libraries (mpi_unit_tests
   PugsUtils
+  PugsAlgebra
   PugsMesh
   kokkos
   ${PARMETIS_LIBRARIES}