Implement the HyKKT Cholesky module

cmake/FindSuiteSparse.cmake

@@ -13,6 +13,7 @@ set(SUITESPARSE_MODULES
   amd
   colamd
   klu
+  cholmod
   suitesparseconfig)
 
 find_library(SUITESPARSE_LIBRARY

resolve/CMakeLists.txt

@@ -164,7 +164,6 @@ endif()
 if(RESOLVE_USE_KLU)
   add_subdirectory(hykkt)
 endif()
-#list(APPEND ReSolve_Targets_List resolve_hykkt)
 
 # Set installable targets
 install(TARGETS ${ReSolve_Targets_List} EXPORT ReSolveTargets)

resolve/MemoryUtils.tpp

@@ -1,82 +1,81 @@
 /**
  * @file MemoryUtils.tpp
- * 
+ *
  * Contains implementation of memory utility functions wrappers.
  * All it does it calls vendor specific functions frm an abstract interface.
- * 
+ *
  * @author Slaven Peles <peless@ornl.gov>
  */
 
 #pragma once
 
-
 namespace ReSolve
 {
-    template <class Policy>
-    void MemoryUtils<Policy>::deviceSynchronize()
-    {
-      Policy::deviceSynchronize();
-    }
+  template <class Policy>
+  void MemoryUtils<Policy>::deviceSynchronize()
+  {
+    Policy::deviceSynchronize();
+  }
+
+  template <class Policy>
+  int MemoryUtils<Policy>::getLastDeviceError()
+  {
+    return Policy::getLastDeviceError();
+  }
+
+  template <class Policy>
+  int MemoryUtils<Policy>::deleteOnDevice(void* v)
+  {
+    return Policy::deleteOnDevice(v);
+  }
+
+  template <class Policy>
+  template <typename I, typename T>
+  int MemoryUtils<Policy>::allocateArrayOnDevice(T** v, I n)
+  {
+    return Policy::template allocateArrayOnDevice<I, T>(v, n);
+  }
+
+  template <class Policy>
+  template <typename I, typename T>
+  int MemoryUtils<Policy>::allocateBufferOnDevice(T** v, I n)
+  {
+    return Policy::template allocateBufferOnDevice<I, T>(v, n);
+  }
+
+  template <class Policy>
+  template <typename I, typename T>
+  int MemoryUtils<Policy>::setZeroArrayOnDevice(T* v, I n)
+  {
+    return Policy::template setZeroArrayOnDevice<I, T>(v, n);
+  }
+
+  template <class Policy>
+  template <typename I, typename T>
+  int MemoryUtils<Policy>::setArrayToConstOnDevice(T* v, T c, I n)
+  {
+    return Policy::template setArrayToConstOnDevice<I, T>(v, c, n);
+  }
+
+  template <class Policy>
+  template <typename I, typename T>
+  int MemoryUtils<Policy>::copyArrayDeviceToHost(T* dst, const T* src, I n)
+  {
+    return Policy::template copyArrayDeviceToHost<I, T>(dst, src, n);
+  }
 
-    template <class Policy>
-    int MemoryUtils<Policy>::getLastDeviceError()
-    {
-      return Policy::getLastDeviceError();
-    }
+  template <class Policy>
+  template <typename I, typename T>
+  int MemoryUtils<Policy>::copyArrayDeviceToDevice(T* dst, const T* src, I n)
+  {
+    return Policy::template copyArrayDeviceToDevice<I, T>(dst, src, n);
+  }
 
-    template <class Policy>
-    int MemoryUtils<Policy>::deleteOnDevice(void* v)
-    {
-      return Policy::deleteOnDevice(v);
-    }
-
-    template <class Policy>
-    template <typename I, typename T>
-    int MemoryUtils<Policy>::allocateArrayOnDevice(T** v, I n)
-    {
-      return Policy::template allocateArrayOnDevice<I, T>(v, n);
-    }
-
-    template <class Policy>
-    template <typename I, typename T>
-    int MemoryUtils<Policy>::allocateBufferOnDevice(T** v, I n)
-    {
-      return Policy::template allocateBufferOnDevice<I, T>(v, n);
-    }
-
-    template <class Policy>
-    template <typename I, typename T>
-    int MemoryUtils<Policy>::setZeroArrayOnDevice(T* v, I n)
-    {
-      return Policy::template setZeroArrayOnDevice<I, T>(v, n);
-    }
-
-    template <class Policy>
-    template <typename I, typename T>
-    int MemoryUtils<Policy>::setArrayToConstOnDevice(T* v, T c, I n)
-    {
-      return Policy::template setArrayToConstOnDevice<I, T>(v, c, n);
-    }
-
-    template <class Policy>
-    template <typename I, typename T>
-    int MemoryUtils<Policy>::copyArrayDeviceToHost(T* dst, const T* src, I n)
-    {
-      return Policy::template copyArrayDeviceToHost<I, T>(dst, src, n);
-    }
-
-    template <class Policy>
-    template <typename I, typename T>
-    int MemoryUtils<Policy>::copyArrayDeviceToDevice(T* dst, const T* src, I n)
-    {
-      return Policy::template copyArrayDeviceToDevice<I, T>(dst, src, n);
-    }
-
-    template <class Policy>
-    template <typename I, typename T>
-    int MemoryUtils<Policy>::copyArrayHostToDevice(T* dst, const T* src, I n)
-    {
-      return Policy::template copyArrayHostToDevice<I, T>(dst, src, n);
-    }
+  template <class Policy>
+  template <typename I, typename T>
+  int MemoryUtils<Policy>::copyArrayHostToDevice(T* dst, const T* src, I n)
+  {
+    return Policy::template copyArrayHostToDevice<I, T>(dst, src, n);
+  }
 
 } // namespace ReSolve

resolve/hykkt/CMakeLists.txt

@@ -1,2 +1,3 @@
 add_subdirectory(permutation)
 add_subdirectory(ruiz)
+add_subdirectory(cholesky)

resolve/hykkt/cholesky/CMakeLists.txt

@@ -0,0 +1,65 @@
+#[[
+
+@brief Build ReSolve matrix module
+
+@author Slaven Peles <peless@ornl.gov>
+
+]]
+
+# C++ code
+set(HyKKT_CHOL_SRC
+    CholeskySolver.cpp
+    CholeskySolverCpu.cpp
+)
+
+# C++ code that depends on CUDA SDK libraries
+set(HyKKT_CUDASDK_SRC
+    CholeskySolverCuda.cpp
+)
+
+# and on HIP
+set(HyKKT_ROCM_SRC
+    CholeskySolverHip.cpp
+)
+
+# Header files to be installed
+set(HyKKT_CHOL_HEADER_INSTALL
+    CholeskySolver.hpp
+    CholeskySolverImpl.hpp
+    CholeskySolverCpu.hpp
+    CholeskySolverCuda.hpp
+    CholeskySolverHip.hpp
+)
+
+add_library(resolve_hykkt_chol SHARED ${HyKKT_CHOL_SRC})
+
+target_link_libraries(resolve_hykkt_chol PUBLIC ${suitesparse_cholmod})
+target_include_directories(resolve_hykkt_chol PUBLIC ${SUITESPARSE_INCLUDE_DIR})
+
+# Link to CUDA ReSolve backend if CUDA is support enabled
+if (RESOLVE_USE_CUDA)
+  target_sources(resolve_hykkt_chol PRIVATE ${HyKKT_CUDASDK_SRC})
+  target_link_libraries(resolve_hykkt_chol PUBLIC resolve_backend_cuda)
+endif()
+
+if (RESOLVE_USE_HIP)
+  target_sources(resolve_hykkt_chol PRIVATE ${HyKKT_ROCM_SRC})
+  target_link_libraries(resolve_hykkt_chol PUBLIC resolve_backend_hip)
+endif()
+
+# Link to dummy device backend if GPU support is not enabled
+if (NOT RESOLVE_USE_GPU)
+  target_link_libraries(resolve_hykkt_chol PUBLIC resolve_backend_cpu)
+endif()
+
+target_link_libraries(resolve_hykkt_chol PUBLIC resolve_workspace
+                                                resolve_vector
+                                                resolve_matrix
+                                                resolve_logger)
+
+target_include_directories(resolve_hykkt_chol INTERFACE
+    $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}>
+    $<INSTALL_INTERFACE:include>
+)
+
+install(FILES ${HyKKT_CHOL_HEADER_INSTALL} DESTINATION include/resolve/hykkt)

