cuDSS implementations of Cholesky solver and refactorization

CHANGELOG.md

@@ -48,6 +48,7 @@
 
 - Fixed GPU refactorization allocation and synchronization errors, uninitialized GMRES example initial guesses, and CUDA ILU0 failures on stored numerical zero pivots.
 
+- Added cuDSS implementation of Cholesky solver (HyKKT) and refactorization
 - Optimized CGS2 coefficient accumulation in `GramSchmidt` by eliminating an unnecessary device-to-host synchronization.
 
 ## Changes to Re::Solve in release 0.99.2

CMakeLists.txt

@@ -44,6 +44,7 @@ option(RESOLVE_USE_UBSAN "Enable the undefined behavior sanitizer" OFF)
 option(RESOLVE_USE_DOXYGEN "Use Doxygen to generate Re::Solve documentation"
        OFF
 )
+option(RESOLVE_USE_CUDSS "Use the CUDA cuDSS libary" OFF)
 set(RESOLVE_CTEST_OUTPUT_DIR
     ${PROJECT_BINARY_DIR}
     CACHE PATH "Directory where CTest outputs are saved"
@@ -183,6 +184,19 @@ else()
   message(STATUS "Not using HIP")
 endif(RESOLVE_USE_HIP)
 
+if(RESOLVE_USE_CUDSS)
+  find_package(cudss)
+  if(NOT (cudss_FOUND AND RESOLVE_USE_CUDA))
+    message(STATUS "Cannot find cuDSS, disabling it ...")
+    set(RESOLVE_USE_CUDSS
+        OFF
+        CACHE BOOL FORCE
+    )
+  endif()
+else()
+  message(STATUS "Not using CUDA cuDSS")
+endif()
+
 # The binary dir is already a global include directory
 configure_file(
   ${CMAKE_SOURCE_DIR}/resolve/resolve_defs.hpp.in

cmake/ReSolveFindCudaLibraries.cmake

@@ -10,6 +10,10 @@ target_link_libraries(
                          CUDA::cudart
 )
 
+if(RESOLVE_USE_CUDSS)
+  target_link_libraries(resolve_cuda INTERFACE cudss)
+endif()
+
 if(RESOLVE_USE_PROFILING)
   target_link_libraries(resolve_cuda INTERFACE CUDA::nvToolsExt)
 endif()

examples/CMakeLists.txt

@@ -41,6 +41,13 @@ if(RESOLVE_USE_KLU)
     add_executable(gluRefactor.exe gluRefactor.cpp)
     target_link_libraries(gluRefactor.exe PRIVATE ReSolve)
 
+    if(RESOLVE_USE_CUDSS)
+      # Build an example with a configurable and portable system solver
+      # implemented with cuDSS
+      add_executable(cudssRefactor.exe cudssRefactor.cpp)
+      target_link_libraries(cudssRefactor.exe PRIVATE ReSolve)
+    endif(RESOLVE_USE_CUDSS)
+
   endif(RESOLVE_USE_CUDA)
 
 endif(RESOLVE_USE_KLU)

examples/cudssRefactor.cpp

@@ -0,0 +1,301 @@
+#include <chrono>
+#include <cmath>
+#include <cuda_runtime.h>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+#include <string>
+
+#include "ExampleHelper.hpp"
+#include <resolve/LinSolverDirectKLU.hpp>
+#include <resolve/Profiling.hpp>
+#include <resolve/SystemSolver.hpp>
+#include <resolve/matrix/Csr.hpp>
+#include <resolve/matrix/MatrixHandler.hpp>
+#include <resolve/matrix/io.hpp>
+#include <resolve/utilities/params/CliOptions.hpp>
+#include <resolve/vector/Vector.hpp>
+#include <resolve/vector/VectorHandler.hpp>
+#include <resolve/workspace/LinAlgWorkspace.hpp>
+
+/// Prints help message describing system usage.
+void printHelpInfo()
+{
+  std::cout << "\ncudssRefactor.exe loads from files and solves a series of linear systems.\n\n";
+  std::cout << "System matrices are in files with names <pathname>XX.mtx, where XX are\n";
+  std::cout << "consecutive integer numbers 00, 01, 02, ...\n\n";
+  std::cout << "System right hand side vectors are stored in files with matching numbering\n";
+  std::cout << "and file extension.\n\n";
+  std::cout << "Usage:\n\t./";
+  std::cout << "cudssRefactor.exe -m <matrix pathname> -r <rhs pathname> -n <number of systems>\n\n";
+  std::cout << "Optional features:\n";
+  std::cout << "\t-e <ext> \tSelects custom extension for input files (default 'mtx').\n";
+  std::cout << "\t-h\tPrints this message.\n";
+  std::cout << "\t-i\tEnables iterative refinement.\n";
+  std::cout << "\t-t, --time\tPrint solve timing for each linear system.\n\n";
+}
+
+using namespace ReSolve::constants;
+
+/// Prototype of the example function
+template <class workspace_type>
+static int cudssRefactor(int argc, char* argv[]);
+
+/// Main function selects example to be run.
+int main(int argc, char* argv[])
+{
+  ReSolve::CliOptions options(argc, argv);
+
+  // If help flag is passed, print help message and return
+  bool is_help = options.hasKey("-h");
+  if (is_help)
+  {
+    printHelpInfo();
+    return 0;
+  }
+
+  cudssRefactor<ReSolve::LinAlgWorkspaceCUDA>(argc, argv);
+
+  return 0;
+}
+
+/**
+ * @brief Example of using system solvers on GPU
+ *
+ * @tparam workspace_type - Type of the workspace to use
+ * @param[in] argc - Number of command line arguments
+ * @param[in] argv - Command line arguments
+ * @return 0 if the example ran successfully, -1 otherwise
+ */
+template <class workspace_type>
+int cudssRefactor(int argc, char* argv[])
+{
+  // Use the same data types as those you specified in ReSolve build.
+  using namespace ReSolve::examples;
+  using namespace ReSolve;
+  using index_type  = ReSolve::index_type;
+  using vector_type = ReSolve::vector::Vector;
+
+  CliOptions options(argc, argv);
+
+  bool is_help = options.hasKey("-h");
+  if (is_help)
+  {
+    printHelpInfo();
+    return 0;
+  }
+
+  bool is_iterative_refinement = options.hasKey("-i");
+  bool is_timing               = options.hasKey("-t") || options.hasKey("--time");
+
+  index_type num_systems = 0;
+  auto       opt         = options.getParamFromKey("-n");
+  if (opt)
+  {
+    num_systems = atoi((opt->second).c_str());
+  }
+  else
+  {
+    std::cout << "Incorrect input!\n";
+    printHelpInfo();
+  }
+
+  std::string matrix_pathname("");
+  opt = options.getParamFromKey("-m");
+  if (opt)
+  {
+    matrix_pathname = opt->second;
+  }
+  else
+  {
+    std::cout << "Incorrect input!\n";
+    printHelpInfo();
+  }
+
+  std::string rhs_pathname("");
+  opt = options.getParamFromKey("-r");
+  if (opt)
+  {
+    rhs_pathname = opt->second;
+  }
+  else
+  {
+    std::cout << "Incorrect input!\n";
+    printHelpInfo();
+  }
+
+  std::string file_extension("");
+  opt = options.getParamFromKey("-e");
+  if (opt)
+  {
+    file_extension = opt->second;
+  }
+  else
+  {
+    file_extension = "mtx";
+  }
+
+  std::cout << "Family matrix file name: " << matrix_pathname
+            << ", total number of matrices: " << num_systems << "\n"
+            << "Family rhs file name: " << rhs_pathname
+            << ", total number of RHSes: " << num_systems << "\n";
+
+  int status = 0;
+
+  workspace_type workspace;
+  workspace.initializeHandles();
+
+  // Create a helper object (for computing errors, printing summaries, etc.)
+  ExampleHelper<workspace_type> helper(workspace);
+  std::cout << "cudssRefactor with CUDA backend\n";
+
+  MatrixHandler matrix_handler(&workspace);
+  VectorHandler vector_handler(&workspace);
+
+  // Pointers to the linear system
+  matrix::Csr* A       = nullptr; // Pointer to the system matrix
+  vector_type* vec_rhs = nullptr; // Pointer for the right hand side vector
+  vector_type* vec_x   = nullptr; // Pointer for the solution vector
+
+  // Create system solver
+  ReSolve::SystemSolver solver(&workspace,
+                               "klu",     // factorization
+                               "cudssrf", // refactorization
+                               "cudssrf", // triangular solve
+                               "none",    // preconditioner (always 'none' here)
+                               "none");   // iterative refinement
+
+  if (is_iterative_refinement)
+  {
+    solver.setRefinementMethod("fgmres", "cgs2");
+    solver.getIterativeSolver().setCliParam("restart", "100");
+    solver.getIterativeSolver().setTol(1e-17);
+  }
+
+  RESOLVE_RANGE_PUSH(__FUNCTION__);
+  for (int i = 1; i <= num_systems; ++i)
+  {
+    std::cout << "System " << i << ":\n";
+    RESOLVE_RANGE_PUSH("File input");
+    std::ostringstream matname;
+    std::ostringstream rhsname;
+    matname << matrix_pathname << std::setfill('0') << std::setw(2) << i << "." << file_extension;
+    rhsname << rhs_pathname << std::setfill('0') << std::setw(2) << i << "." << file_extension;
+    std::string matrix_pathname_full = matname.str();
+    std::string rhs_pathname_full    = rhsname.str();
+
+    // Read matrix and right-hand-side vector
+    std::ifstream mat_file(matrix_pathname_full);
+    if (!mat_file.is_open())
+    {
+      std::cout << "Failed to open file " << matrix_pathname_full << "\n";
+      return 1;
+    }
+    std::ifstream rhs_file(rhs_pathname_full);
+    if (!rhs_file.is_open())
+    {
+      std::cout << "Failed to open file " << rhs_pathname_full << "\n";
+      return 1;
+    }
+
+    // Refactorization is LU-based, so need to expand symmetric matrices
+    bool is_expand_symmetric = true;
+    if (i == 1)
+    {
+      A       = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric);
+      vec_rhs = ReSolve::io::createVectorFromFile(rhs_file);
+      vec_x   = new vector_type(A->getNumRows());
+      vec_x->allocateAll(memory::DEVICE);
+    }
+    else
+    {
+      ReSolve::io::updateMatrixFromFile(mat_file, A);
+      ReSolve::io::updateVectorFromFile(rhs_file, vec_rhs);
+    }
+
+    mat_file.close();
+    rhs_file.close();
+
+    // Ensure matrix data is synced to the device before any GPU operations
+    if (i == 1)
+    {
+      A->allocateMatrixData(memory::DEVICE);
+      vec_rhs->allocate(memory::DEVICE);
+    }
+    A->syncData(memory::DEVICE);
+    vec_rhs->syncData(memory::DEVICE);
+
+    std::cout << "CSR matrix loaded. Expanded NNZ: " << A->getNnz() << std::endl;
+    printSystemInfo(matrix_pathname_full, A);
+
+    double solve_time_ms = 0.0;
+    cudaDeviceSynchronize();
+    auto solve_start = std::chrono::high_resolution_clock::now();
+
+    // Now call direct solver
+    if (i == 1)
+    {
+      // Set matrix in solver after the initial matrix is loaded
+      status = solver.setMatrix(A);
+      if (status != 0)
+      {
+        std::cout << "Failed to set matrix in solver. Status: " << status << std::endl;
+        return 1;
+      }
+
+      // Analysis (symbolic factorization)
+      status = solver.analyze();
+      std::cout << "Analysis on the host status: " << status << std::endl;
+
+      // Numeric factorization on the host
+      status = solver.factorize();
+      std::cout << "Numeric factorization on the host status: " << status << std::endl;
+      // Set up refactorization solver
+      status = solver.refactorizationSetup();
+      std::cout << "Refactorization setup status: " << status << std::endl;
+    }
+    else
+    {
+      // Refactorize on the device
+      status = solver.refactorize();
+      std::cout << "Refactorization on the device status: " << status << std::endl;
+    }
+
+    status = solver.solve(vec_rhs, vec_x);
+
+    if (i == 1)
+    {
+      vec_x->syncData(memory::DEVICE);
+    }
+
+    cudaDeviceSynchronize();
+    auto solve_end = std::chrono::high_resolution_clock::now();
+    solve_time_ms  = std::chrono::duration<double, std::milli>(solve_end - solve_start).count();
+    std::cout << "Triangular solve status: " << status << std::endl;
+
+    // Print summary of results
+    helper.printShortSummary(A, vec_rhs, vec_x);
+    if ((i > 1) && is_iterative_refinement)
+    {
+      helper.printIrSummary(&(solver.getIterativeSolver()));
+    }
+
+    if (is_timing)
+    {
+      std::cout << "TIMING,"
+                << "cudssRefactor,"
+                << "CUDA" << ","
+                << is_iterative_refinement << ","
+                << i << ","
+                << solve_time_ms
+                << std::endl;
+    }
+  }
+
+  // Delete objects created on heap
+  delete A;
+  delete vec_x;   // Delete the solution vector
+  delete vec_rhs; // Delete the RHS vector
+
+  return 0;
+}

resolve/CMakeLists.txt

@@ -33,6 +33,9 @@ set(ReSolve_CUDASDK_SRC
     LinSolverDirectCuSolverGLU.cpp LinSolverDirectCuSolverRf.cpp
     LinSolverDirectCuSparseILU0.cpp
 )
+if(RESOLVE_USE_CUDSS)
+  list(APPEND ReSolve_CUDASDK_SRC LinSolverDirectCuDssRf.cpp)
+endif()
 
 # C++ code that links to ROCm libraries
 set(ReSolve_ROCM_SRC LinSolverDirectRocSolverRf.cpp
@@ -64,6 +67,9 @@ set(ReSolve_CUDA_HEADER_INSTALL
     LinSolverDirectCuSolverGLU.hpp LinSolverDirectCuSolverRf.hpp
     LinSolverDirectCuSparseILU0.hpp
 )
+if(RESOLVE_USE_CUDSS)
+  list(APPEND ReSolve_CUDA_HEADER_INSTALL LinSolverDirectCuDssRf.hpp)
+endif()
 
 set(ReSolve_ROCM_HEADER_INSTALL LinSolverDirectRocSolverRf.hpp
                                 LinSolverDirectRocSparseILU0.hpp