resolve/hykkt/cholesky/CholeskySolver.cpp

@@ -0,0 +1,109 @@
+/**
+ * @file CholeskySolver.cpp
+ * @author Adham Ibrahim (ibrahimas@ornl.gov)
+ * @brief Cholesky decomposition solver implementation
+ */
+
+#include "CholeskySolver.hpp"
+
+#include "CholeskySolverCpu.hpp"
+#ifdef RESOLVE_USE_CUDA
+#include "CholeskySolverCuda.hpp"
+#elif defined(RESOLVE_USE_HIP)
+#include "CholeskySolverHip.hpp"
+#endif
+
+namespace ReSolve
+{
+  using real_type = ReSolve::real_type;
+  using out       = ReSolve::io::Logger;
+
+  namespace hykkt
+  {
+    /**
+     * @brief Cholesky Solver constructor
+     * @param[in] memspace - memory space to use for computations
+     */
+    CholeskySolver::CholeskySolver(memory::MemorySpace memspace)
+      : memspace_(memspace)
+    {
+      if (memspace_ == memory::HOST)
+      {
+        impl_ = new CholeskySolverCpu();
+      }
+      else
+      {
+#ifdef RESOLVE_USE_CUDA
+        impl_ = new CholeskySolverCuda();
+#elif defined(RESOLVE_USE_HIP)
+        impl_ = new CholeskySolverHip();
+#else
+        out::error() << "Memory space set to DEVICE, though no GPU support enabled. Must enable RESOLVE_USE_CUDA or RESOLVE_USE_HIP.\n";
+        exit(1);
+#endif
+      }
+    }
+
+    /**
+     * @brief Cholesky Solver destructor
+     */
+    CholeskySolver::~CholeskySolver()
+    {
+      delete impl_;
+    }
+
+    /**
+     * @brief Loads or reloads matrix pointer to the solver
+     * @param[in] A - pointer to the matrix in CSR format
+     */
+    void CholeskySolver::addMatrixInfo(matrix::Csr* A)
+    {
+      A_ = A;
+      impl_->addMatrixInfo(A);
+    }
+
+    /**
+     * @brief Performs symbolic analysis. This need only be called once
+     *        as long as the sparsity pattern does not change.
+     */
+    void CholeskySolver::symbolicAnalysis()
+    {
+      impl_->symbolicAnalysis();
+    }
+
+    /**
+     * @brief Sets the pivot tolerance for the solver.
+     *
+     * This is only used in the CUDA implementation. For other backends,
+     * it is ignored.
+     *
+     * @param[in] tol - pivot tolerance value
+     */
+    void CholeskySolver::setPivotTolerance(real_type tol)
+    {
+      tol_ = tol;
+    }
+
+    /**
+     * @brief Performs numerical factorization.
+     */
+    void CholeskySolver::numericalFactorization()
+    {
+      impl_->numericalFactorization(tol_);
+    }
+
+    /**
+     * @brief Solves the linear system Ax = b and stores the result in x.
+     *
+     * @pre The vector x is allocated in the given memspace.
+     *
+     * @param[out] x - pointer to the solution vector
+     * @param[in] b - pointer to the right-hand side vector
+     */
+    void CholeskySolver::solve(vector::Vector* x, vector::Vector* b)
+    {
+      impl_->solve(x, b);
+      x->setDataUpdated(memspace_);
+    }
+  } // namespace hykkt
+} // namespace ReSolve

resolve/hykkt/cholesky/CholeskySolver.hpp

@@ -0,0 +1,37 @@
+/**
+ * @file CholeskySolver.hpp
+ * @author Adham Ibrahim (ibrahimas@ornl.gov)
+ * @brief Cholesky decomposition solver header
+ */
+
+#pragma once
+#include "CholeskySolverImpl.hpp"
+#include <resolve/MemoryUtils.hpp>
+#include <resolve/matrix/Csr.hpp>
+#include <resolve/vector/Vector.hpp>
+
+namespace ReSolve
+{
+  namespace hykkt
+  {
+    class CholeskySolver
+    {
+    public:
+      CholeskySolver(memory::MemorySpace memspace);
+      ~CholeskySolver();
+
+      void addMatrixInfo(matrix::Csr* A);
+      void symbolicAnalysis();
+      void setPivotTolerance(real_type tol);
+      void numericalFactorization();
+      void solve(vector::Vector* x, vector::Vector* b);
+
+    private:
+      memory::MemorySpace memspace_;
+
+      matrix::Csr*        A_;
+      real_type           tol_ = 1e-12;
+      CholeskySolverImpl* impl_;
+    };
+  } // namespace hykkt
+} // namespace ReSolve