From 490439db18cc3ff1c801c53508dacb0df680c26f Mon Sep 17 00:00:00 2001 From: Andrew Xu Date: Tue, 21 Jul 2026 16:19:45 +0000 Subject: [PATCH 01/18] cuDSS implementation of Cholesky --- cmake/ReSolveFindCudaLibraries.cmake | 3 +- resolve/hykkt/cholesky/CholeskySolverCuda.cpp | 239 +++++++++++++----- resolve/hykkt/cholesky/CholeskySolverCuda.hpp | 19 +- 3 files changed, 197 insertions(+), 64 deletions(-) diff --git a/cmake/ReSolveFindCudaLibraries.cmake b/cmake/ReSolveFindCudaLibraries.cmake index dc9e8aedc..be5137deb 100644 --- a/cmake/ReSolveFindCudaLibraries.cmake +++ b/cmake/ReSolveFindCudaLibraries.cmake @@ -4,10 +4,11 @@ add_library(resolve_cuda INTERFACE) find_package(CUDAToolkit REQUIRED) +find_package(cudss REQUIRED) target_link_libraries( resolve_cuda INTERFACE CUDA::cusolver CUDA::cublas CUDA::cusparse - CUDA::cudart + CUDA::cudart cudss ) if(RESOLVE_USE_PROFILING) diff --git a/resolve/hykkt/cholesky/CholeskySolverCuda.cpp b/resolve/hykkt/cholesky/CholeskySolverCuda.cpp index 13cf29514..6d1e7a822 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCuda.cpp +++ b/resolve/hykkt/cholesky/CholeskySolverCuda.cpp @@ -13,87 +13,167 @@ namespace ReSolve namespace hykkt { - CholeskySolverCuda::CholeskySolverCuda() + CholeskySolverCuda::CholeskySolverCuda(bool use_cudss) + : use_cudss_(use_cudss) { - cusolverSpCreate(&cusolverHandle_); - cusparseCreateMatDescr(&descrA_); - cusolverSpCreateCsrcholInfo(&factorizationInfo_); - buffer_ = nullptr; + if (use_cudss_) + { + cudssCreate(&cudss_handle_); + cudssConfigCreate(&cudss_config_); + cudssDataCreate(cudss_handle_, &cudss_data_); + } + else + { + cusolverSpCreate(&cusolverHandle_); + cusparseCreateMatDescr(&descr_A_cusolver_); + cusolverSpCreateCsrcholInfo(&factorizationInfo_); + buffer_ = nullptr; + } } CholeskySolverCuda::~CholeskySolverCuda() { - cusolverSpDestroy(cusolverHandle_); - cusparseDestroyMatDescr(descrA_); - cusolverSpDestroyCsrcholInfo(factorizationInfo_); - mem_.deleteOnDevice(buffer_); + if (use_cudss_) + { + cudssDataDestroy(cudss_handle_, cudss_data_); + cudssConfigDestroy(cudss_config_); + cudssDestroy(cudss_handle_); + cudssMatrixDestroy(descr_A_cudss_); + cudssMatrixDestroy(descr_b_); + cudssMatrixDestroy(descr_x_); + } + else + { + cusolverSpDestroy(cusolverHandle_); + cusparseDestroyMatDescr(descr_A_cusolver_); + cusolverSpDestroyCsrcholInfo(factorizationInfo_); + mem_.deleteOnDevice(buffer_); + } } void CholeskySolverCuda::addMatrixInfo(matrix::Csr* A) { A_ = A; + + if (use_cudss_) + { + cudssMatrixCreateCsr(&descr_A_cudss_, + A_->getNumRows(), + A_->getNumColumns(), + A_->getNnz(), + A_->getRowData(memory::DEVICE), + nullptr, // Row end offsets (null for standard CSR) + A_->getColData(memory::DEVICE), + A_->getValues(memory::DEVICE), + CUDA_R_32I, + CUDA_R_64F, + CUDSS_MTYPE_SPD, + CUDSS_MVIEW_LOWER, + CUDSS_BASE_ZERO); + } } /** * @brief Perform symbolic analysis for the Cholesky factorization - * - * Uses the `cusolverSpXcsrcholAnalysis` routine. */ void CholeskySolverCuda::symbolicAnalysis() { - cusolverSpXcsrcholAnalysis(cusolverHandle_, - A_->getNumRows(), - A_->getNnz(), - descrA_, - A_->getRowData(memory::DEVICE), - A_->getColData(memory::DEVICE), - factorizationInfo_); - // Calculate size of buffer needed - size_t internalDataBytes = 0; - size_t workspaceBytes = 0; - cusolverSpDcsrcholBufferInfo(cusolverHandle_, - A_->getNumRows(), - A_->getNnz(), - descrA_, - A_->getValues(memory::DEVICE), - A_->getRowData(memory::DEVICE), - A_->getColData(memory::DEVICE), - factorizationInfo_, - &internalDataBytes, - &workspaceBytes); - if (buffer_ != nullptr) + if (use_cudss_) { - mem_.deleteOnDevice(buffer_); + cudssMatrixCreateDn(&descr_b_, + A_->getNumRows(), + 1, + A_->getNumRows(), + nullptr, + CUDA_R_64F, + CUDSS_LAYOUT_COL_MAJOR); + cudssMatrixCreateDn(&descr_x_, + A_->getNumRows(), + 1, + A_->getNumRows(), + nullptr, + CUDA_R_64F, + CUDSS_LAYOUT_COL_MAJOR); + cudssExecute(cudss_handle_, + CUDSS_PHASE_ANALYSIS, + cudss_config_, + cudss_data_, + descr_A_cudss_, + descr_x_, + descr_b_); + } + else + { + cusolverSpXcsrcholAnalysis(cusolverHandle_, + A_->getNumRows(), + A_->getNnz(), + descr_A_cusolver_, + A_->getRowData(memory::DEVICE), + A_->getColData(memory::DEVICE), + factorizationInfo_); + // Calculate size of buffer needed + size_t internalDataBytes = 0; + size_t workspaceBytes = 0; + cusolverSpDcsrcholBufferInfo(cusolverHandle_, + A_->getNumRows(), + A_->getNnz(), + descr_A_cusolver_, + A_->getValues(memory::DEVICE), + A_->getRowData(memory::DEVICE), + A_->getColData(memory::DEVICE), + factorizationInfo_, + &internalDataBytes, + &workspaceBytes); + if (buffer_ != nullptr) + { + mem_.deleteOnDevice(buffer_); + } + mem_.allocateBufferOnDevice(&buffer_, workspaceBytes); } - mem_.allocateBufferOnDevice(&buffer_, workspaceBytes); } /** * @brief Perform numerical factorization for the Cholesky factorization - * - * Uses the `cusolverSpDcsrcholFactor` routine. - * + * * @param[in] tol - Tolerance for zero pivot detection. */ void CholeskySolverCuda::numericalFactorization(real_type tol) { - int singularity = 0; - cusolverSpDcsrcholFactor(cusolverHandle_, - A_->getNumRows(), - A_->getNnz(), - descrA_, - A_->getValues(memory::DEVICE), - A_->getRowData(memory::DEVICE), - A_->getColData(memory::DEVICE), - factorizationInfo_, - buffer_); - cusolverSpDcsrcholZeroPivot(cusolverHandle_, - factorizationInfo_, - tol, - &singularity); - if (singularity >= 0) + if (use_cudss_) { - out::error() << "Cholesky factorization failed with singularity at index: " << singularity << "\n"; + cudssConfigSet(cudss_config_, CUDSS_CONFIG_PIVOT_EPSILON, &tol, sizeof(real_type)); + cudssStatus_t status = cudssExecute(cudss_handle_, + CUDSS_PHASE_FACTORIZATION, + cudss_config_, + cudss_data_, + descr_A_cudss_, + descr_x_, + descr_b_); + if (status != CUDSS_STATUS_SUCCESS) + { + out::error() << "Cholesky factorization failed with status: " << status << "\n"; + } + } + else + { + int singularity = 0; + cusolverSpDcsrcholFactor(cusolverHandle_, + A_->getNumRows(), + A_->getNnz(), + descr_A_cusolver_, + A_->getValues(memory::DEVICE), + A_->getRowData(memory::DEVICE), + A_->getColData(memory::DEVICE), + factorizationInfo_, + buffer_); + cusolverSpDcsrcholZeroPivot(cusolverHandle_, + factorizationInfo_, + tol, + &singularity); + if (singularity >= 0) + { + out::error() << "Cholesky factorization failed with singularity at index: " << singularity << "\n"; + } } } @@ -107,13 +187,52 @@ namespace ReSolve */ void CholeskySolverCuda::solve(vector::Vector* x, vector::Vector* b) { - cusolverSpDcsrcholSolve(cusolverHandle_, - A_->getNumRows(), - b->getData(memory::DEVICE), - x->getData(memory::DEVICE), - factorizationInfo_, - buffer_); + if (use_cudss_) + { + if (descr_b_) + { + cudssMatrixDestroy(descr_b_); + } + if (descr_x_) + { + cudssMatrixDestroy(descr_x_); + } + + cudssMatrixCreateDn(&descr_b_, + b->getSize(), + b->getNumVectors(), + b->getSize(), + b->getData(memory::DEVICE), + CUDA_R_64F, + CUDSS_LAYOUT_COL_MAJOR); + cudssMatrixCreateDn(&descr_x_, + b->getSize(), + b->getNumVectors(), + b->getSize(), + x->getData(memory::DEVICE), + CUDA_R_64F, + CUDSS_LAYOUT_COL_MAJOR); + + cudssStatus_t status = cudssExecute(cudss_handle_, CUDSS_PHASE_SOLVE, cudss_config_, cudss_data_, descr_A_cudss_, descr_x_, descr_b_); + if (status != CUDSS_STATUS_SUCCESS) + { + out::error() << "cuDSS triangular solve failed with status: " << status << "\n"; + } + } + else + { + for (index_type i = 0; i < b->getNumVectors(); i++) + { + cusolverSpDcsrcholSolve(cusolverHandle_, + A_->getNumRows(), + b->getData(i, memory::DEVICE), + x->getData(i, memory::DEVICE), + factorizationInfo_, + buffer_); + } + } + x->setDataUpdated(memory::DEVICE); } } // namespace hykkt -} // namespace ReSolve +} // namespace ReSolve \ No newline at end of file diff --git a/resolve/hykkt/cholesky/CholeskySolverCuda.hpp b/resolve/hykkt/cholesky/CholeskySolverCuda.hpp index 29cce95de..0488ffcad 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCuda.hpp +++ b/resolve/hykkt/cholesky/CholeskySolverCuda.hpp @@ -1,9 +1,13 @@ /** * @file CholeskySolverCuda.hpp * @author Adham Ibrahim (ibrahimas@ornl.gov) - * @brief Header for CUDA implementation of Cholesky Solver + * @brief Header for CUDA implementation of Cholesky Solver using cuDSS */ +#pragma once + +#include + #include #include #include @@ -18,7 +22,7 @@ namespace ReSolve class CholeskySolverCuda : public CholeskySolverImpl { public: - CholeskySolverCuda(); + CholeskySolverCuda(bool use_cudss = true); ~CholeskySolverCuda(); void addMatrixInfo(matrix::Csr* A); @@ -31,9 +35,18 @@ namespace ReSolve matrix::Csr* A_; // pointer to the input matrix + bool use_cudss_; + + cudssHandle_t cudss_handle_; + cudssConfig_t cudss_config_; + cudssData_t cudss_data_; + cudssMatrix_t descr_A_cudss_; + cudssMatrix_t descr_b_; + cudssMatrix_t descr_x_; + // handle to the cuSPARSE library context cusolverSpHandle_t cusolverHandle_; - cusparseMatDescr_t descrA_; // descriptor for matrix A + cusparseMatDescr_t descr_A_cusolver_; // descriptor for matrix A csrcholInfo_t factorizationInfo_; // stores Cholesky factorization void* buffer_; // buffer for Cholesky factorization }; From 2e9ea510dc7df6c3de047370647d414692d45e69 Mon Sep 17 00:00:00 2001 From: Andrew Xu Date: Tue, 21 Jul 2026 20:40:51 +0000 Subject: [PATCH 02/18] cuDSS implementation of refactorization (rf) --- CHANGELOG.md | 2 + resolve/CMakeLists.txt | 4 +- resolve/LinSolverDirectCuDssRf.cpp | 464 ++++++++++++++++++++++++ resolve/LinSolverDirectCuDssRf.hpp | 83 +++++ resolve/SystemSolver.cpp | 20 +- tests/functionality/testRefactor.cpp | 8 +- tests/functionality/testSysRefactor.cpp | 19 +- 7 files changed, 588 insertions(+), 12 deletions(-) create mode 100644 resolve/LinSolverDirectCuDssRf.cpp create mode 100644 resolve/LinSolverDirectCuDssRf.hpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 126abed0a..c909fb0a0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -46,6 +46,8 @@ - Added line to testlog to tell the user to expect warnings as part of normal testing. - 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 refactorization + ## Changes to Re::Solve in release 0.99.2 ### Major Features diff --git a/resolve/CMakeLists.txt b/resolve/CMakeLists.txt index d54644ca8..21c6111e3 100644 --- a/resolve/CMakeLists.txt +++ b/resolve/CMakeLists.txt @@ -31,7 +31,7 @@ set(ReSolve_LUSOL_SRC LinSolverDirectLUSOL.cpp) # C++ code that links to CUDA SDK libraries set(ReSolve_CUDASDK_SRC LinSolverDirectCuSolverGLU.cpp LinSolverDirectCuSolverRf.cpp - LinSolverDirectCuSparseILU0.cpp + LinSolverDirectCuDssRf.cpp inSolverDirectCuSparseILU0.cpp ) # C++ code that links to ROCm libraries @@ -62,7 +62,7 @@ set(ReSolve_LUSOL_HEADER_INSTALL LinSolverDirectLUSOL.hpp) set(ReSolve_CUDA_HEADER_INSTALL LinSolverDirectCuSolverGLU.hpp LinSolverDirectCuSolverRf.hpp - LinSolverDirectCuSparseILU0.hpp + LinSolverDirectCuDssRf.hpp LinSolverDirectCuSparseILU0.hpp ) set(ReSolve_ROCM_HEADER_INSTALL LinSolverDirectRocSolverRf.hpp diff --git a/resolve/LinSolverDirectCuDssRf.cpp b/resolve/LinSolverDirectCuDssRf.cpp new file mode 100644 index 000000000..2d73e4b27 --- /dev/null +++ b/resolve/LinSolverDirectCuDssRf.cpp @@ -0,0 +1,464 @@ +#include "LinSolverDirectCuDssRf.hpp" + +#include +#include + +#include +#include +#include + +namespace ReSolve +{ + using out = io::Logger; + + /** + * @brief Constructor for LinSolverDirectCuDssRf + * + * @param workspace - pointer to the LinAlgWorkspaceCUDA object (not used) + */ + LinSolverDirectCuDssRf::LinSolverDirectCuDssRf(LinAlgWorkspaceCUDA* /* workspace */) + { + cudssCreate(&handle_cudss_); + cudssConfigCreate(&config_cudss_); + cudssDataCreate(handle_cudss_, &data_cudss_); + setup_completed_ = false; + initParamList(); + } + + /** + * @brief Destructor for LinSolverDirectCuDssRf + * + * Destroys the cuDssRf handle, config, and data, then deletes the permutation vectors + * from the device memory. + * + * @pre The cuDssRf handle, config, and data have been created. + * @post The cuDssRf handle, config, and data have destroyed. + * + * @pre The permutation vectors are allocated on the device. + * @post The permutation vectors are deleted from the device. + * + * @note The permutation vectors are not deleted from the host memory. + */ + LinSolverDirectCuDssRf::~LinSolverDirectCuDssRf() + { + cudssMatrixDestroy(descr_A_); + cudssDataDestroy(handle_cudss_, data_cudss_); + cudssConfigDestroy(config_cudss_); + cudssDestroy(handle_cudss_); + + mem_.deleteOnDevice(d_P_); + } + + /** + * @brief For backward compatibility. Only A and P are needed. + */ + int LinSolverDirectCuDssRf::setup(matrix::Sparse* A, + matrix::Sparse*, + matrix::Sparse*, + index_type* P, + index_type*, + vector_type*) + { + return setup(A, P); + } + + /** + * @brief Setup the cuDssRf factorization with factors already in CSR + * + * Sets up the cuDssRf factorization for the given matrix A. The + * permutation vector P is also set up. + * + * @param[in] A - pointer to the matrix A + * @param[in] P - pointer to the permutation vector P + * + * @pre The matrix A is in CSR format. + * @post Device storage for L and U is allocated and synchronized. + */ + int LinSolverDirectCuDssRf::setup(matrix::Sparse* A, + index_type* P) + { + assert(A->getSparseFormat() == matrix::Sparse::COMPRESSED_SPARSE_ROW && "Matrix A has to be in CSR format for cuDssRf input.\n"); + int error_sum = 0; + this->A_ = A; + index_type n = A_->getNumRows(); + index_type nnz = A_->getNnz(); + + if (setup_completed_) + { + cudssMatrixDestroy(descr_A_); + cudssDataDestroy(handle_cudss_, data_cudss_); + cudssConfigDestroy(config_cudss_); + + descr_A_ = nullptr; + cudssConfigCreate(&config_cudss_); + cudssDataCreate(handle_cudss_, &data_cudss_); + } + + if (d_P_ == nullptr) + { + mem_.allocateArrayOnDevice(&d_P_, n); + } + + mem_.copyArrayHostToDevice(d_P_, P, n); + + if (d_P_ != nullptr) + { + cudssAlgType_t reorder_alg = CUDSS_ALG_1; + cudssConfigSet(config_cudss_, CUDSS_CONFIG_REORDERING_ALG, &reorder_alg, sizeof(cudssAlgType_t)); + cudssDataSet(handle_cudss_, data_cudss_, CUDSS_DATA_USER_PERM, d_P_, n * sizeof(index_type)); + } + + cudssStatus_t status = cudssMatrixCreateCsr(&descr_A_, + n, + n, + nnz, + A_->getRowData(memory::DEVICE), + nullptr, + A_->getColData(memory::DEVICE), + A_->getValues(memory::DEVICE), + CUDA_R_32I, + CUDA_R_64F, + CUDSS_MTYPE_GENERAL, + CUDSS_MVIEW_FULL, + CUDSS_BASE_ZERO); + error_sum += status; + + status = cudssExecute(handle_cudss_, CUDSS_PHASE_ANALYSIS, config_cudss_, data_cudss_, descr_A_, nullptr, nullptr); + error_sum += status; + + setup_completed_ = true; + return error_sum; + } + + /** + * @brief Refactorizes the matrix A + * + * Refactorizes the matrix A using the cuDssRf handle. + * + * @pre The cuDssRf handle has been created. + * @pre The matrix A is in CSR format. + * @pre The permutation vectors P is allocated on the device. + * @pre Matrix A's data is on the device. + * + * @post The matrix A is refactorized. + * + * @return 0 if successful, 1 otherwise + */ + int LinSolverDirectCuDssRf::refactorize() + { + assert(A_ != nullptr && "Matrix A is null!"); + assert(A_->getNumRows() > 0 && "Matrix A must have positive row count!"); + assert(A_->getNnz() > 0 && "Matrix A must have positive nonzero count!"); + + return cudssExecute(handle_cudss_, + CUDSS_PHASE_FACTORIZATION, + config_cudss_, + data_cudss_, + descr_A_, + nullptr, + nullptr); + } + + /** + * @brief Solves the system of equations Ax=rhs + * + * Solves the system of equations Ax=rhs using the cuDssRf handle. + * The solution overwrites the right-hand side vector. + * + * @param[in,out] rhs - pointer to right-hand side vector, changes to solution + * + * @return 0 if successful, 1 otherwise + */ + int LinSolverDirectCuDssRf::solve(vector_type* rhs) + { + return solve(rhs, rhs); + } + + /** + * @brief solves the system of equations Ax=rhs + * + * Solves the system of equations Ax=rhs using the cuDssRf handle. + * The solution is stored in x. + * + * @param[in] rhs - pointer to right-hand side vector + * @param[out] x - pointer to solution vector + * + * @return 0 if successful, 1 otherwise + */ + int LinSolverDirectCuDssRf::solve(vector_type* rhs, vector_type* x) + { + index_type n = A_->getNumRows(); + + int error_sum = 0; + + cudssMatrix_t rhs_descr; + cudssMatrix_t x_descr; + error_sum += cudssMatrixCreateDn(&rhs_descr, n, 1, n, rhs->getData(memory::DEVICE), CUDA_R_64F, CUDSS_LAYOUT_COL_MAJOR); + error_sum += cudssMatrixCreateDn(&x_descr, n, 1, n, x->getData(memory::DEVICE), CUDA_R_64F, CUDSS_LAYOUT_COL_MAJOR); + error_sum += cudssExecute(handle_cudss_, CUDSS_PHASE_SOLVE, config_cudss_, data_cudss_, descr_A_, x_descr, rhs_descr); + + x->setDataUpdated(memory::DEVICE); + + cudssMatrixDestroy(rhs_descr); + cudssMatrixDestroy(x_descr); + + return error_sum; + } + + // For backward compatibility + int LinSolverDirectCuDssRf::setNumericalProperties(real_type nzero, real_type) + { + return setNumericalProperties(nzero); + } + + /** + * @brief Sets a flag threshold for zero pivots and a boost factor + * + * Sets the zero flagging threshold and boost factor for the cuDssRf handle. + * Must be called before setup() + * + * @param[in] nzero - zero flagging threshold + * @param[in] nboost - boost factor () + * + * @return 0 if successful, 1 otherwise + */ + int LinSolverDirectCuDssRf::setNumericalProperties(real_type nzero) + { + zero_pivot_ = nzero; + return cudssConfigSet(config_cudss_, + CUDSS_CONFIG_PIVOT_EPSILON, + &zero_pivot_, + sizeof(real_type)); + } + + /** + * @brief Sets the paramters from Cli to the cuDssRf handle + * + * @param[in] id - string ID for the parameter to set + * @param[in] value - string value for the parameter to set + * + * @return 0 if successful, 1 otherwise + */ + int LinSolverDirectCuDssRf::setCliParam(const std::string id, const std::string value) + { + switch (getParamId(id)) + { + case ZERO_PIVOT: + zero_pivot_ = atof(value.c_str()); + setNumericalProperties(zero_pivot_); + break; + case PIVOT_BOOST: + out::warning() << "Pivot boost is not implemented for cuDSS refactor.\n"; + break; + default: + std::cout << "Setting parameter failed!\n"; + } + return 0; + } + + /** + * @brief Placeholder function for now. + * + * The following switch (getParamId(Id)) cases always run the default and + * are currently redundant code (like an if (true)). + * In the future, they will be expanded to include more options. + * + * @param id - string ID for parameter to get. + * @return std::string Value of the string parameter to return. + */ + std::string LinSolverDirectCuDssRf::getCliParamString(const std::string id) const + { + switch (getParamId(id)) + { + default: + out::error() << "Trying to get unknown string parameter " << id << "\n"; + } + return ""; + } + + /** + * @brief Placeholder function for now. + * + * The following switch (getParamId(Id)) cases always run the default and + * are currently redundant code (like an if (true)). + * In the future, they will be expanded to include more options. + * + * @param id - string ID for parameter to get. + * @return int Value of the int parameter to return. + */ + index_type LinSolverDirectCuDssRf::getCliParamInt(const std::string id) const + { + switch (getParamId(id)) + { + default: + out::error() << "Trying to get unknown integer parameter " << id << "\n"; + } + return -1; + } + + real_type LinSolverDirectCuDssRf::getCliParamReal(const std::string id) const + { + switch (getParamId(id)) + { + case ZERO_PIVOT: + return zero_pivot_; + case PIVOT_BOOST: + out::warning() << "Pivot boost is not implemented for cuDSS refactor.\n"; + break; + default: + out::error() << "Trying to get unknown real parameter " << id << "\n"; + } + return std::numeric_limits::quiet_NaN(); + } + + /** + * @brief Placeholder function for now. + * + * The following switch (getParamId(Id)) cases always run the default and + * are currently redundant code (like an if (true)). + * In the future, they will be expanded to include more options. + * + * @param id - string ID for parameter to get. + * @return bool Value of the bool parameter to return. + */ + bool LinSolverDirectCuDssRf::getCliParamBool(const std::string id) const + { + switch (getParamId(id)) + { + default: + out::error() << "Trying to get unknown boolean parameter " << id << "\n"; + } + return false; + } + + /** + * @brief Prints the parameters from Cli to the console + */ + int LinSolverDirectCuDssRf::printCliParam(const std::string id) const + { + switch (getParamId(id)) + { + case ZERO_PIVOT: + std::cout << zero_pivot_ << "\n"; + break; + case PIVOT_BOOST: + out::warning() << "Pivot boost is not implemented for cuDSS refactor.\n"; + break; + default: + out::error() << "Trying to print unknown parameter " << id << "\n"; + return 1; + } + return 0; + } + + // + // Private methods + // + + /** + * @brief Set the zero pivot and pivot boost parameters + */ + void LinSolverDirectCuDssRf::initParamList() + { + params_list_["zero_pivot"] = ZERO_PIVOT; + params_list_["pivot_boost"] = PIVOT_BOOST; + } + + /** + * @brief Convert CSC to CSR matrix on the host + * + * @authors Slaven Peles , Daniel Reynolds (SMU), and + * David Gardner and Carol Woodward (LLNL) + * + * @param[in] A_csc - pointer to the CSC matrix + * @param[out] A_csr - pointer to an empty CSR matrix + * + * @return 0 if successful, 1 otherwise + */ + int LinSolverDirectCuDssRf::csc2csr(matrix::Csc* A_csc, matrix::Csr* A_csr) + { + // int error_sum = 0; TODO: Collect error output! + assert(A_csc->getNnz() == A_csr->getNnz()); + assert(A_csc->getNumRows() == A_csr->getNumRows()); + assert(A_csr->getNumColumns() == A_csc->getNumColumns()); + + A_csr->allocateMatrixData(memory::HOST); + + index_type nnz = A_csc->getNnz(); + index_type n = A_csc->getNumColumns(); + + index_type* rowIdxCsc = A_csc->getRowData(memory::HOST); + index_type* colPtrCsc = A_csc->getColData(memory::HOST); + real_type* valuesCsc = A_csc->getValues(memory::HOST); + + index_type* rowPtrCsr = A_csr->getRowData(memory::HOST); + index_type* colIdxCsr = A_csr->getColData(memory::HOST); + real_type* valuesCsr = A_csr->getValues(memory::HOST); + + // Set all CSR row pointers to zero + for (index_type i = 0; i <= n; ++i) + { + rowPtrCsr[i] = 0; + } + + // Set all CSR values and column indices to zero + for (index_type i = 0; i < nnz; ++i) + { + colIdxCsr[i] = 0; + valuesCsr[i] = 0.0; + } + + // Compute number of entries per row + for (index_type i = 0; i < nnz; ++i) + { + rowPtrCsr[rowIdxCsc[i]]++; + } + + // Compute cumualtive sum of nnz per row + for (index_type row = 0, rowsum = 0; row < n; ++row) + { + // Store value in row pointer to temp + index_type temp = rowPtrCsr[row]; + + // Copy cumulative sum to the row pointer + rowPtrCsr[row] = rowsum; + + // Update row sum + rowsum += temp; + } + rowPtrCsr[n] = nnz; + + for (index_type col = 0; col < n; ++col) + { + // Compute positions of column indices and values in CSR matrix and store them there + // Overwrites CSR row pointers in the process + // adding to them the number of elements in that row + for (index_type jj = colPtrCsc[col]; jj < colPtrCsc[col + 1]; jj++) + { + index_type row = rowIdxCsc[jj]; + index_type dest = rowPtrCsr[row]; + + colIdxCsr[dest] = col; + valuesCsr[dest] = valuesCsc[jj]; + + rowPtrCsr[row]++; + } + } + + // Restore CSR row pointer values + // All values in rowPtrCsr have shifted by the number of elements in that row + // for i>=1: new rowPtrCsr[i] = old rowPtrCsr[i-1] and new rowPtrCsr[0]=0 + for (index_type row = 0, last = 0; row <= n; row++) + { + index_type temp = rowPtrCsr[row]; + rowPtrCsr[row] = last; + last = temp; + } + + // Mark data on the host as updated + A_csr->setUpdated(memory::HOST); + + return 0; + } + +} // namespace ReSolve diff --git a/resolve/LinSolverDirectCuDssRf.hpp b/resolve/LinSolverDirectCuDssRf.hpp new file mode 100644 index 000000000..b3c0a5492 --- /dev/null +++ b/resolve/LinSolverDirectCuDssRf.hpp @@ -0,0 +1,83 @@ +#pragma once + +#include "Common.hpp" +#include +#include +#include + +namespace ReSolve +{ + // Forward declaration of vector::Vector class + namespace vector + { + class Vector; + } + + // Forward declaration of matrix::Sparse class + namespace matrix + { + class Sparse; + class Csr; + class Csc; + } // namespace matrix + + // Forward declaration of ReSolve handlers workspace + class LinAlgWorkspaceCUDA; + + class LinSolverDirectCuDssRf : public LinSolverDirect + { + using vector_type = vector::Vector; + + public: + LinSolverDirectCuDssRf(LinAlgWorkspaceCUDA* workspace = nullptr); + ~LinSolverDirectCuDssRf(); + + // For backward compatibility. + int setup(matrix::Sparse* A, + matrix::Sparse*, + matrix::Sparse*, + index_type* P, + index_type* Q, + vector_type* = nullptr) override; + + int setup(matrix::Sparse* A, + index_type* P); + + int refactorize() override; + int solve(vector_type* rhs, vector_type* x) override; + int solve(vector_type* rhs) override; // rhs overwritten by solution + + int setNumericalProperties(real_type nzero, real_type); // For backward compatibility + int setNumericalProperties(real_type nzero); + + int setCliParam(const std::string id, const std::string value) override; + std::string getCliParamString(const std::string id) const override; + index_type getCliParamInt(const std::string id) const override; + real_type getCliParamReal(const std::string id) const override; + bool getCliParamBool(const std::string id) const override; + int printCliParam(const std::string id) const override; + + private: + void initParamList(); + int csc2csr(matrix::Csc* A_csc, matrix::Csr* A_csr); + + private: + enum ParamaterIDs + { + ZERO_PIVOT = 0, + PIVOT_BOOST + }; + + real_type zero_pivot_{0.0}; ///< The value below which zero pivot is flagged. + + cudssHandle_t handle_cudss_{nullptr}; + cudssConfig_t config_cudss_{nullptr}; + cudssData_t data_cudss_{nullptr}; + cudssMatrix_t descr_A_{nullptr}; + + index_type* d_P_{nullptr}; + bool setup_completed_{false}; + + MemoryHandler mem_; ///< Device memory manager object + }; +} // namespace ReSolve \ No newline at end of file diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 0bd9fc43e..cec7f70f4 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -18,6 +18,7 @@ #ifdef RESOLVE_USE_CUDA #include #include +#include #include #include #endif @@ -316,6 +317,10 @@ namespace ReSolve else if (refactorizationMethod_ == "cusolverrf") { refactorizationSolver_ = new ReSolve::LinSolverDirectCuSolverRf(); + } + else if (refactorizationMethod_ == "cudssrf") + { + refactorizationSolver_ = new ReSolve::LinSolverDirectCuDssRf(); #endif #ifdef RESOLVE_USE_HIP } @@ -452,7 +457,7 @@ namespace ReSolve return factorizationSolver_->refactorize(); } - if (refactorizationMethod_ == "glu" || refactorizationMethod_ == "cusolverrf" || refactorizationMethod_ == "rocsolverrf") + if (refactorizationMethod_ == "glu" || refactorizationMethod_ == "cusolverrf" || refactorizationMethod_ == "cudssrf" || refactorizationMethod_ == "rocsolverrf") { is_solve_on_device_ = true; return refactorizationSolver_->refactorize(); @@ -499,13 +504,22 @@ namespace ReSolve is_solve_on_device_ = true; status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_); } - if (refactorizationMethod_ == "cusolverrf") + else if (refactorizationMethod_ == "cusolverrf") { status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_); LinSolverDirectCuSolverRf* Rf = dynamic_cast(refactorizationSolver_); Rf->setNumericalProperties(1e-14, 1e-1); + is_solve_on_device_ = false; + } + else if (refactorizationMethod_ == "cudssrf") + { + LinSolverDirectCuDssRf* Rf = dynamic_cast(refactorizationSolver_); + Rf->setNumericalProperties(1e-14, 1e-1); + + status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_); + is_solve_on_device_ = false; } #endif @@ -564,7 +578,7 @@ namespace ReSolve status += factorizationSolver_->solve(rhs, x); } - if (solveMethod_ == "glu" || solveMethod_ == "cusolverrf" || solveMethod_ == "rocsolverrf") + if (solveMethod_ == "glu" || solveMethod_ == "cusolverrf" || solveMethod_ == "cudssrf" || solveMethod_ == "rocsolverrf") { if (is_solve_on_device_) { diff --git a/tests/functionality/testRefactor.cpp b/tests/functionality/testRefactor.cpp index 2ba8e3fa1..61fdbd119 100644 --- a/tests/functionality/testRefactor.cpp +++ b/tests/functionality/testRefactor.cpp @@ -29,6 +29,7 @@ #ifdef RESOLVE_USE_CUDA #include #include +#include #endif #include "TestHelper.hpp" @@ -61,9 +62,12 @@ int main(int argc, char* argv[]) } if (rf_solver == "rf") { - std::string solver_name("cusolverRf"); + std::string cusolver_solver_name("cusolverRf"); error_sum += runTest(argc, argv, solver_name); + LinSolverDirectCuSolverRf>(argc, argv, cusolver_solver_name); + std::string cudss_solver_name("cusdssRf"); + error_sum += runTest(argc, argv, cudss_solver_name); } else { diff --git a/tests/functionality/testSysRefactor.cpp b/tests/functionality/testSysRefactor.cpp index 86d629057..b775964e1 100644 --- a/tests/functionality/testSysRefactor.cpp +++ b/tests/functionality/testSysRefactor.cpp @@ -27,6 +27,7 @@ #ifdef RESOLVE_USE_CUDA #include +#include #endif #ifdef RESOLVE_USE_HIP @@ -36,7 +37,7 @@ #include "TestHelper.hpp" template -static int runTest(int argc, char* argv[], std::string backend); +static int runTest(int argc, char* argv[], std::string backend, bool use_cudss); int main(int argc, char* argv[]) { @@ -46,18 +47,19 @@ int main(int argc, char* argv[]) // error_sum += runTest(argc, argv, "cpu"); #ifdef RESOLVE_USE_CUDA - error_sum += runTest(argc, argv, "cuda"); + error_sum += runTest(argc, argv, "cuda", false); + error_sum += runTest(argc, argv, "cuda", true); #endif #ifdef RESOLVE_USE_HIP - error_sum += runTest(argc, argv, "hip"); + error_sum += runTest(argc, argv, "hip", false); #endif return error_sum; } template -static int runTest(int argc, char* argv[], std::string backend) +static int runTest(int argc, char* argv[], std::string backend, bool use_cudss) { // Use ReSolve data types. using namespace ReSolve; @@ -118,7 +120,14 @@ static int runTest(int argc, char* argv[], std::string backend) std::string refactor("none"); if (backend == "cuda") { - refactor = "cusolverrf"; + if (use_cudss) + { + refactor = "cudssrf"; + } + else + { + refactor = "cusolverrf"; + } } else if (backend == "hip") { From ed653745729ef5d00790d95ed2ee2831a22c21b5 Mon Sep 17 00:00:00 2001 From: Andrew Xu Date: Tue, 21 Jul 2026 20:56:49 +0000 Subject: [PATCH 03/18] Update changelog and fix typo --- CHANGELOG.md | 2 +- resolve/CMakeLists.txt | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c909fb0a0..1467a75f9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -46,7 +46,7 @@ - Added line to testlog to tell the user to expect warnings as part of normal testing. - 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 refactorization +- Added cuDSS implementation of Cholesky solver (HyKKT) and refactorization ## Changes to Re::Solve in release 0.99.2 diff --git a/resolve/CMakeLists.txt b/resolve/CMakeLists.txt index 21c6111e3..bf8c008c7 100644 --- a/resolve/CMakeLists.txt +++ b/resolve/CMakeLists.txt @@ -31,7 +31,7 @@ set(ReSolve_LUSOL_SRC LinSolverDirectLUSOL.cpp) # C++ code that links to CUDA SDK libraries set(ReSolve_CUDASDK_SRC LinSolverDirectCuSolverGLU.cpp LinSolverDirectCuSolverRf.cpp - LinSolverDirectCuDssRf.cpp inSolverDirectCuSparseILU0.cpp + LinSolverDirectCuDssRf.cpp LinSolverDirectCuSparseILU0.cpp ) # C++ code that links to ROCm libraries From 7f33a4b6eedc31011849ab3508113418595b13ff Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Tue, 21 Jul 2026 20:58:24 +0000 Subject: [PATCH 04/18] Apply pre-commmit fixes --- resolve/LinSolverDirectCuDssRf.cpp | 24 +++--- resolve/LinSolverDirectCuDssRf.hpp | 11 +-- resolve/SystemSolver.cpp | 2 +- resolve/hykkt/cholesky/CholeskySolverCuda.cpp | 74 +++++++++---------- resolve/hykkt/cholesky/CholeskySolverCuda.hpp | 5 +- tests/functionality/testRefactor.cpp | 2 +- tests/functionality/testSysRefactor.cpp | 2 +- 7 files changed, 60 insertions(+), 60 deletions(-) diff --git a/resolve/LinSolverDirectCuDssRf.cpp b/resolve/LinSolverDirectCuDssRf.cpp index 2d73e4b27..c55f32aee 100644 --- a/resolve/LinSolverDirectCuDssRf.cpp +++ b/resolve/LinSolverDirectCuDssRf.cpp @@ -55,7 +55,7 @@ namespace ReSolve int LinSolverDirectCuDssRf::setup(matrix::Sparse* A, matrix::Sparse*, matrix::Sparse*, - index_type* P, + index_type* P, index_type*, vector_type*) { @@ -78,11 +78,11 @@ namespace ReSolve index_type* P) { assert(A->getSparseFormat() == matrix::Sparse::COMPRESSED_SPARSE_ROW && "Matrix A has to be in CSR format for cuDssRf input.\n"); - int error_sum = 0; - this->A_ = A; - index_type n = A_->getNumRows(); + int error_sum = 0; + this->A_ = A; + index_type n = A_->getNumRows(); index_type nnz = A_->getNnz(); - + if (setup_completed_) { cudssMatrixDestroy(descr_A_); @@ -151,12 +151,12 @@ namespace ReSolve assert(A_->getNnz() > 0 && "Matrix A must have positive nonzero count!"); return cudssExecute(handle_cudss_, - CUDSS_PHASE_FACTORIZATION, - config_cudss_, - data_cudss_, - descr_A_, - nullptr, - nullptr); + CUDSS_PHASE_FACTORIZATION, + config_cudss_, + data_cudss_, + descr_A_, + nullptr, + nullptr); } /** @@ -224,7 +224,7 @@ namespace ReSolve */ int LinSolverDirectCuDssRf::setNumericalProperties(real_type nzero) { - zero_pivot_ = nzero; + zero_pivot_ = nzero; return cudssConfigSet(config_cudss_, CUDSS_CONFIG_PIVOT_EPSILON, &zero_pivot_, diff --git a/resolve/LinSolverDirectCuDssRf.hpp b/resolve/LinSolverDirectCuDssRf.hpp index b3c0a5492..08032789c 100644 --- a/resolve/LinSolverDirectCuDssRf.hpp +++ b/resolve/LinSolverDirectCuDssRf.hpp @@ -1,7 +1,8 @@ #pragma once -#include "Common.hpp" #include + +#include "Common.hpp" #include #include @@ -36,8 +37,8 @@ namespace ReSolve int setup(matrix::Sparse* A, matrix::Sparse*, matrix::Sparse*, - index_type* P, - index_type* Q, + index_type* P, + index_type* Q, vector_type* = nullptr) override; int setup(matrix::Sparse* A, @@ -68,7 +69,7 @@ namespace ReSolve PIVOT_BOOST }; - real_type zero_pivot_{0.0}; ///< The value below which zero pivot is flagged. + real_type zero_pivot_{0.0}; ///< The value below which zero pivot is flagged. cudssHandle_t handle_cudss_{nullptr}; cudssConfig_t config_cudss_{nullptr}; @@ -80,4 +81,4 @@ namespace ReSolve MemoryHandler mem_; ///< Device memory manager object }; -} // namespace ReSolve \ No newline at end of file +} // namespace ReSolve diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index cec7f70f4..1e5895677 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -16,9 +16,9 @@ #include #ifdef RESOLVE_USE_CUDA +#include #include #include -#include #include #include #endif diff --git a/resolve/hykkt/cholesky/CholeskySolverCuda.cpp b/resolve/hykkt/cholesky/CholeskySolverCuda.cpp index 6d1e7a822..79b7f49c9 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCuda.cpp +++ b/resolve/hykkt/cholesky/CholeskySolverCuda.cpp @@ -58,18 +58,18 @@ namespace ReSolve if (use_cudss_) { cudssMatrixCreateCsr(&descr_A_cudss_, - A_->getNumRows(), - A_->getNumColumns(), - A_->getNnz(), - A_->getRowData(memory::DEVICE), - nullptr, // Row end offsets (null for standard CSR) - A_->getColData(memory::DEVICE), - A_->getValues(memory::DEVICE), - CUDA_R_32I, - CUDA_R_64F, - CUDSS_MTYPE_SPD, - CUDSS_MVIEW_LOWER, - CUDSS_BASE_ZERO); + A_->getNumRows(), + A_->getNumColumns(), + A_->getNnz(), + A_->getRowData(memory::DEVICE), + nullptr, // Row end offsets (null for standard CSR) + A_->getColData(memory::DEVICE), + A_->getValues(memory::DEVICE), + CUDA_R_32I, + CUDA_R_64F, + CUDSS_MTYPE_SPD, + CUDSS_MVIEW_LOWER, + CUDSS_BASE_ZERO); } } @@ -105,25 +105,25 @@ namespace ReSolve else { cusolverSpXcsrcholAnalysis(cusolverHandle_, - A_->getNumRows(), - A_->getNnz(), - descr_A_cusolver_, - A_->getRowData(memory::DEVICE), - A_->getColData(memory::DEVICE), - factorizationInfo_); + A_->getNumRows(), + A_->getNnz(), + descr_A_cusolver_, + A_->getRowData(memory::DEVICE), + A_->getColData(memory::DEVICE), + factorizationInfo_); // Calculate size of buffer needed size_t internalDataBytes = 0; size_t workspaceBytes = 0; cusolverSpDcsrcholBufferInfo(cusolverHandle_, - A_->getNumRows(), - A_->getNnz(), - descr_A_cusolver_, - A_->getValues(memory::DEVICE), - A_->getRowData(memory::DEVICE), - A_->getColData(memory::DEVICE), - factorizationInfo_, - &internalDataBytes, - &workspaceBytes); + A_->getNumRows(), + A_->getNnz(), + descr_A_cusolver_, + A_->getValues(memory::DEVICE), + A_->getRowData(memory::DEVICE), + A_->getColData(memory::DEVICE), + factorizationInfo_, + &internalDataBytes, + &workspaceBytes); if (buffer_ != nullptr) { mem_.deleteOnDevice(buffer_); @@ -134,7 +134,7 @@ namespace ReSolve /** * @brief Perform numerical factorization for the Cholesky factorization - * + * * @param[in] tol - Tolerance for zero pivot detection. */ void CholeskySolverCuda::numericalFactorization(real_type tol) @@ -158,14 +158,14 @@ namespace ReSolve { int singularity = 0; cusolverSpDcsrcholFactor(cusolverHandle_, - A_->getNumRows(), - A_->getNnz(), - descr_A_cusolver_, - A_->getValues(memory::DEVICE), - A_->getRowData(memory::DEVICE), - A_->getColData(memory::DEVICE), - factorizationInfo_, - buffer_); + A_->getNumRows(), + A_->getNnz(), + descr_A_cusolver_, + A_->getValues(memory::DEVICE), + A_->getRowData(memory::DEVICE), + A_->getColData(memory::DEVICE), + factorizationInfo_, + buffer_); cusolverSpDcsrcholZeroPivot(cusolverHandle_, factorizationInfo_, tol, @@ -235,4 +235,4 @@ namespace ReSolve x->setDataUpdated(memory::DEVICE); } } // namespace hykkt -} // namespace ReSolve \ No newline at end of file +} // namespace ReSolve diff --git a/resolve/hykkt/cholesky/CholeskySolverCuda.hpp b/resolve/hykkt/cholesky/CholeskySolverCuda.hpp index 0488ffcad..c5a4302e8 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCuda.hpp +++ b/resolve/hykkt/cholesky/CholeskySolverCuda.hpp @@ -6,9 +6,8 @@ #pragma once -#include - #include +#include #include #include #include @@ -46,7 +45,7 @@ namespace ReSolve // handle to the cuSPARSE library context cusolverSpHandle_t cusolverHandle_; - cusparseMatDescr_t descr_A_cusolver_; // descriptor for matrix A + cusparseMatDescr_t descr_A_cusolver_; // descriptor for matrix A csrcholInfo_t factorizationInfo_; // stores Cholesky factorization void* buffer_; // buffer for Cholesky factorization }; diff --git a/tests/functionality/testRefactor.cpp b/tests/functionality/testRefactor.cpp index 61fdbd119..76de7c518 100644 --- a/tests/functionality/testRefactor.cpp +++ b/tests/functionality/testRefactor.cpp @@ -27,9 +27,9 @@ #endif #ifdef RESOLVE_USE_CUDA +#include #include #include -#include #endif #include "TestHelper.hpp" diff --git a/tests/functionality/testSysRefactor.cpp b/tests/functionality/testSysRefactor.cpp index b775964e1..f203282ca 100644 --- a/tests/functionality/testSysRefactor.cpp +++ b/tests/functionality/testSysRefactor.cpp @@ -26,8 +26,8 @@ #include #ifdef RESOLVE_USE_CUDA -#include #include +#include #endif #ifdef RESOLVE_USE_HIP From fd1d82529edd118fa56648e79d893407156de0b4 Mon Sep 17 00:00:00 2001 From: Shaked Regev <35384901+shakedregev@users.noreply.github.com> Date: Wed, 22 Jul 2026 09:19:22 -0400 Subject: [PATCH 05/18] Apply suggestion from @shakedregev --- resolve/LinSolverDirectCuDssRf.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/LinSolverDirectCuDssRf.cpp b/resolve/LinSolverDirectCuDssRf.cpp index c55f32aee..ee88b08cf 100644 --- a/resolve/LinSolverDirectCuDssRf.cpp +++ b/resolve/LinSolverDirectCuDssRf.cpp @@ -32,7 +32,7 @@ namespace ReSolve * from the device memory. * * @pre The cuDssRf handle, config, and data have been created. - * @post The cuDssRf handle, config, and data have destroyed. + * @post The cuDssRf handle, config, and data have been destroyed. * * @pre The permutation vectors are allocated on the device. * @post The permutation vectors are deleted from the device. From c4c24bba45163d1268f099251c051b5e41c86b6a Mon Sep 17 00:00:00 2001 From: Shaked Regev <35384901+shakedregev@users.noreply.github.com> Date: Wed, 22 Jul 2026 09:24:04 -0400 Subject: [PATCH 06/18] Apply suggestion from @shakedregev --- resolve/LinSolverDirectCuDssRf.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/LinSolverDirectCuDssRf.cpp b/resolve/LinSolverDirectCuDssRf.cpp index ee88b08cf..666137791 100644 --- a/resolve/LinSolverDirectCuDssRf.cpp +++ b/resolve/LinSolverDirectCuDssRf.cpp @@ -148,7 +148,7 @@ namespace ReSolve { assert(A_ != nullptr && "Matrix A is null!"); assert(A_->getNumRows() > 0 && "Matrix A must have positive row count!"); - assert(A_->getNnz() > 0 && "Matrix A must have positive nonzero count!"); + assert(A_->getNnz() > 0 && "Matrix A must have positive number of nonzeros!"); return cudssExecute(handle_cudss_, CUDSS_PHASE_FACTORIZATION, From d7336cd5746fc84749370efaa6c00facfedf74fe Mon Sep 17 00:00:00 2001 From: Andrew Xu Date: Wed, 22 Jul 2026 14:44:41 +0000 Subject: [PATCH 07/18] Make cuDSS an optional dependency enabled via the RESOLVE_USE_CUDSS macro --- CMakeLists.txt | 13 ++++++++++++ cmake/ReSolveFindCudaLibraries.cmake | 7 +++++-- resolve/CMakeLists.txt | 10 ++++++++-- resolve/SystemSolver.cpp | 20 ++++++++++++++++--- resolve/hykkt/cholesky/CholeskySolver.cpp | 12 ++++++++++- resolve/hykkt/cholesky/CholeskySolver.hpp | 1 + resolve/hykkt/cholesky/CholeskySolverCpu.cpp | 3 ++- resolve/hykkt/cholesky/CholeskySolverCpu.hpp | 2 +- resolve/hykkt/cholesky/CholeskySolverCuda.cpp | 19 +++++++++++++++++- resolve/hykkt/cholesky/CholeskySolverCuda.hpp | 6 ++++++ resolve/hykkt/cholesky/CholeskySolverHip.cpp | 3 ++- resolve/hykkt/cholesky/CholeskySolverHip.hpp | 2 +- tests/functionality/testRefactor.cpp | 6 +++++- tests/functionality/testSysRefactor.cpp | 6 +++++- 14 files changed, 95 insertions(+), 15 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2873ea1c8..040ebf5bb 100644 --- a/CMakeLists.txt +++ b/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,18 @@ else() message(STATUS "Not using HIP") endif(RESOLVE_USE_HIP) +if(RESOLVE_USE_CUDSS) + find_package(cudss) + if(NOT cudss_FOUND) + 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 diff --git a/cmake/ReSolveFindCudaLibraries.cmake b/cmake/ReSolveFindCudaLibraries.cmake index be5137deb..e03ea264a 100644 --- a/cmake/ReSolveFindCudaLibraries.cmake +++ b/cmake/ReSolveFindCudaLibraries.cmake @@ -4,13 +4,16 @@ add_library(resolve_cuda INTERFACE) find_package(CUDAToolkit REQUIRED) -find_package(cudss REQUIRED) target_link_libraries( resolve_cuda INTERFACE CUDA::cusolver CUDA::cublas CUDA::cusparse - CUDA::cudart cudss + 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() diff --git a/resolve/CMakeLists.txt b/resolve/CMakeLists.txt index bf8c008c7..e64085db3 100644 --- a/resolve/CMakeLists.txt +++ b/resolve/CMakeLists.txt @@ -31,8 +31,11 @@ set(ReSolve_LUSOL_SRC LinSolverDirectLUSOL.cpp) # C++ code that links to CUDA SDK libraries set(ReSolve_CUDASDK_SRC LinSolverDirectCuSolverGLU.cpp LinSolverDirectCuSolverRf.cpp - LinSolverDirectCuDssRf.cpp LinSolverDirectCuSparseILU0.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 @@ -62,8 +65,11 @@ set(ReSolve_LUSOL_HEADER_INSTALL LinSolverDirectLUSOL.hpp) set(ReSolve_CUDA_HEADER_INSTALL LinSolverDirectCuSolverGLU.hpp LinSolverDirectCuSolverRf.hpp - LinSolverDirectCuDssRf.hpp LinSolverDirectCuSparseILU0.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 diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 1e5895677..8bf3ababd 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -16,11 +16,13 @@ #include #ifdef RESOLVE_USE_CUDA -#include #include #include #include #include +#ifdef RESOLVE_USE_CUDSS +#include +#endif #endif #ifdef RESOLVE_USE_HIP @@ -317,10 +319,12 @@ namespace ReSolve else if (refactorizationMethod_ == "cusolverrf") { refactorizationSolver_ = new ReSolve::LinSolverDirectCuSolverRf(); + #ifdef RESOLVE_USE_CUDSS } else if (refactorizationMethod_ == "cudssrf") { refactorizationSolver_ = new ReSolve::LinSolverDirectCuDssRf(); + #endif #endif #ifdef RESOLVE_USE_HIP } @@ -457,7 +461,11 @@ namespace ReSolve return factorizationSolver_->refactorize(); } - if (refactorizationMethod_ == "glu" || refactorizationMethod_ == "cusolverrf" || refactorizationMethod_ == "cudssrf" || refactorizationMethod_ == "rocsolverrf") + if (refactorizationMethod_ == "glu" || refactorizationMethod_ == "cusolverrf" || refactorizationMethod_ == "rocsolverrf" +#ifdef RESOLVE_USE_CUDSS + || refactorizationMethod_ == "cudssrf" +#endif + ) { is_solve_on_device_ = true; return refactorizationSolver_->refactorize(); @@ -513,6 +521,7 @@ namespace ReSolve is_solve_on_device_ = false; } + #ifdef RESOLVE_USE_CUDSS else if (refactorizationMethod_ == "cudssrf") { LinSolverDirectCuDssRf* Rf = dynamic_cast(refactorizationSolver_); @@ -522,6 +531,7 @@ namespace ReSolve is_solve_on_device_ = false; } + #endif #endif #ifdef RESOLVE_USE_HIP @@ -578,7 +588,11 @@ namespace ReSolve status += factorizationSolver_->solve(rhs, x); } - if (solveMethod_ == "glu" || solveMethod_ == "cusolverrf" || solveMethod_ == "cudssrf" || solveMethod_ == "rocsolverrf") + if (solveMethod_ == "glu" || solveMethod_ == "cusolverrf" || solveMethod_ == "rocsolverrf" +#ifdef RESOLVE_USE_CUDSS + || solveMethod_ == "cudssrf" +#endif + ) { if (is_solve_on_device_) { diff --git a/resolve/hykkt/cholesky/CholeskySolver.cpp b/resolve/hykkt/cholesky/CholeskySolver.cpp index 5063180cf..98603275a 100644 --- a/resolve/hykkt/cholesky/CholeskySolver.cpp +++ b/resolve/hykkt/cholesky/CholeskySolver.cpp @@ -25,6 +25,16 @@ namespace ReSolve * @param[in] memspace - memory space to use for computations */ CholeskySolver::CholeskySolver(memory::MemorySpace memspace) + : CholeskySolver(true, memspace) + { + } + + /** + * @brief Cholesky Solver constructor that allows specifying whether to use cuDSS in the CUDA implementation + * @param[in] memspace - memory space to use for computations + * @param[in] use_cudss - whether to use cuDSS (true) or cuSolver in CUDA + */ + CholeskySolver::CholeskySolver(bool use_cudss, memory::MemorySpace memspace) : memspace_(memspace) { if (memspace_ == memory::HOST) @@ -34,7 +44,7 @@ namespace ReSolve else { #ifdef RESOLVE_USE_CUDA - impl_ = new CholeskySolverCuda(); + impl_ = new CholeskySolverCuda(use_cudss); #elif defined(RESOLVE_USE_HIP) impl_ = new CholeskySolverHip(); #else diff --git a/resolve/hykkt/cholesky/CholeskySolver.hpp b/resolve/hykkt/cholesky/CholeskySolver.hpp index 0881ff56d..4be1aa3a6 100644 --- a/resolve/hykkt/cholesky/CholeskySolver.hpp +++ b/resolve/hykkt/cholesky/CholeskySolver.hpp @@ -18,6 +18,7 @@ namespace ReSolve { public: CholeskySolver(memory::MemorySpace memspace); + CholeskySolver(bool use_cudss, memory::MemorySpace memspace); ~CholeskySolver(); void addMatrixInfo(matrix::Csr* A); diff --git a/resolve/hykkt/cholesky/CholeskySolverCpu.cpp b/resolve/hykkt/cholesky/CholeskySolverCpu.cpp index 296585341..490a2e2bb 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCpu.cpp +++ b/resolve/hykkt/cholesky/CholeskySolverCpu.cpp @@ -13,7 +13,8 @@ namespace ReSolve namespace hykkt { - CholeskySolverCpu::CholeskySolverCpu() + // The bool argument is use_cudss, Used by CUDA version. It's defaulted to true and unused here. + CholeskySolverCpu::CholeskySolverCpu(bool) { Common_.nmethods = 1; // Use natural ordering diff --git a/resolve/hykkt/cholesky/CholeskySolverCpu.hpp b/resolve/hykkt/cholesky/CholeskySolverCpu.hpp index fcfda9e31..4b8731911 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCpu.hpp +++ b/resolve/hykkt/cholesky/CholeskySolverCpu.hpp @@ -15,7 +15,7 @@ namespace ReSolve class CholeskySolverCpu : public CholeskySolverImpl { public: - CholeskySolverCpu(); + CholeskySolverCpu(bool = true); // The bool argument is use_cudss, Used by CUDA version. It's defaulted to true and unused here. ~CholeskySolverCpu(); void addMatrixInfo(matrix::Csr* A); diff --git a/resolve/hykkt/cholesky/CholeskySolverCuda.cpp b/resolve/hykkt/cholesky/CholeskySolverCuda.cpp index 79b7f49c9..3b54e8f17 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCuda.cpp +++ b/resolve/hykkt/cholesky/CholeskySolverCuda.cpp @@ -14,8 +14,14 @@ namespace ReSolve namespace hykkt { CholeskySolverCuda::CholeskySolverCuda(bool use_cudss) - : use_cudss_(use_cudss) { +#ifdef RESOLVE_USE_CUDSS + use_cudss_ = use_cudss; +#else + use_cudss_ = false; +#endif + +#ifdef RESOLVE_USE_CUDSS if (use_cudss_) { cudssCreate(&cudss_handle_); @@ -23,6 +29,7 @@ namespace ReSolve cudssDataCreate(cudss_handle_, &cudss_data_); } else +#endif { cusolverSpCreate(&cusolverHandle_); cusparseCreateMatDescr(&descr_A_cusolver_); @@ -33,6 +40,7 @@ namespace ReSolve CholeskySolverCuda::~CholeskySolverCuda() { +#ifdef RESOLVE_USE_CUDSS if (use_cudss_) { cudssDataDestroy(cudss_handle_, cudss_data_); @@ -43,6 +51,7 @@ namespace ReSolve cudssMatrixDestroy(descr_x_); } else +#endif { cusolverSpDestroy(cusolverHandle_); cusparseDestroyMatDescr(descr_A_cusolver_); @@ -55,6 +64,7 @@ namespace ReSolve { A_ = A; +#ifdef RESOLVE_USE_CUDSS if (use_cudss_) { cudssMatrixCreateCsr(&descr_A_cudss_, @@ -71,6 +81,7 @@ namespace ReSolve CUDSS_MVIEW_LOWER, CUDSS_BASE_ZERO); } +#endif } /** @@ -78,6 +89,7 @@ namespace ReSolve */ void CholeskySolverCuda::symbolicAnalysis() { +#ifdef RESOLVE_USE_CUDSS if (use_cudss_) { cudssMatrixCreateDn(&descr_b_, @@ -103,6 +115,7 @@ namespace ReSolve descr_b_); } else +#endif { cusolverSpXcsrcholAnalysis(cusolverHandle_, A_->getNumRows(), @@ -139,6 +152,7 @@ namespace ReSolve */ void CholeskySolverCuda::numericalFactorization(real_type tol) { +#ifdef RESOLVE_USE_CUDSS if (use_cudss_) { cudssConfigSet(cudss_config_, CUDSS_CONFIG_PIVOT_EPSILON, &tol, sizeof(real_type)); @@ -155,6 +169,7 @@ namespace ReSolve } } else +#endif { int singularity = 0; cusolverSpDcsrcholFactor(cusolverHandle_, @@ -187,6 +202,7 @@ namespace ReSolve */ void CholeskySolverCuda::solve(vector::Vector* x, vector::Vector* b) { +#ifdef RESOLVE_USE_CUDSS if (use_cudss_) { if (descr_b_) @@ -220,6 +236,7 @@ namespace ReSolve } } else +#endif { for (index_type i = 0; i < b->getNumVectors(); i++) { diff --git a/resolve/hykkt/cholesky/CholeskySolverCuda.hpp b/resolve/hykkt/cholesky/CholeskySolverCuda.hpp index c5a4302e8..2d95237ed 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCuda.hpp +++ b/resolve/hykkt/cholesky/CholeskySolverCuda.hpp @@ -7,7 +7,11 @@ #pragma once #include + +#ifdef RESOLVE_USE_CUDSS #include +#endif + #include #include #include @@ -36,12 +40,14 @@ namespace ReSolve bool use_cudss_; +#ifdef RESOLVE_USE_CUDSS cudssHandle_t cudss_handle_; cudssConfig_t cudss_config_; cudssData_t cudss_data_; cudssMatrix_t descr_A_cudss_; cudssMatrix_t descr_b_; cudssMatrix_t descr_x_; +#endif // handle to the cuSPARSE library context cusolverSpHandle_t cusolverHandle_; diff --git a/resolve/hykkt/cholesky/CholeskySolverHip.cpp b/resolve/hykkt/cholesky/CholeskySolverHip.cpp index 358ebea9d..e4e0ee56b 100644 --- a/resolve/hykkt/cholesky/CholeskySolverHip.cpp +++ b/resolve/hykkt/cholesky/CholeskySolverHip.cpp @@ -13,7 +13,8 @@ namespace ReSolve namespace hykkt { - CholeskySolverHip::CholeskySolverHip() + // The bool argument is use_cudss, Used by CUDA version. It's defaulted to true and unused here. + CholeskySolverHip::CholeskySolverHip(bool) { rocblas_create_handle(&handle_); rocsolver_create_rfinfo(&rfinfo_, handle_); diff --git a/resolve/hykkt/cholesky/CholeskySolverHip.hpp b/resolve/hykkt/cholesky/CholeskySolverHip.hpp index 3661bcaa6..d84076b42 100644 --- a/resolve/hykkt/cholesky/CholeskySolverHip.hpp +++ b/resolve/hykkt/cholesky/CholeskySolverHip.hpp @@ -17,7 +17,7 @@ namespace ReSolve class CholeskySolverHip : public CholeskySolverImpl { public: - CholeskySolverHip(); + CholeskySolverHip(bool = true); // The bool argument is use_cudss, Used by CUDA version. It's defaulted to true and unused here. ~CholeskySolverHip(); void addMatrixInfo(matrix::Csr* A); diff --git a/tests/functionality/testRefactor.cpp b/tests/functionality/testRefactor.cpp index 76de7c518..7799000b8 100644 --- a/tests/functionality/testRefactor.cpp +++ b/tests/functionality/testRefactor.cpp @@ -27,9 +27,11 @@ #endif #ifdef RESOLVE_USE_CUDA -#include #include #include +#ifdef RESOLVE_USE_CUDSS +#include +#endif #endif #include "TestHelper.hpp" @@ -65,9 +67,11 @@ int main(int argc, char* argv[]) std::string cusolver_solver_name("cusolverRf"); error_sum += runTest(argc, argv, cusolver_solver_name); +#ifdef RESOLVE_USE_CUDSS std::string cudss_solver_name("cusdssRf"); error_sum += runTest(argc, argv, cudss_solver_name); +#endif } else { diff --git a/tests/functionality/testSysRefactor.cpp b/tests/functionality/testSysRefactor.cpp index f203282ca..6613736e9 100644 --- a/tests/functionality/testSysRefactor.cpp +++ b/tests/functionality/testSysRefactor.cpp @@ -26,8 +26,10 @@ #include #ifdef RESOLVE_USE_CUDA -#include #include +#ifdef RESOLVE_USE_CUDSS +#include +#endif #endif #ifdef RESOLVE_USE_HIP @@ -48,8 +50,10 @@ int main(int argc, char* argv[]) #ifdef RESOLVE_USE_CUDA error_sum += runTest(argc, argv, "cuda", false); +#ifdef RESOLVE_USE_CUDSS error_sum += runTest(argc, argv, "cuda", true); #endif +#endif #ifdef RESOLVE_USE_HIP error_sum += runTest(argc, argv, "hip", false); From 7dce176a84ef2de20ef4290871ce4a632c8966f1 Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Wed, 22 Jul 2026 14:45:26 +0000 Subject: [PATCH 08/18] Apply pre-commmit fixes --- CMakeLists.txt | 3 ++- resolve/SystemSolver.cpp | 12 ++++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 040ebf5bb..bc90a18d2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -189,7 +189,8 @@ if(RESOLVE_USE_CUDSS) if(NOT cudss_FOUND) message(STATUS "Cannot find cuDSS, disabling it ...") set(RESOLVE_USE_CUDSS - OFF CACHE BOOL FORCE + OFF + CACHE BOOL FORCE ) endif() else() diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 8bf3ababd..789285d9f 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -319,12 +319,12 @@ namespace ReSolve else if (refactorizationMethod_ == "cusolverrf") { refactorizationSolver_ = new ReSolve::LinSolverDirectCuSolverRf(); - #ifdef RESOLVE_USE_CUDSS +#ifdef RESOLVE_USE_CUDSS } else if (refactorizationMethod_ == "cudssrf") { refactorizationSolver_ = new ReSolve::LinSolverDirectCuDssRf(); - #endif +#endif #endif #ifdef RESOLVE_USE_HIP } @@ -463,7 +463,7 @@ namespace ReSolve if (refactorizationMethod_ == "glu" || refactorizationMethod_ == "cusolverrf" || refactorizationMethod_ == "rocsolverrf" #ifdef RESOLVE_USE_CUDSS - || refactorizationMethod_ == "cudssrf" + || refactorizationMethod_ == "cudssrf" #endif ) { @@ -521,7 +521,7 @@ namespace ReSolve is_solve_on_device_ = false; } - #ifdef RESOLVE_USE_CUDSS +#ifdef RESOLVE_USE_CUDSS else if (refactorizationMethod_ == "cudssrf") { LinSolverDirectCuDssRf* Rf = dynamic_cast(refactorizationSolver_); @@ -531,7 +531,7 @@ namespace ReSolve is_solve_on_device_ = false; } - #endif +#endif #endif #ifdef RESOLVE_USE_HIP @@ -590,7 +590,7 @@ namespace ReSolve if (solveMethod_ == "glu" || solveMethod_ == "cusolverrf" || solveMethod_ == "rocsolverrf" #ifdef RESOLVE_USE_CUDSS - || solveMethod_ == "cudssrf" + || solveMethod_ == "cudssrf" #endif ) { From 1b93dabb5e5af481c35f46ff3f6093ceb1cb20ab Mon Sep 17 00:00:00 2001 From: Andrew Xu Date: Wed, 22 Jul 2026 17:53:18 +0000 Subject: [PATCH 09/18] Add cudssRefactor example --- examples/CMakeLists.txt | 6 + examples/cudssRefactor.cpp | 302 ++++++++++++++++++++++++++++++++++++ resolve/resolve_defs.hpp.in | 1 + 3 files changed, 309 insertions(+) create mode 100644 examples/cudssRefactor.cpp diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index a8676ddb3..39c0b34ab 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -41,6 +41,12 @@ 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) diff --git a/examples/cudssRefactor.cpp b/examples/cudssRefactor.cpp new file mode 100644 index 000000000..4b70a1c49 --- /dev/null +++ b/examples/cudssRefactor.cpp @@ -0,0 +1,302 @@ +#include +#include +#include +#include +#include +#include + +#include "ExampleHelper.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +/// 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 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 -r -n \n\n"; + std::cout << "Optional features:\n"; + std::cout << "\t-e \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 +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(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 +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 (computing errors, printing summaries, etc.) + ExampleHelper 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(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; +} diff --git a/resolve/resolve_defs.hpp.in b/resolve/resolve_defs.hpp.in index 2d6b3b986..c4b409b04 100644 --- a/resolve/resolve_defs.hpp.in +++ b/resolve/resolve_defs.hpp.in @@ -12,6 +12,7 @@ #cmakedefine RESOLVE_USE_RAJA #cmakedefine RESOLVE_USE_EIGEN #cmakedefine RESOLVE_USE_KLU +#cmakedefine RESOLVE_USE_CUDSS #cmakedefine RESOLVE_USE_PROFILING #define RESOLVE_VERSION "@PROJECT_VERSION@" From a39a6897c9ae05774591b4f5915ea2e637b45aaf Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Wed, 22 Jul 2026 17:53:38 +0000 Subject: [PATCH 10/18] Apply pre-commmit fixes --- examples/CMakeLists.txt | 3 ++- examples/cudssRefactor.cpp | 9 ++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 39c0b34ab..508b9a944 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -42,7 +42,8 @@ if(RESOLVE_USE_KLU) target_link_libraries(gluRefactor.exe PRIVATE ReSolve) if(RESOLVE_USE_CUDSS) - # Build an example with a configurable and portable system solver implemented with 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) diff --git a/examples/cudssRefactor.cpp b/examples/cudssRefactor.cpp index 4b70a1c49..c57dd5480 100644 --- a/examples/cudssRefactor.cpp +++ b/examples/cudssRefactor.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -17,8 +18,6 @@ #include #include -#include - /// Prints help message describing system usage. void printHelpInfo() { @@ -160,11 +159,11 @@ int cudssRefactor(int argc, char* argv[]) // Create system solver ReSolve::SystemSolver solver(&workspace, - "klu", // factorization + "klu", // factorization "cudssrf", // refactorization "cudssrf", // triangular solve - "none", // preconditioner (always 'none' here) - "none"); // iterative refinement + "none", // preconditioner (always 'none' here) + "none"); // iterative refinement if (is_iterative_refinement) { From b024bca410c76dbe7ad0c888aaf42f96b3457a7b Mon Sep 17 00:00:00 2001 From: Andrew Xu Date: Thu, 23 Jul 2026 15:33:35 +0000 Subject: [PATCH 11/18] Make Cholesky unit tests test both the CuSolver and the cuDSS implementations --- tests/unit/hykkt/HykktCholeskyTests.hpp | 7 ++++--- tests/unit/hykkt/runHykktCholeskyTests.cpp | 23 ++++++++++++++++------ 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/tests/unit/hykkt/HykktCholeskyTests.hpp b/tests/unit/hykkt/HykktCholeskyTests.hpp index e04d39aeb..ba604869b 100644 --- a/tests/unit/hykkt/HykktCholeskyTests.hpp +++ b/tests/unit/hykkt/HykktCholeskyTests.hpp @@ -27,8 +27,8 @@ namespace ReSolve class HykktCholeskyTests : public TestBase { public: - HykktCholeskyTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler, std::mt19937& generator) - : memspace_(memspace), matrixHandler_(matrixHandler), generator_(generator) + HykktCholeskyTests(bool use_cudss, memory::MemorySpace memspace, MatrixHandler& matrixHandler, std::mt19937& generator) + : use_cudss_(use_cudss), memspace_(memspace), matrixHandler_(matrixHandler), generator_(generator) { cholmod_start(&Common); } @@ -60,7 +60,7 @@ namespace ReSolve A->syncData(memory::DEVICE); } - ReSolve::hykkt::CholeskySolver solver(memspace_); + ReSolve::hykkt::CholeskySolver solver(use_cudss_, memspace_); solver.addMatrixInfo(A); solver.symbolicAnalysis(); solver.setPivotTolerance(1e-12); @@ -273,6 +273,7 @@ namespace ReSolve ReSolve::memory::MemorySpace memspace_; MatrixHandler& matrixHandler_; std::mt19937& generator_; + bool use_cudss_; cholmod_common Common; diff --git a/tests/unit/hykkt/runHykktCholeskyTests.cpp b/tests/unit/hykkt/runHykktCholeskyTests.cpp index 55ebeb7d6..dfe4dc1be 100644 --- a/tests/unit/hykkt/runHykktCholeskyTests.cpp +++ b/tests/unit/hykkt/runHykktCholeskyTests.cpp @@ -20,15 +20,23 @@ * @param result - test results */ template -void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace, ReSolve::tests::TestingResults& result) +void runTests(bool use_cudss, const std::string& backend, ReSolve::memory::MemorySpace memspace, ReSolve::tests::TestingResults& result) { - std::cout << "Running tests on " << backend << " device:\n"; + if (backend == "CUDA") + { + std::string impl = use_cudss ? "cuDSS" : "CuSolver"; + std::cout << "Running tests on " << backend << " device with " << impl << " implementation:\n"; + } + else + { + std::cout << "Running tests on " << backend << " device:\n"; + } WorkspaceType workspace; workspace.initializeHandles(); ReSolve::MatrixHandler handler(&workspace); std::mt19937 generator(ReSolve::constants::SEED); // set random seed for reproducibility - ReSolve::tests::HykktCholeskyTests test(memspace, handler, generator); + ReSolve::tests::HykktCholeskyTests test(use_cudss, memspace, handler, generator); result += test.minimalCorrectness(); handler.setValuesChanged(true, memspace); @@ -49,14 +57,17 @@ void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace, int main(int, char**) { ReSolve::tests::TestingResults result; - runTests("CPU", ReSolve::memory::HOST, result); + runTests(false, "CPU", ReSolve::memory::HOST, result); #ifdef RESOLVE_USE_CUDA - runTests("CUDA", ReSolve::memory::DEVICE, result); + runTests(false, "CUDA", ReSolve::memory::DEVICE, result); +#ifdef RESOLVE_USE_CUDSS +#endif + runTests(true, "CUDA", ReSolve::memory::DEVICE, result); #endif #ifdef RESOLVE_USE_HIP - runTests("HIP", ReSolve::memory::DEVICE, result); + runTests(false, "HIP", ReSolve::memory::DEVICE, result); #endif return result.summary(); From df242ed972a77490f3a23705aac5f7a20a9b3621 Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Thu, 23 Jul 2026 15:33:55 +0000 Subject: [PATCH 12/18] Apply pre-commmit fixes --- tests/unit/hykkt/HykktCholeskyTests.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/hykkt/HykktCholeskyTests.hpp b/tests/unit/hykkt/HykktCholeskyTests.hpp index ba604869b..99eca2618 100644 --- a/tests/unit/hykkt/HykktCholeskyTests.hpp +++ b/tests/unit/hykkt/HykktCholeskyTests.hpp @@ -273,7 +273,7 @@ namespace ReSolve ReSolve::memory::MemorySpace memspace_; MatrixHandler& matrixHandler_; std::mt19937& generator_; - bool use_cudss_; + bool use_cudss_; cholmod_common Common; From f07044b9d5c4f2b54cd9a1a0adcb829dfc270ecc Mon Sep 17 00:00:00 2001 From: Andrew Xu Date: Fri, 24 Jul 2026 16:16:43 +0000 Subject: [PATCH 13/18] Split CholeskySolverCuda into two versions and reverted changes to the other Cholesky classes --- CMakeLists.txt | 2 +- resolve/hykkt/CMakeLists.txt | 4 +- resolve/hykkt/cholesky/CholeskySolver.cpp | 12 +- resolve/hykkt/cholesky/CholeskySolverCpu.cpp | 3 +- resolve/hykkt/cholesky/CholeskySolverCpu.hpp | 2 +- resolve/hykkt/cholesky/CholeskySolverCuda.cpp | 246 +++---------- resolve/hykkt/cholesky/CholeskySolverCuda.hpp | 24 +- resolve/hykkt/cholesky/CholeskySolverHip.cpp | 3 +- resolve/hykkt/cholesky/CholeskySolverHip.hpp | 2 +- resolve/hykkt/cholesky_cudss/CMakeLists.txt | 35 ++ .../cholesky_cudss/CholeskySolverCuDss.cpp | 92 +++++ .../cholesky_cudss/CholeskySolverCuDss.hpp | 37 ++ .../CholeskySolverCuDssCuda.cpp | 142 ++++++++ .../CholeskySolverCuDssCuda.hpp | 44 +++ tests/functionality/testSysRefactor.cpp | 9 +- tests/unit/hykkt/CMakeLists.txt | 14 + tests/unit/hykkt/HykktCholeskyCuDssTests.hpp | 337 ++++++++++++++++++ tests/unit/hykkt/HykktCholeskyTests.hpp | 7 +- .../unit/hykkt/runHykktCholeskyCuDssTests.cpp | 55 +++ tests/unit/hykkt/runHykktCholeskyTests.cpp | 23 +- 20 files changed, 833 insertions(+), 260 deletions(-) create mode 100644 resolve/hykkt/cholesky_cudss/CMakeLists.txt create mode 100644 resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp create mode 100644 resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.hpp create mode 100644 resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.cpp create mode 100644 resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.hpp create mode 100644 tests/unit/hykkt/HykktCholeskyCuDssTests.hpp create mode 100644 tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index bc90a18d2..c366ee3c1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -186,7 +186,7 @@ endif(RESOLVE_USE_HIP) if(RESOLVE_USE_CUDSS) find_package(cudss) - if(NOT cudss_FOUND) + if(NOT (cudss_FOUND AND RESOLVE_USE_CUDA)) message(STATUS "Cannot find cuDSS, disabling it ...") set(RESOLVE_USE_CUDSS OFF diff --git a/resolve/hykkt/CMakeLists.txt b/resolve/hykkt/CMakeLists.txt index e9adad2d0..265d579b3 100644 --- a/resolve/hykkt/CMakeLists.txt +++ b/resolve/hykkt/CMakeLists.txt @@ -11,6 +11,9 @@ add_subdirectory(ruiz) add_subdirectory(cholesky) add_subdirectory(spgemm) add_subdirectory(sccg) +if(RESOLVE_USE_CUDSS) + add_subdirectory(cholesky_cudss) +endif() # Build shared library ReSolve::hykkt set(HYKKT_SOLVER_SRC HyKKTSolver.cpp) @@ -30,7 +33,6 @@ target_link_libraries( resolve_vector resolve_matrix ) - # Link to CUDA ReSolve backend if CUDA is support enabled. # target_include_directories(resolve_hykkt_solver PUBLIC diff --git a/resolve/hykkt/cholesky/CholeskySolver.cpp b/resolve/hykkt/cholesky/CholeskySolver.cpp index 98603275a..5063180cf 100644 --- a/resolve/hykkt/cholesky/CholeskySolver.cpp +++ b/resolve/hykkt/cholesky/CholeskySolver.cpp @@ -25,16 +25,6 @@ namespace ReSolve * @param[in] memspace - memory space to use for computations */ CholeskySolver::CholeskySolver(memory::MemorySpace memspace) - : CholeskySolver(true, memspace) - { - } - - /** - * @brief Cholesky Solver constructor that allows specifying whether to use cuDSS in the CUDA implementation - * @param[in] memspace - memory space to use for computations - * @param[in] use_cudss - whether to use cuDSS (true) or cuSolver in CUDA - */ - CholeskySolver::CholeskySolver(bool use_cudss, memory::MemorySpace memspace) : memspace_(memspace) { if (memspace_ == memory::HOST) @@ -44,7 +34,7 @@ namespace ReSolve else { #ifdef RESOLVE_USE_CUDA - impl_ = new CholeskySolverCuda(use_cudss); + impl_ = new CholeskySolverCuda(); #elif defined(RESOLVE_USE_HIP) impl_ = new CholeskySolverHip(); #else diff --git a/resolve/hykkt/cholesky/CholeskySolverCpu.cpp b/resolve/hykkt/cholesky/CholeskySolverCpu.cpp index 490a2e2bb..296585341 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCpu.cpp +++ b/resolve/hykkt/cholesky/CholeskySolverCpu.cpp @@ -13,8 +13,7 @@ namespace ReSolve namespace hykkt { - // The bool argument is use_cudss, Used by CUDA version. It's defaulted to true and unused here. - CholeskySolverCpu::CholeskySolverCpu(bool) + CholeskySolverCpu::CholeskySolverCpu() { Common_.nmethods = 1; // Use natural ordering diff --git a/resolve/hykkt/cholesky/CholeskySolverCpu.hpp b/resolve/hykkt/cholesky/CholeskySolverCpu.hpp index 4b8731911..fcfda9e31 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCpu.hpp +++ b/resolve/hykkt/cholesky/CholeskySolverCpu.hpp @@ -15,7 +15,7 @@ namespace ReSolve class CholeskySolverCpu : public CholeskySolverImpl { public: - CholeskySolverCpu(bool = true); // The bool argument is use_cudss, Used by CUDA version. It's defaulted to true and unused here. + CholeskySolverCpu(); ~CholeskySolverCpu(); void addMatrixInfo(matrix::Csr* A); diff --git a/resolve/hykkt/cholesky/CholeskySolverCuda.cpp b/resolve/hykkt/cholesky/CholeskySolverCuda.cpp index 3b54e8f17..13cf29514 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCuda.cpp +++ b/resolve/hykkt/cholesky/CholeskySolverCuda.cpp @@ -13,182 +13,87 @@ namespace ReSolve namespace hykkt { - CholeskySolverCuda::CholeskySolverCuda(bool use_cudss) + CholeskySolverCuda::CholeskySolverCuda() { -#ifdef RESOLVE_USE_CUDSS - use_cudss_ = use_cudss; -#else - use_cudss_ = false; -#endif - -#ifdef RESOLVE_USE_CUDSS - if (use_cudss_) - { - cudssCreate(&cudss_handle_); - cudssConfigCreate(&cudss_config_); - cudssDataCreate(cudss_handle_, &cudss_data_); - } - else -#endif - { - cusolverSpCreate(&cusolverHandle_); - cusparseCreateMatDescr(&descr_A_cusolver_); - cusolverSpCreateCsrcholInfo(&factorizationInfo_); - buffer_ = nullptr; - } + cusolverSpCreate(&cusolverHandle_); + cusparseCreateMatDescr(&descrA_); + cusolverSpCreateCsrcholInfo(&factorizationInfo_); + buffer_ = nullptr; } CholeskySolverCuda::~CholeskySolverCuda() { -#ifdef RESOLVE_USE_CUDSS - if (use_cudss_) - { - cudssDataDestroy(cudss_handle_, cudss_data_); - cudssConfigDestroy(cudss_config_); - cudssDestroy(cudss_handle_); - cudssMatrixDestroy(descr_A_cudss_); - cudssMatrixDestroy(descr_b_); - cudssMatrixDestroy(descr_x_); - } - else -#endif - { - cusolverSpDestroy(cusolverHandle_); - cusparseDestroyMatDescr(descr_A_cusolver_); - cusolverSpDestroyCsrcholInfo(factorizationInfo_); - mem_.deleteOnDevice(buffer_); - } + cusolverSpDestroy(cusolverHandle_); + cusparseDestroyMatDescr(descrA_); + cusolverSpDestroyCsrcholInfo(factorizationInfo_); + mem_.deleteOnDevice(buffer_); } void CholeskySolverCuda::addMatrixInfo(matrix::Csr* A) { A_ = A; - -#ifdef RESOLVE_USE_CUDSS - if (use_cudss_) - { - cudssMatrixCreateCsr(&descr_A_cudss_, - A_->getNumRows(), - A_->getNumColumns(), - A_->getNnz(), - A_->getRowData(memory::DEVICE), - nullptr, // Row end offsets (null for standard CSR) - A_->getColData(memory::DEVICE), - A_->getValues(memory::DEVICE), - CUDA_R_32I, - CUDA_R_64F, - CUDSS_MTYPE_SPD, - CUDSS_MVIEW_LOWER, - CUDSS_BASE_ZERO); - } -#endif } /** * @brief Perform symbolic analysis for the Cholesky factorization + * + * Uses the `cusolverSpXcsrcholAnalysis` routine. */ void CholeskySolverCuda::symbolicAnalysis() { -#ifdef RESOLVE_USE_CUDSS - if (use_cudss_) - { - cudssMatrixCreateDn(&descr_b_, - A_->getNumRows(), - 1, - A_->getNumRows(), - nullptr, - CUDA_R_64F, - CUDSS_LAYOUT_COL_MAJOR); - cudssMatrixCreateDn(&descr_x_, - A_->getNumRows(), - 1, - A_->getNumRows(), - nullptr, - CUDA_R_64F, - CUDSS_LAYOUT_COL_MAJOR); - cudssExecute(cudss_handle_, - CUDSS_PHASE_ANALYSIS, - cudss_config_, - cudss_data_, - descr_A_cudss_, - descr_x_, - descr_b_); - } - else -#endif - { - cusolverSpXcsrcholAnalysis(cusolverHandle_, + cusolverSpXcsrcholAnalysis(cusolverHandle_, + A_->getNumRows(), + A_->getNnz(), + descrA_, + A_->getRowData(memory::DEVICE), + A_->getColData(memory::DEVICE), + factorizationInfo_); + // Calculate size of buffer needed + size_t internalDataBytes = 0; + size_t workspaceBytes = 0; + cusolverSpDcsrcholBufferInfo(cusolverHandle_, A_->getNumRows(), A_->getNnz(), - descr_A_cusolver_, + descrA_, + A_->getValues(memory::DEVICE), A_->getRowData(memory::DEVICE), A_->getColData(memory::DEVICE), - factorizationInfo_); - // Calculate size of buffer needed - size_t internalDataBytes = 0; - size_t workspaceBytes = 0; - cusolverSpDcsrcholBufferInfo(cusolverHandle_, - A_->getNumRows(), - A_->getNnz(), - descr_A_cusolver_, - A_->getValues(memory::DEVICE), - A_->getRowData(memory::DEVICE), - A_->getColData(memory::DEVICE), - factorizationInfo_, - &internalDataBytes, - &workspaceBytes); - if (buffer_ != nullptr) - { - mem_.deleteOnDevice(buffer_); - } - mem_.allocateBufferOnDevice(&buffer_, workspaceBytes); + factorizationInfo_, + &internalDataBytes, + &workspaceBytes); + if (buffer_ != nullptr) + { + mem_.deleteOnDevice(buffer_); } + mem_.allocateBufferOnDevice(&buffer_, workspaceBytes); } /** * @brief Perform numerical factorization for the Cholesky factorization * + * Uses the `cusolverSpDcsrcholFactor` routine. + * * @param[in] tol - Tolerance for zero pivot detection. */ void CholeskySolverCuda::numericalFactorization(real_type tol) { -#ifdef RESOLVE_USE_CUDSS - if (use_cudss_) + int singularity = 0; + cusolverSpDcsrcholFactor(cusolverHandle_, + A_->getNumRows(), + A_->getNnz(), + descrA_, + A_->getValues(memory::DEVICE), + A_->getRowData(memory::DEVICE), + A_->getColData(memory::DEVICE), + factorizationInfo_, + buffer_); + cusolverSpDcsrcholZeroPivot(cusolverHandle_, + factorizationInfo_, + tol, + &singularity); + if (singularity >= 0) { - cudssConfigSet(cudss_config_, CUDSS_CONFIG_PIVOT_EPSILON, &tol, sizeof(real_type)); - cudssStatus_t status = cudssExecute(cudss_handle_, - CUDSS_PHASE_FACTORIZATION, - cudss_config_, - cudss_data_, - descr_A_cudss_, - descr_x_, - descr_b_); - if (status != CUDSS_STATUS_SUCCESS) - { - out::error() << "Cholesky factorization failed with status: " << status << "\n"; - } - } - else -#endif - { - int singularity = 0; - cusolverSpDcsrcholFactor(cusolverHandle_, - A_->getNumRows(), - A_->getNnz(), - descr_A_cusolver_, - A_->getValues(memory::DEVICE), - A_->getRowData(memory::DEVICE), - A_->getColData(memory::DEVICE), - factorizationInfo_, - buffer_); - cusolverSpDcsrcholZeroPivot(cusolverHandle_, - factorizationInfo_, - tol, - &singularity); - if (singularity >= 0) - { - out::error() << "Cholesky factorization failed with singularity at index: " << singularity << "\n"; - } + out::error() << "Cholesky factorization failed with singularity at index: " << singularity << "\n"; } } @@ -202,53 +107,12 @@ namespace ReSolve */ void CholeskySolverCuda::solve(vector::Vector* x, vector::Vector* b) { -#ifdef RESOLVE_USE_CUDSS - if (use_cudss_) - { - if (descr_b_) - { - cudssMatrixDestroy(descr_b_); - } - if (descr_x_) - { - cudssMatrixDestroy(descr_x_); - } - - cudssMatrixCreateDn(&descr_b_, - b->getSize(), - b->getNumVectors(), - b->getSize(), - b->getData(memory::DEVICE), - CUDA_R_64F, - CUDSS_LAYOUT_COL_MAJOR); - cudssMatrixCreateDn(&descr_x_, - b->getSize(), - b->getNumVectors(), - b->getSize(), - x->getData(memory::DEVICE), - CUDA_R_64F, - CUDSS_LAYOUT_COL_MAJOR); - - cudssStatus_t status = cudssExecute(cudss_handle_, CUDSS_PHASE_SOLVE, cudss_config_, cudss_data_, descr_A_cudss_, descr_x_, descr_b_); - if (status != CUDSS_STATUS_SUCCESS) - { - out::error() << "cuDSS triangular solve failed with status: " << status << "\n"; - } - } - else -#endif - { - for (index_type i = 0; i < b->getNumVectors(); i++) - { - cusolverSpDcsrcholSolve(cusolverHandle_, - A_->getNumRows(), - b->getData(i, memory::DEVICE), - x->getData(i, memory::DEVICE), - factorizationInfo_, - buffer_); - } - } - + cusolverSpDcsrcholSolve(cusolverHandle_, + A_->getNumRows(), + b->getData(memory::DEVICE), + x->getData(memory::DEVICE), + factorizationInfo_, + buffer_); x->setDataUpdated(memory::DEVICE); } } // namespace hykkt diff --git a/resolve/hykkt/cholesky/CholeskySolverCuda.hpp b/resolve/hykkt/cholesky/CholeskySolverCuda.hpp index 2d95237ed..29cce95de 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCuda.hpp +++ b/resolve/hykkt/cholesky/CholeskySolverCuda.hpp @@ -1,17 +1,10 @@ /** * @file CholeskySolverCuda.hpp * @author Adham Ibrahim (ibrahimas@ornl.gov) - * @brief Header for CUDA implementation of Cholesky Solver using cuDSS + * @brief Header for CUDA implementation of Cholesky Solver */ -#pragma once - #include - -#ifdef RESOLVE_USE_CUDSS -#include -#endif - #include #include #include @@ -25,7 +18,7 @@ namespace ReSolve class CholeskySolverCuda : public CholeskySolverImpl { public: - CholeskySolverCuda(bool use_cudss = true); + CholeskySolverCuda(); ~CholeskySolverCuda(); void addMatrixInfo(matrix::Csr* A); @@ -38,20 +31,9 @@ namespace ReSolve matrix::Csr* A_; // pointer to the input matrix - bool use_cudss_; - -#ifdef RESOLVE_USE_CUDSS - cudssHandle_t cudss_handle_; - cudssConfig_t cudss_config_; - cudssData_t cudss_data_; - cudssMatrix_t descr_A_cudss_; - cudssMatrix_t descr_b_; - cudssMatrix_t descr_x_; -#endif - // handle to the cuSPARSE library context cusolverSpHandle_t cusolverHandle_; - cusparseMatDescr_t descr_A_cusolver_; // descriptor for matrix A + cusparseMatDescr_t descrA_; // descriptor for matrix A csrcholInfo_t factorizationInfo_; // stores Cholesky factorization void* buffer_; // buffer for Cholesky factorization }; diff --git a/resolve/hykkt/cholesky/CholeskySolverHip.cpp b/resolve/hykkt/cholesky/CholeskySolverHip.cpp index e4e0ee56b..358ebea9d 100644 --- a/resolve/hykkt/cholesky/CholeskySolverHip.cpp +++ b/resolve/hykkt/cholesky/CholeskySolverHip.cpp @@ -13,8 +13,7 @@ namespace ReSolve namespace hykkt { - // The bool argument is use_cudss, Used by CUDA version. It's defaulted to true and unused here. - CholeskySolverHip::CholeskySolverHip(bool) + CholeskySolverHip::CholeskySolverHip() { rocblas_create_handle(&handle_); rocsolver_create_rfinfo(&rfinfo_, handle_); diff --git a/resolve/hykkt/cholesky/CholeskySolverHip.hpp b/resolve/hykkt/cholesky/CholeskySolverHip.hpp index d84076b42..3661bcaa6 100644 --- a/resolve/hykkt/cholesky/CholeskySolverHip.hpp +++ b/resolve/hykkt/cholesky/CholeskySolverHip.hpp @@ -17,7 +17,7 @@ namespace ReSolve class CholeskySolverHip : public CholeskySolverImpl { public: - CholeskySolverHip(bool = true); // The bool argument is use_cudss, Used by CUDA version. It's defaulted to true and unused here. + CholeskySolverHip(); ~CholeskySolverHip(); void addMatrixInfo(matrix::Csr* A); diff --git a/resolve/hykkt/cholesky_cudss/CMakeLists.txt b/resolve/hykkt/cholesky_cudss/CMakeLists.txt new file mode 100644 index 000000000..ef26c395e --- /dev/null +++ b/resolve/hykkt/cholesky_cudss/CMakeLists.txt @@ -0,0 +1,35 @@ +#[[ + +@brief Build ReSolve matrix module + +@author Slaven Peles + +]] + +set(HyKKT_CHOL_CUDSS_SRC CholeskySolverCuDss.cpp CholeskySolverCuDssCuda.cpp) + +# Header files to be installed +set(HyKKT_CHOL_CUDSS_HEADER_INSTALL + CholeskySolverCuDss.hpp CholeskySolverCuDssCuda.hpp +) + +add_library(resolve_hykkt_chol_cudss SHARED ${HyKKT_CHOL_CUDSS_SRC}) + +target_link_libraries(resolve_hykkt_chol_cudss PUBLIC ${suitesparse_cholmod}) +target_include_directories(resolve_hykkt_chol_cudss PUBLIC ${SUITESPARSE_INCLUDE_DIR}) + +# Link to CUDA ReSolve backend if CUDA is support enabled +target_sources(resolve_hykkt_chol_cudss PRIVATE ${HyKKT_CHOL_CUDSS_SRC}) +target_link_libraries(resolve_hykkt_chol_cudss PUBLIC resolve_backend_cuda) + +target_link_libraries( + resolve_hykkt_chol_cudss PUBLIC resolve_workspace resolve_vector resolve_matrix + resolve_logger +) + +target_include_directories( + resolve_hykkt_chol_cudss INTERFACE $ + $ +) + +install(FILES ${HyKKT_CHOL_CUDSS_HEADER_INSTALL} DESTINATION include/resolve/hykkt) diff --git a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp new file mode 100644 index 000000000..e0b2ae632 --- /dev/null +++ b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp @@ -0,0 +1,92 @@ +/** + * @file CholeskySolverCuDss.cpp + * @author Adham Ibrahim (ibrahimas@ornl.gov) + * @brief Cholesky decomposition solver CuDSS implementation. This is a CUDA-only variant of CholeskySolver. + */ + +#include "CholeskySolverCuDss.hpp" + +#include "CholeskySolverCuDssCuda.hpp" + +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 + */ + CholeskySolverCuDss::CholeskySolverCuDss(memory::MemorySpace memspace) + : memspace_(memspace), + impl_(new CholeskySolverCuDssCuda()) + { + } + + /** + * @brief Cholesky Solver destructor + */ + CholeskySolverCuDss::~CholeskySolverCuDss() + { + delete impl_; + } + + /** + * @brief Loads or reloads matrix pointer to the solver + * @param[in] A - pointer to the matrix in CSR format + */ + void CholeskySolverCuDss::addMatrixInfo(matrix::Csr* A) + { + A_ = A; + impl_->addMatrixInfo(A); + } + + /** + * @brief Performs symbolic analysis. Determines the sparsity pattern of + * the factor L. Values will be computed by numerical analysis. + * This need only be called once as long as the sparsity pattern does not change. + */ + void CholeskySolverCuDss::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 CholeskySolverCuDss::setPivotTolerance(real_type tol) + { + tol_ = tol; + } + + /** + * @brief Performs numerical factorization. Fills in the values of the factor L such + * that LL^T = A. + */ + void CholeskySolverCuDss::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 CholeskySolverCuDss::solve(vector::Vector* x, vector::Vector* b) + { + impl_->solve(x, b); + x->setDataUpdated(memspace_); + } + } // namespace hykkt +} // namespace ReSolve diff --git a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.hpp b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.hpp new file mode 100644 index 000000000..4e3a11110 --- /dev/null +++ b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.hpp @@ -0,0 +1,37 @@ +/** + * @file CholeskySolverCuDss.hpp + * @author Adham Ibrahim (ibrahimas@ornl.gov) + * @brief Cholesky decomposition (cuDSS implementation) solver header. This is a CUDA-only variant of CholeskySolver. + */ + +#pragma once +#include "CholeskySolverCuDssCuda.hpp" +#include +#include +#include + +namespace ReSolve +{ + namespace hykkt + { + class CholeskySolverCuDss + { + public: + CholeskySolverCuDss(memory::MemorySpace memspace); + ~CholeskySolverCuDss(); + + 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; + CholeskySolverCuDssCuda* impl_; + }; + } // namespace hykkt +} // namespace ReSolve diff --git a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.cpp b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.cpp new file mode 100644 index 000000000..fa52c6eb0 --- /dev/null +++ b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.cpp @@ -0,0 +1,142 @@ +/** + * @file CholeskySolverCuDssCuda.cpp + * @author Andrew Xu (xua1@ornl.gov) + * @brief CUDA cuDSS implementation of Cholesky Solver + */ + +#include "CholeskySolverCuDssCuda.hpp" + +namespace ReSolve +{ + using real_type = ReSolve::real_type; + using out = ReSolve::io::Logger; + + namespace hykkt + { + CholeskySolverCuDssCuda::CholeskySolverCuDssCuda() + { + cudssCreate(&cudss_handle_); + cudssConfigCreate(&cudss_config_); + cudssDataCreate(cudss_handle_, &cudss_data_); + } + + CholeskySolverCuDssCuda::~CholeskySolverCuDssCuda() + { + cudssDataDestroy(cudss_handle_, cudss_data_); + cudssConfigDestroy(cudss_config_); + cudssDestroy(cudss_handle_); + cudssMatrixDestroy(descr_A_cudss_); + cudssMatrixDestroy(descr_b_); + cudssMatrixDestroy(descr_x_); + } + + void CholeskySolverCuDssCuda::addMatrixInfo(matrix::Csr* A) + { + A_ = A; + + cudssMatrixCreateCsr(&descr_A_cudss_, + A_->getNumRows(), + A_->getNumColumns(), + A_->getNnz(), + A_->getRowData(memory::DEVICE), + nullptr, // Row end offsets (null for standard CSR) + A_->getColData(memory::DEVICE), + A_->getValues(memory::DEVICE), + CUDA_R_32I, + CUDA_R_64F, + CUDSS_MTYPE_SPD, + CUDSS_MVIEW_LOWER, + CUDSS_BASE_ZERO); + } + + /** + * @brief Perform symbolic analysis for the Cholesky factorization + */ + void CholeskySolverCuDssCuda::symbolicAnalysis() + { + cudssMatrixCreateDn(&descr_b_, + A_->getNumRows(), + 1, + A_->getNumRows(), + nullptr, + CUDA_R_64F, + CUDSS_LAYOUT_COL_MAJOR); + cudssMatrixCreateDn(&descr_x_, + A_->getNumRows(), + 1, + A_->getNumRows(), + nullptr, + CUDA_R_64F, + CUDSS_LAYOUT_COL_MAJOR); + cudssExecute(cudss_handle_, + CUDSS_PHASE_ANALYSIS, + cudss_config_, + cudss_data_, + descr_A_cudss_, + descr_x_, + descr_b_); + } + + /** + * @brief Perform numerical factorization for the Cholesky factorization + * + * @param[in] tol - Tolerance for zero pivot detection. + */ + void CholeskySolverCuDssCuda::numericalFactorization(real_type tol) + { + cudssConfigSet(cudss_config_, CUDSS_CONFIG_PIVOT_EPSILON, &tol, sizeof(real_type)); + cudssStatus_t status = cudssExecute(cudss_handle_, + CUDSS_PHASE_FACTORIZATION, + cudss_config_, + cudss_data_, + descr_A_cudss_, + descr_x_, + descr_b_); + if (status != CUDSS_STATUS_SUCCESS) + { + out::error() << "Cholesky factorization failed with status: " << status << "\n"; + } + } + + /** + * @brief Solve the linear system Ax = b + * + * Uses the `cusolverSpDcsrcholSolve` routine. + * + * @param[out] x - Solution vector. + * @param[in] b - Right-hand side vector. + */ + void CholeskySolverCuDssCuda::solve(vector::Vector* x, vector::Vector* b) + { + if (descr_b_) + { + cudssMatrixDestroy(descr_b_); + } + if (descr_x_) + { + cudssMatrixDestroy(descr_x_); + } + + cudssMatrixCreateDn(&descr_b_, + b->getSize(), + b->getNumVectors(), + b->getSize(), + b->getData(memory::DEVICE), + CUDA_R_64F, + CUDSS_LAYOUT_COL_MAJOR); + cudssMatrixCreateDn(&descr_x_, + b->getSize(), + b->getNumVectors(), + b->getSize(), + x->getData(memory::DEVICE), + CUDA_R_64F, + CUDSS_LAYOUT_COL_MAJOR); + + cudssStatus_t status = cudssExecute(cudss_handle_, CUDSS_PHASE_SOLVE, cudss_config_, cudss_data_, descr_A_cudss_, descr_x_, descr_b_); + if (status != CUDSS_STATUS_SUCCESS) + { + out::error() << "cuDSS triangular solve failed with status: " << status << "\n"; + } + } + } // namespace hykkt +} // namespace ReSolve diff --git a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.hpp b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.hpp new file mode 100644 index 000000000..7e055e156 --- /dev/null +++ b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.hpp @@ -0,0 +1,44 @@ +/** + * @file CholeskySolverCuDssCuda.hpp + * @author Andrew Xu (xua1@ornl.gov) + * @brief Header for CUDA cuDSS implementation of Cholesky Solver using + */ + +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace ReSolve +{ + namespace hykkt + { + class CholeskySolverCuDssCuda + { + public: + CholeskySolverCuDssCuda(); + ~CholeskySolverCuDssCuda(); + + void addMatrixInfo(matrix::Csr* A); + void symbolicAnalysis(); + void numericalFactorization(real_type tol); + void solve(vector::Vector* x, vector::Vector* b); + + private: + MemoryHandler mem_; + + matrix::Csr* A_; // pointer to the input matrix + + cudssHandle_t cudss_handle_; + cudssConfig_t cudss_config_; + cudssData_t cudss_data_; + cudssMatrix_t descr_A_cudss_; + cudssMatrix_t descr_b_; + cudssMatrix_t descr_x_; + }; + } // namespace hykkt +} // namespace ReSolve diff --git a/tests/functionality/testSysRefactor.cpp b/tests/functionality/testSysRefactor.cpp index 6613736e9..96cb285d7 100644 --- a/tests/functionality/testSysRefactor.cpp +++ b/tests/functionality/testSysRefactor.cpp @@ -124,14 +124,7 @@ static int runTest(int argc, char* argv[], std::string backend, bool use_cudss) std::string refactor("none"); if (backend == "cuda") { - if (use_cudss) - { - refactor = "cudssrf"; - } - else - { - refactor = "cusolverrf"; - } + refactor = "cusolverrf"; } else if (backend == "hip") { diff --git a/tests/unit/hykkt/CMakeLists.txt b/tests/unit/hykkt/CMakeLists.txt index f05552d0f..90d825a68 100644 --- a/tests/unit/hykkt/CMakeLists.txt +++ b/tests/unit/hykkt/CMakeLists.txt @@ -31,6 +31,14 @@ target_link_libraries( resolve_vector ) +if(RESOLVE_USE_CUDSS) + add_executable(runHykktCholeskyCuDssTests.exe runHykktCholeskyCuDssTests.cpp) + target_link_libraries( + runHykktCholeskyCuDssTests.exe PRIVATE resolve_hykkt_chol_cudss resolve_matrix + resolve_vector + ) +endif() + add_executable(runHykktSCCGTests.exe runHykktSCCGTests.cpp) target_link_libraries( runHykktSCCGTests.exe @@ -62,6 +70,9 @@ set(installable_tests runHykktPermutationTests.exe runHykktRuizScalingTests.exe runHykktCholeskyTests.exe runHykktSpGEMMTests.exe runHykktSolverTests.exe ) +if(RESOLVE_USE_CUDSS) + list(APPEND installable_tests runHykktCholeskyCuDssTests.exe) +endif() install(TARGETS ${installable_tests} RUNTIME DESTINATION bin/resolve/tests/unit) @@ -75,3 +86,6 @@ add_test(NAME hykkt_chol_test COMMAND $) add_test(NAME hykkt_spgemm_test COMMAND $) add_test(NAME hykkt_sccg_test COMMAND $) add_test(NAME hykkt_solver_test COMMAND $) +if(RESOLVE_USE_CUDSS) + add_test(NAME hykkt_chol_cudss_test COMMAND $) +endif() \ No newline at end of file diff --git a/tests/unit/hykkt/HykktCholeskyCuDssTests.hpp b/tests/unit/hykkt/HykktCholeskyCuDssTests.hpp new file mode 100644 index 000000000..79810aec5 --- /dev/null +++ b/tests/unit/hykkt/HykktCholeskyCuDssTests.hpp @@ -0,0 +1,337 @@ + +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace ReSolve +{ + namespace tests + { + /** + * @brief Tests for class hykkt::CholeskySolverCuDss + * + */ + class HykktCholeskyCuDssTests : public TestBase + { + public: + HykktCholeskyCuDssTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler, std::mt19937& generator) + : memspace_(memspace), matrixHandler_(matrixHandler), generator_(generator) + { + cholmod_start(&Common); + } + + virtual ~HykktCholeskyCuDssTests() + { + cholmod_finish(&Common); + } + + /** + * @brief Test the solver on a dense 3x3 matrix with known solution. + * + * @return TestOutcome the outcome of the test + */ + TestOutcome minimalCorrectness() + { + TestStatus status; + std::string testname(__func__); + + index_type n = 3; + matrix::Csr* A = new matrix::Csr(n, n, 9); + index_type A_row_data[4] = {0, 3, 6, 9}; + index_type A_col_data[9] = {0, 1, 2, 0, 1, 2, 0, 1, 2}; + real_type A_values[9] = {4.0, 12.0, -16.0, 12.0, 37.0, -43.0, -16.0, -43.0, 98.0}; + A->allocateAll(memspace_); + A->copyFromExternal(A_row_data, A_col_data, A_values, memory::HOST, memory::HOST); + if (memspace_ == memory::DEVICE) + { + A->syncData(memory::DEVICE); + } + + ReSolve::hykkt::CholeskySolverCuDss solver(memspace_); + solver.addMatrixInfo(A); + solver.symbolicAnalysis(); + solver.setPivotTolerance(1e-12); + solver.numericalFactorization(); + vector::Vector* x = new vector::Vector(3); + x->allocateAll(memspace_); + x->setDataUpdated(memory::DEVICE); + vector::Vector* b = new vector::Vector(3); + real_type b_data[3] = {-6.0, -17.25, 30.0}; + b->allocate(memspace_); + b->copyFromExternal(b_data, memory::HOST, memspace_); + solver.solve(x, b); + + if (memspace_ == memory::DEVICE) + { + x->syncData(memory::HOST); + } + + real_type expected_x[3] = {1.0, -0.5, 0.25}; + + real_type tol = 1e-12; + for (index_type i = 0; i < n; ++i) + { + if (fabs(x->getData(memory::HOST)[i] - expected_x[i]) > tol) + { + std::cout << "Test failed at index " << i << ": expected " + << expected_x[i] << ", got " + << x->getData(memory::HOST)[i] << "\n"; + status *= false; + } + } + + delete A; + delete x; + delete b; + + return status.report(testname.c_str()); + } + + TestOutcome randomized(index_type n) + { + TestStatus status; + std::string testname(__func__); + testname += " n = " + std::to_string(n); + + cholmod_sparse* L = randomSparseLowerTriangular((size_t) n); + cholmod_sparse* L_tr = cholmod_transpose(L, 1, &Common); + cholmod_sparse* L_times_L_tr = cholmod_ssmult(L, L_tr, 0, 1, 0, &Common); + + matrix::Csr* A = new matrix::Csr((index_type) L_times_L_tr->nrow, + (index_type) L_times_L_tr->ncol, + (index_type) L_times_L_tr->nzmax); + A->allocateAll(memspace_); + A->copyFromExternal( + static_cast(L_times_L_tr->p), + static_cast(L_times_L_tr->i), + static_cast(L_times_L_tr->x), + memory::HOST, + memspace_); + + if (memspace_ == memory::DEVICE) + { + A->syncData(memory::HOST); + } + + ReSolve::hykkt::CholeskySolverCuDss solver(memspace_); + + // Add A to the solver, symbolic analysis, and numerical factorization + solver.addMatrixInfo(A); + solver.symbolicAnalysis(); + solver.numericalFactorization(); + // Generate a random vector x_expected and compute b = A * x_expected + vector::Vector* x_expected = randomVector(n); + + vector::Vector* b = new vector::Vector(n); + b->allocate(memspace_); + real_type alpha = 1.0; + real_type beta = 0.0; + matrixHandler_.matvec(A, x_expected, b, &alpha, &beta, memspace_); + + // Solve the system A * x = b + vector::Vector* x = new vector::Vector(n); + x->allocate(memspace_); + solver.solve(x, b); + + if (memspace_ == memory::DEVICE) + { + x->syncData(memory::HOST); + } + + // Verify result + real_type tol = 1e-12; + for (index_type j = 0; j < n; ++j) + { + if (fabs(x->getData(memory::HOST)[j] - x_expected->getData(memory::HOST)[j]) > tol) + { + printf("Test failed at index %d: expected %.12f, got %.12f\n, difference %.12f\n", + j, + x_expected->getData(memory::HOST)[j], + x->getData(memory::HOST)[j], + fabs(x->getData(memory::HOST)[j] - x_expected->getData(memory::HOST)[j])); + status *= false; + } + } + + cholmod_free_sparse(&L, &Common); + cholmod_free_sparse(&L_tr, &Common); + cholmod_free_sparse(&L_times_L_tr, &Common); + + delete A; + delete x_expected; + delete b; + delete x; + + return status.report(testname.c_str()); + } + + TestOutcome randomizedReuseSparsityPattern(index_type n, index_type trials) + { + TestStatus status; + std::string testname(__func__); + testname += " n = " + std::to_string(n) + ", trials = " + std::to_string(trials); + + cholmod_sparse* L = randomSparseLowerTriangular((size_t) n); + cholmod_sparse* L_tr = cholmod_transpose(L, 1, &Common); + cholmod_sparse* A_chol = cholmod_ssmult(L, L_tr, 0, 1, 0, &Common); + + matrix::Csr* A = new matrix::Csr((index_type) A_chol->nrow, (index_type) A_chol->ncol, (index_type) A_chol->nzmax); + A->allocateAll(memspace_); + A->copyFromExternal( + static_cast(A_chol->p), static_cast(A_chol->i), static_cast(A_chol->x), memory::HOST, memspace_); + + if (memspace_ == memory::DEVICE) + { + A->syncData(memory::HOST); + } + + ReSolve::hykkt::CholeskySolverCuDss solver(memspace_); + for (index_type i = 0; i < trials; ++i) + { + // Only do symbolic analysis the first iteration + solver.addMatrixInfo(A); + if (i == 0) + { + solver.symbolicAnalysis(); + } + solver.numericalFactorization(); + + // Generate a random vector x_expected and compute b = A * x_expected + vector::Vector* x_expected = randomVector(n); + + vector::Vector* b = new vector::Vector(n); + b->allocate(memspace_); + b->setToZero(memspace_); + real_type alpha = 1.0; + real_type beta = 0.0; + matrixHandler_.matvec(A, x_expected, b, &alpha, &beta, memspace_); + + // Solve the system A * x = b + vector::Vector* x = new vector::Vector(n); + x->allocate(memspace_); + solver.solve(x, b); + + if (memspace_ == memory::DEVICE) + { + x->syncData(memory::HOST); + } + + // Verify result + real_type tol = 1e-12; + for (index_type j = 0; j < n; ++j) + { + if (fabs(x->getData(memory::HOST)[j] - x_expected->getData(memory::HOST)[j]) > tol) + { + printf("Test failed at index %d: expected %.12f, got %.12f\n, difference %.12f\n", + j, + x_expected->getData(memory::HOST)[j], + x->getData(memory::HOST)[j], + fabs(x->getData(memory::HOST)[j] - x_expected->getData(memory::HOST)[j])); + status *= false; + } + } + + // reset values + for (size_t j = 0; j < L->nzmax; j++) + { + std::uniform_real_distribution distribution(-1.0, 1.0); + static_cast(L->x)[j] = 2.0 * distribution(generator_) - 1.0; + } + cholmod_free_sparse(&L_tr, &Common); + cholmod_free_sparse(&A_chol, &Common); + L_tr = cholmod_transpose(L, 1, &Common); + A_chol = cholmod_ssmult(L, L_tr, 0, 1, 0, &Common); + A->copyValues(static_cast(A_chol->x), memory::HOST, memspace_); + A->setUpdated(memspace_); + matrixHandler_.setValuesChanged(true, memspace_); + + delete b; + delete x; + delete x_expected; + } + + delete A; + cholmod_free_sparse(&L, &Common); + + return status.report(testname.c_str()); + } + + private: + ReSolve::memory::MemorySpace memspace_; + MatrixHandler& matrixHandler_; + std::mt19937& generator_; + + cholmod_common Common; + + cholmod_sparse* randomSparseLowerTriangular(size_t n) + { + double density = 2.0 / (double) n; + size_t nnz = 0; + std::vector L_p(n + 1, 0); + std::vector L_i; + std::vector L_x; + for (size_t i = 0; i < n; ++i) + { + L_p[i + 1] = L_p[i]; + for (size_t j = i; j < n; ++j) + { + std::uniform_real_distribution rand_dist(0.0, 1.0); + if (i == j || rand_dist(generator_) < density) + { + L_i.push_back((int) j); + // force diagonal entry to be non-zero + double value = 2.0; + if (i != j) + { + std::uniform_real_distribution value_dist(-1.0, 1.0); + value = 2.0 * value_dist(generator_); + value = 2.0 * static_cast(rand()) / RAND_MAX - 1.0; + } + L_x.push_back(value); + L_p[i + 1]++; + nnz++; + } + } + } + + cholmod_sparse* L = cholmod_allocate_sparse( + n, n, nnz, 1, 1, 0, CHOLMOD_REAL, &Common); + std::copy(L_p.begin(), L_p.end(), static_cast(L->p)); + std::copy(L_i.begin(), L_i.end(), static_cast(L->i)); + std::copy(L_x.begin(), L_x.end(), static_cast(L->x)); + + return L; + } + + vector::Vector* randomVector(index_type n) + { + vector::Vector* v = new vector::Vector(n); + v->allocateAll(memspace_); + std::uniform_real_distribution distribution(0.0, 1.0); + for (index_type i = 0; i < n; ++i) + { + v->getData(memory::HOST)[i] = distribution(generator_); + } + v->setDataUpdated(memory::HOST); + if (memspace_ == memory::DEVICE) + { + v->syncData(memory::DEVICE); + } + return v; + } + }; // class HykktCholeskyCuDssTests + } // namespace tests +} // namespace ReSolve diff --git a/tests/unit/hykkt/HykktCholeskyTests.hpp b/tests/unit/hykkt/HykktCholeskyTests.hpp index 99eca2618..e04d39aeb 100644 --- a/tests/unit/hykkt/HykktCholeskyTests.hpp +++ b/tests/unit/hykkt/HykktCholeskyTests.hpp @@ -27,8 +27,8 @@ namespace ReSolve class HykktCholeskyTests : public TestBase { public: - HykktCholeskyTests(bool use_cudss, memory::MemorySpace memspace, MatrixHandler& matrixHandler, std::mt19937& generator) - : use_cudss_(use_cudss), memspace_(memspace), matrixHandler_(matrixHandler), generator_(generator) + HykktCholeskyTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler, std::mt19937& generator) + : memspace_(memspace), matrixHandler_(matrixHandler), generator_(generator) { cholmod_start(&Common); } @@ -60,7 +60,7 @@ namespace ReSolve A->syncData(memory::DEVICE); } - ReSolve::hykkt::CholeskySolver solver(use_cudss_, memspace_); + ReSolve::hykkt::CholeskySolver solver(memspace_); solver.addMatrixInfo(A); solver.symbolicAnalysis(); solver.setPivotTolerance(1e-12); @@ -273,7 +273,6 @@ namespace ReSolve ReSolve::memory::MemorySpace memspace_; MatrixHandler& matrixHandler_; std::mt19937& generator_; - bool use_cudss_; cholmod_common Common; diff --git a/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp b/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp new file mode 100644 index 000000000..ffca10b2e --- /dev/null +++ b/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp @@ -0,0 +1,55 @@ +/** + * @file runHykktCholeskyCuDssTests.hpp + * @author Shaked Regev (regevs@ornl.gov) + * @author Adham Ibrahim (ibrahimas@ornl.gov) + * @brief Tests for class hykkt::CholeskySolver + * + */ +#include +#include +#include +#include + +#include "resolve/Common.hpp" +#include "tests/unit/hykkt/HykktCholeskyCuDssTests.hpp" + +/** + * @brief Run tests with a given backend + * + * @param backend - string name of the hardware backend + * @param result - test results + */ +template +void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace, ReSolve::tests::TestingResults& result) +{ + std::cout << "Running tests on " << backend << " device:\n"; + + WorkspaceType workspace; + workspace.initializeHandles(); + ReSolve::MatrixHandler handler(&workspace); + std::mt19937 generator(ReSolve::constants::SEED); // set random seed for reproducibility + ReSolve::tests::HykktCholeskyCuDssTests test(memspace, handler, generator); + + result += test.minimalCorrectness(); + handler.setValuesChanged(true, memspace); + workspace.resetLinAlgWorkspace(); // reset is necessary due to different sparsity. + + for (int size : {3, 10, 100, 1000}) + { + result += test.randomized(size); + handler.setValuesChanged(true, memspace); + workspace.resetLinAlgWorkspace(); + } + + result += test.randomizedReuseSparsityPattern(3, 10); + + std::cout << "\n"; +} + +int main(int, char**) +{ + ReSolve::tests::TestingResults result; + runTests("CUDA", ReSolve::memory::DEVICE, result); + + return result.summary(); +} diff --git a/tests/unit/hykkt/runHykktCholeskyTests.cpp b/tests/unit/hykkt/runHykktCholeskyTests.cpp index dfe4dc1be..55ebeb7d6 100644 --- a/tests/unit/hykkt/runHykktCholeskyTests.cpp +++ b/tests/unit/hykkt/runHykktCholeskyTests.cpp @@ -20,23 +20,15 @@ * @param result - test results */ template -void runTests(bool use_cudss, const std::string& backend, ReSolve::memory::MemorySpace memspace, ReSolve::tests::TestingResults& result) +void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace, ReSolve::tests::TestingResults& result) { - if (backend == "CUDA") - { - std::string impl = use_cudss ? "cuDSS" : "CuSolver"; - std::cout << "Running tests on " << backend << " device with " << impl << " implementation:\n"; - } - else - { - std::cout << "Running tests on " << backend << " device:\n"; - } + std::cout << "Running tests on " << backend << " device:\n"; WorkspaceType workspace; workspace.initializeHandles(); ReSolve::MatrixHandler handler(&workspace); std::mt19937 generator(ReSolve::constants::SEED); // set random seed for reproducibility - ReSolve::tests::HykktCholeskyTests test(use_cudss, memspace, handler, generator); + ReSolve::tests::HykktCholeskyTests test(memspace, handler, generator); result += test.minimalCorrectness(); handler.setValuesChanged(true, memspace); @@ -57,17 +49,14 @@ void runTests(bool use_cudss, const std::string& backend, ReSolve::memory::Memor int main(int, char**) { ReSolve::tests::TestingResults result; - runTests(false, "CPU", ReSolve::memory::HOST, result); + runTests("CPU", ReSolve::memory::HOST, result); #ifdef RESOLVE_USE_CUDA - runTests(false, "CUDA", ReSolve::memory::DEVICE, result); -#ifdef RESOLVE_USE_CUDSS -#endif - runTests(true, "CUDA", ReSolve::memory::DEVICE, result); + runTests("CUDA", ReSolve::memory::DEVICE, result); #endif #ifdef RESOLVE_USE_HIP - runTests(false, "HIP", ReSolve::memory::DEVICE, result); + runTests("HIP", ReSolve::memory::DEVICE, result); #endif return result.summary(); From 3380b906ab6111f8d1c07177fcc02c88f47b3630 Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Fri, 24 Jul 2026 16:18:32 +0000 Subject: [PATCH 14/18] Apply pre-commmit fixes --- resolve/hykkt/cholesky_cudss/CMakeLists.txt | 18 ++++++---- .../cholesky_cudss/CholeskySolverCuDss.cpp | 2 +- .../cholesky_cudss/CholeskySolverCuDss.hpp | 4 +-- .../CholeskySolverCuDssCuda.cpp | 36 +++++++++---------- .../CholeskySolverCuDssCuda.hpp | 3 +- tests/unit/hykkt/CMakeLists.txt | 10 +++--- .../unit/hykkt/runHykktCholeskyCuDssTests.cpp | 4 +-- 7 files changed, 42 insertions(+), 35 deletions(-) diff --git a/resolve/hykkt/cholesky_cudss/CMakeLists.txt b/resolve/hykkt/cholesky_cudss/CMakeLists.txt index ef26c395e..05db376c5 100644 --- a/resolve/hykkt/cholesky_cudss/CMakeLists.txt +++ b/resolve/hykkt/cholesky_cudss/CMakeLists.txt @@ -9,27 +9,31 @@ set(HyKKT_CHOL_CUDSS_SRC CholeskySolverCuDss.cpp CholeskySolverCuDssCuda.cpp) # Header files to be installed -set(HyKKT_CHOL_CUDSS_HEADER_INSTALL - CholeskySolverCuDss.hpp CholeskySolverCuDssCuda.hpp +set(HyKKT_CHOL_CUDSS_HEADER_INSTALL CholeskySolverCuDss.hpp + CholeskySolverCuDssCuda.hpp ) add_library(resolve_hykkt_chol_cudss SHARED ${HyKKT_CHOL_CUDSS_SRC}) target_link_libraries(resolve_hykkt_chol_cudss PUBLIC ${suitesparse_cholmod}) -target_include_directories(resolve_hykkt_chol_cudss PUBLIC ${SUITESPARSE_INCLUDE_DIR}) +target_include_directories( + resolve_hykkt_chol_cudss PUBLIC ${SUITESPARSE_INCLUDE_DIR} +) # Link to CUDA ReSolve backend if CUDA is support enabled target_sources(resolve_hykkt_chol_cudss PRIVATE ${HyKKT_CHOL_CUDSS_SRC}) target_link_libraries(resolve_hykkt_chol_cudss PUBLIC resolve_backend_cuda) target_link_libraries( - resolve_hykkt_chol_cudss PUBLIC resolve_workspace resolve_vector resolve_matrix - resolve_logger + resolve_hykkt_chol_cudss PUBLIC resolve_workspace resolve_vector + resolve_matrix resolve_logger ) target_include_directories( resolve_hykkt_chol_cudss INTERFACE $ - $ + $ ) -install(FILES ${HyKKT_CHOL_CUDSS_HEADER_INSTALL} DESTINATION include/resolve/hykkt) +install(FILES ${HyKKT_CHOL_CUDSS_HEADER_INSTALL} + DESTINATION include/resolve/hykkt +) diff --git a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp index e0b2ae632..87bfc879c 100644 --- a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp +++ b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp @@ -21,7 +21,7 @@ namespace ReSolve */ CholeskySolverCuDss::CholeskySolverCuDss(memory::MemorySpace memspace) : memspace_(memspace), - impl_(new CholeskySolverCuDssCuda()) + impl_(new CholeskySolverCuDssCuda()) { } diff --git a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.hpp b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.hpp index 4e3a11110..2dd90d512 100644 --- a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.hpp +++ b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.hpp @@ -29,8 +29,8 @@ namespace ReSolve private: memory::MemorySpace memspace_; - matrix::Csr* A_; - real_type tol_ = 1e-12; + matrix::Csr* A_; + real_type tol_ = 1e-12; CholeskySolverCuDssCuda* impl_; }; } // namespace hykkt diff --git a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.cpp b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.cpp index fa52c6eb0..853c467a8 100644 --- a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.cpp +++ b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.cpp @@ -35,18 +35,18 @@ namespace ReSolve A_ = A; cudssMatrixCreateCsr(&descr_A_cudss_, - A_->getNumRows(), - A_->getNumColumns(), - A_->getNnz(), - A_->getRowData(memory::DEVICE), - nullptr, // Row end offsets (null for standard CSR) - A_->getColData(memory::DEVICE), - A_->getValues(memory::DEVICE), - CUDA_R_32I, - CUDA_R_64F, - CUDSS_MTYPE_SPD, - CUDSS_MVIEW_LOWER, - CUDSS_BASE_ZERO); + A_->getNumRows(), + A_->getNumColumns(), + A_->getNnz(), + A_->getRowData(memory::DEVICE), + nullptr, // Row end offsets (null for standard CSR) + A_->getColData(memory::DEVICE), + A_->getValues(memory::DEVICE), + CUDA_R_32I, + CUDA_R_64F, + CUDSS_MTYPE_SPD, + CUDSS_MVIEW_LOWER, + CUDSS_BASE_ZERO); } /** @@ -69,12 +69,12 @@ namespace ReSolve CUDA_R_64F, CUDSS_LAYOUT_COL_MAJOR); cudssExecute(cudss_handle_, - CUDSS_PHASE_ANALYSIS, - cudss_config_, - cudss_data_, - descr_A_cudss_, - descr_x_, - descr_b_); + CUDSS_PHASE_ANALYSIS, + cudss_config_, + cudss_data_, + descr_A_cudss_, + descr_x_, + descr_b_); } /** diff --git a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.hpp b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.hpp index 7e055e156..efc239992 100644 --- a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.hpp +++ b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDssCuda.hpp @@ -6,12 +6,13 @@ #pragma once +#include + #include #include #include #include #include -#include namespace ReSolve { diff --git a/tests/unit/hykkt/CMakeLists.txt b/tests/unit/hykkt/CMakeLists.txt index 90d825a68..19dc5db53 100644 --- a/tests/unit/hykkt/CMakeLists.txt +++ b/tests/unit/hykkt/CMakeLists.txt @@ -34,8 +34,8 @@ target_link_libraries( if(RESOLVE_USE_CUDSS) add_executable(runHykktCholeskyCuDssTests.exe runHykktCholeskyCuDssTests.cpp) target_link_libraries( - runHykktCholeskyCuDssTests.exe PRIVATE resolve_hykkt_chol_cudss resolve_matrix - resolve_vector + runHykktCholeskyCuDssTests.exe PRIVATE resolve_hykkt_chol_cudss + resolve_matrix resolve_vector ) endif() @@ -87,5 +87,7 @@ add_test(NAME hykkt_spgemm_test COMMAND $) add_test(NAME hykkt_sccg_test COMMAND $) add_test(NAME hykkt_solver_test COMMAND $) if(RESOLVE_USE_CUDSS) - add_test(NAME hykkt_chol_cudss_test COMMAND $) -endif() \ No newline at end of file + add_test(NAME hykkt_chol_cudss_test + COMMAND $ + ) +endif() diff --git a/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp b/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp index ffca10b2e..65fe5e95b 100644 --- a/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp +++ b/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp @@ -26,8 +26,8 @@ void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace, WorkspaceType workspace; workspace.initializeHandles(); - ReSolve::MatrixHandler handler(&workspace); - std::mt19937 generator(ReSolve::constants::SEED); // set random seed for reproducibility + ReSolve::MatrixHandler handler(&workspace); + std::mt19937 generator(ReSolve::constants::SEED); // set random seed for reproducibility ReSolve::tests::HykktCholeskyCuDssTests test(memspace, handler, generator); result += test.minimalCorrectness(); From 25679f531f60f5be9aa997cabe0d66772d886aad Mon Sep 17 00:00:00 2001 From: Shaked Regev <35384901+shakedregev@users.noreply.github.com> Date: Fri, 24 Jul 2026 13:51:14 -0400 Subject: [PATCH 15/18] Minor authorship changes Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com> --- resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp | 1 - tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp index 87bfc879c..b5da436fb 100644 --- a/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp +++ b/resolve/hykkt/cholesky_cudss/CholeskySolverCuDss.cpp @@ -1,6 +1,5 @@ /** * @file CholeskySolverCuDss.cpp - * @author Adham Ibrahim (ibrahimas@ornl.gov) * @brief Cholesky decomposition solver CuDSS implementation. This is a CUDA-only variant of CholeskySolver. */ diff --git a/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp b/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp index 65fe5e95b..b46bfe78c 100644 --- a/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp +++ b/tests/unit/hykkt/runHykktCholeskyCuDssTests.cpp @@ -1,7 +1,7 @@ /** * @file runHykktCholeskyCuDssTests.hpp * @author Shaked Regev (regevs@ornl.gov) - * @author Adham Ibrahim (ibrahimas@ornl.gov) + * @author Andrew Xu (xua1@ornl.gov) * @brief Tests for class hykkt::CholeskySolver * */ From 8c63206dd0efdc947bb8530551c856326414edbc Mon Sep 17 00:00:00 2001 From: Andrew Xu Date: Fri, 24 Jul 2026 18:23:54 +0000 Subject: [PATCH 16/18] Make separate tests for cuDSS refactor --- tests/functionality/CMakeLists.txt | 22 ++ tests/functionality/testCuDssRefactor.cpp | 268 ++++++++++++++++++++++ tests/functionality/testSysRefactor.cpp | 14 +- 3 files changed, 294 insertions(+), 10 deletions(-) create mode 100644 tests/functionality/testCuDssRefactor.cpp diff --git a/tests/functionality/CMakeLists.txt b/tests/functionality/CMakeLists.txt index c36e18a24..2790e509a 100644 --- a/tests/functionality/CMakeLists.txt +++ b/tests/functionality/CMakeLists.txt @@ -37,6 +37,12 @@ if(RESOLVE_USE_KLU) target_link_libraries(sys_glu_test.exe PRIVATE ReSolve) endif(RESOLVE_USE_CUDA) + + if(RESOLVE_USE_CUDSS) + # System solver test with cuDSS rf and iterative refinement + add_executable(cudss_refactor_test.exe testCuDssRefactor.cpp) + target_link_libraries(cudss_refactor_test.exe PRIVATE ReSolve) + endif(RESOLVE_USE_CUDSS) endif(RESOLVE_USE_KLU) if(RESOLVE_USE_LUSOL) @@ -291,6 +297,22 @@ if(RESOLVE_USE_KLU) "${test_data_dir}" "-m" "${asym_matrix_path}" "-r" "${asym_rhs_path}" ) endif(RESOLVE_USE_CUDA) + + # cUDSS specific tests + if(RESOLVE_USE_CUDSS) + add_test( + NAME cudss_refactor_cuda_test + COMMAND + $ "${test_data_dir}" "-d" + "${test_data_dir}" "-m" "${sym_matrix_path}" "-r" "${sym_rhs_path}" + ) + add_test( + NAME cudss_refactor_cuda_asym_test + COMMAND + $ "${test_data_dir}" "-d" + "${test_data_dir}" "-m" "${asym_matrix_path}" "-r" "${asym_rhs_path}" + ) + endif(RESOLVE_USE_CUDSS) # ROCm specific tests if(RESOLVE_USE_HIP) diff --git a/tests/functionality/testCuDssRefactor.cpp b/tests/functionality/testCuDssRefactor.cpp new file mode 100644 index 000000000..868c7a9f0 --- /dev/null +++ b/tests/functionality/testCuDssRefactor.cpp @@ -0,0 +1,268 @@ + +/** + * @file testSysHipRefine.cpp + * @author Kasia Swirydowicz (kasia.swirydowicz@pnnl.gov) + * @author Slaven Peles (peless@ornl.gov) + * @brief Functionality test for SystemSolver class + * @date 2023-12-14 + * + * + */ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "TestHelper.hpp" + +template +static int runTest(int argc, char* argv[], std::string backend); + +int main(int argc, char* argv[]) +{ + int error_sum = 0; + + // Refactorization on CPU not currently supported in SystemSolver class + // error_sum += runTest(argc, argv, "cpu"); + + error_sum += runTest(argc, argv, "cuda"); + + return error_sum; +} + +template +static int runTest(int argc, char* argv[], std::string backend) +{ + // Use ReSolve data types. + using namespace ReSolve; + using real_type = ReSolve::real_type; + using index_type = ReSolve::index_type; + using vector_type = ReSolve::vector::Vector; + + // Error sum needs to be 0 at the end for test to PASS. + // It is a FAIL otheriwse. + int error_sum = 0; + int status = 0; + + memory::MemorySpace memspace = memory::HOST; + if (backend != "cpu") + { + memspace = memory::DEVICE; + } + + // Collect all command line options + ReSolve::CliOptions options(argc, argv); + + // Get directory with input files + auto opt = options.getParamFromKey("-d"); + std::string data_path = opt ? (*opt).second : "."; + + // Get matrix file name + opt = options.getParamFromKey("-m"); + if (!opt) + { + std::cout << "Matrix file name not provided. Use -m .\n"; + return -1; + } + std::string matrix_temp = (*opt).second; + + // Get rhs file name + opt = options.getParamFromKey("-r"); + if (!opt) + { + std::cout << "RHS file name not provided. Use -r .\n"; + return -1; + } + std::string rhs_temp = (*opt).second; + + // Construct matrix and rhs file names from inputs + std::string matrix_file_name_1 = data_path + matrix_temp + "01.mtx"; + std::string matrix_file_name_2 = data_path + matrix_temp + "02.mtx"; + std::string rhs_file_name_1 = data_path + rhs_temp + "01.mtx"; + std::string rhs_file_name_2 = data_path + rhs_temp + "02.mtx"; + + // Create workspace and initialize its handles. + workspace_type workspace; + workspace.initializeHandles(); + + // Create test helper + TestHelper helper(workspace); + + // Create system solver + std::string refactor("none"); + refactor = "cudssrf"; + ReSolve::SystemSolver solver(&workspace, + "klu", // factorization + refactor, // refactorization + refactor, // triangular solve + "none", // preconditioner (always 'none' here) + "none"); // iterative refinement + + // Configure solver (CUDA-based solver needs slightly different + // settings than HIP-based one) + solver.setRefinementMethod("fgmres", "cgs2"); + solver.getIterativeSolver().setCliParam("restart", "100"); + solver.getIterativeSolver().setTol(ReSolve::constants::MACHINE_EPSILON); + solver.getIterativeSolver().setMaxit(400); + + // Read first matrix + std::ifstream mat1(matrix_file_name_1); + if (!mat1.is_open()) + { + std::cout << "Failed to open file " << matrix_file_name_1 << "\n"; + return -1; + } + ReSolve::matrix::Csr* A = ReSolve::io::createCsrFromFile(mat1, true); + if (memspace == memory::DEVICE) + { + A->allocateMatrixData(memory::DEVICE); + A->syncData(memory::DEVICE); + } + mat1.close(); + + // Read first rhs vector + std::ifstream rhs1_file(rhs_file_name_1); + if (!rhs1_file.is_open()) + { + std::cout << "Failed to open file " << rhs_file_name_1 << "\n"; + return -1; + } + real_type* rhs = ReSolve::io::createArrayFromFile(rhs1_file); + rhs1_file.close(); + + // Create and set residual vector + vector_type vec_rhs = (A->getNumRows()); + vec_rhs.allocate(ReSolve::memory::HOST); + vec_rhs.copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); + + // Create and allocate solution vector + vector_type vec_x(A->getNumRows()); + vec_x.allocateAll(memspace); + + // Add system matrix to the solver + status = solver.setMatrix(A); + error_sum += status; + + // Solve the first system using KLU + status = solver.analyze(); + error_sum += status; + + status = solver.factorize(); + error_sum += status; + + status = solver.solve(&vec_rhs, &vec_x); + error_sum += status; + + // Compute error norms for the system + if (memspace == ReSolve::memory::DEVICE) + { + vec_x.syncData(ReSolve::memory::DEVICE); + vec_rhs.syncData(ReSolve::memory::DEVICE); + } + helper.setSystem(A, &vec_rhs, &vec_x); + + // Print result summary and check solution + std::cout << "\nResults (first matrix): \n\n"; + helper.printSummary(); + error_sum += helper.checkResult(ReSolve::constants::MACHINE_EPSILON); + + // Verify norm of scaled residuals calculation in SystemSolver class + real_type nsr_system = solver.getNormOfScaledResiduals(&vec_rhs, &vec_x); + error_sum += helper.checkNormOfScaledResiduals(nsr_system); + + // Verify relative residual norm computation in SystemSolver + real_type rel_residual_norm = solver.getResidualNorm(&vec_rhs, &vec_x); + error_sum += helper.checkRelativeResidualNorm(rel_residual_norm); + + // Now prepare the Rf solver + status = solver.refactorizationSetup(); + error_sum += status; + + // Load the second matrix + std::ifstream mat2(matrix_file_name_2); + if (!mat2.is_open()) + { + std::cout << "Failed to open file " << matrix_file_name_2 << "\n"; + return -1; + } + ReSolve::io::updateMatrixFromFile(mat2, A); + if (memspace != memory::HOST) + { + A->syncData(memspace); + } + mat2.close(); + + // Load the second rhs vector + std::ifstream rhs2_file(rhs_file_name_2); + if (!rhs2_file.is_open()) + { + std::cout << "Failed to open file " << rhs_file_name_2 << "\n"; + return -1; + } + ReSolve::io::updateArrayFromFile(rhs2_file, &rhs); + rhs2_file.close(); + + if (memspace == ReSolve::memory::DEVICE) + { + vec_rhs.copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); + } + vec_rhs.copyFromExternal(rhs, ReSolve::memory::HOST, memspace); + + // Refactorize matrix + status = solver.refactorize(); + error_sum += status; + + // Solve system + status = solver.solve(&vec_rhs, &vec_x); + error_sum += status; + + // Compute error norms for the system + if (memspace == ReSolve::memory::DEVICE) + { + if (!vec_rhs.isUpdated(ReSolve::memory::HOST)) + { + vec_rhs.syncData(ReSolve::memory::HOST); + } + vec_x.syncData(ReSolve::memory::HOST); + } + helper.resetSystem(A, &vec_rhs, &vec_x); + + // Print result summary and check solution + std::cout << "\nResults (second matrix): \n\n"; + helper.printSummary(); + helper.printIrSummary(&(solver.getIterativeSolver())); + error_sum += helper.checkResult(ReSolve::constants::MACHINE_EPSILON); + + // Verify norm of scaled residuals calculation in SystemSolver class + nsr_system = solver.getNormOfScaledResiduals(&vec_rhs, &vec_x); + error_sum += helper.checkNormOfScaledResiduals(nsr_system); + + // Verify relative residual norm computation in SystemSolver + rel_residual_norm = solver.getResidualNorm(&vec_rhs, &vec_x); + error_sum += helper.checkRelativeResidualNorm(rel_residual_norm); + + // Add one output specific to GMRES + index_type restart = solver.getIterativeSolver().getCliParamInt("restart"); + std::cout << "\t IR GMRES restart : " << restart << "\n"; + + isTestPass(error_sum, "Test SystemSolver: KLU with Rf and IR"); + + delete A; + delete[] rhs; + + return error_sum; +} diff --git a/tests/functionality/testSysRefactor.cpp b/tests/functionality/testSysRefactor.cpp index 96cb285d7..86d629057 100644 --- a/tests/functionality/testSysRefactor.cpp +++ b/tests/functionality/testSysRefactor.cpp @@ -27,9 +27,6 @@ #ifdef RESOLVE_USE_CUDA #include -#ifdef RESOLVE_USE_CUDSS -#include -#endif #endif #ifdef RESOLVE_USE_HIP @@ -39,7 +36,7 @@ #include "TestHelper.hpp" template -static int runTest(int argc, char* argv[], std::string backend, bool use_cudss); +static int runTest(int argc, char* argv[], std::string backend); int main(int argc, char* argv[]) { @@ -49,21 +46,18 @@ int main(int argc, char* argv[]) // error_sum += runTest(argc, argv, "cpu"); #ifdef RESOLVE_USE_CUDA - error_sum += runTest(argc, argv, "cuda", false); -#ifdef RESOLVE_USE_CUDSS - error_sum += runTest(argc, argv, "cuda", true); -#endif + error_sum += runTest(argc, argv, "cuda"); #endif #ifdef RESOLVE_USE_HIP - error_sum += runTest(argc, argv, "hip", false); + error_sum += runTest(argc, argv, "hip"); #endif return error_sum; } template -static int runTest(int argc, char* argv[], std::string backend, bool use_cudss) +static int runTest(int argc, char* argv[], std::string backend) { // Use ReSolve data types. using namespace ReSolve; From ca98f196931cde29e898652d9b5b41a6367db1d0 Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Fri, 24 Jul 2026 18:24:27 +0000 Subject: [PATCH 17/18] Apply pre-commmit fixes --- tests/functionality/CMakeLists.txt | 2 +- tests/functionality/testCuDssRefactor.cpp | 6 ++---- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/tests/functionality/CMakeLists.txt b/tests/functionality/CMakeLists.txt index 2790e509a..3770dd3c3 100644 --- a/tests/functionality/CMakeLists.txt +++ b/tests/functionality/CMakeLists.txt @@ -297,7 +297,7 @@ if(RESOLVE_USE_KLU) "${test_data_dir}" "-m" "${asym_matrix_path}" "-r" "${asym_rhs_path}" ) endif(RESOLVE_USE_CUDA) - + # cUDSS specific tests if(RESOLVE_USE_CUDSS) add_test( diff --git a/tests/functionality/testCuDssRefactor.cpp b/tests/functionality/testCuDssRefactor.cpp index 868c7a9f0..d96116f31 100644 --- a/tests/functionality/testCuDssRefactor.cpp +++ b/tests/functionality/testCuDssRefactor.cpp @@ -12,6 +12,8 @@ #include #include +#include "TestHelper.hpp" +#include #include #include #include @@ -25,10 +27,6 @@ #include #include -#include - -#include "TestHelper.hpp" - template static int runTest(int argc, char* argv[], std::string backend); From 196d4df88df90de4c681fa9b140b4f33710a883c Mon Sep 17 00:00:00 2001 From: Shaked Regev <35384901+shakedregev@users.noreply.github.com> Date: Fri, 24 Jul 2026 14:53:47 -0400 Subject: [PATCH 18/18] Apply suggestion from @shakedregev --- examples/cudssRefactor.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/cudssRefactor.cpp b/examples/cudssRefactor.cpp index c57dd5480..416fe78c2 100644 --- a/examples/cudssRefactor.cpp +++ b/examples/cudssRefactor.cpp @@ -145,7 +145,7 @@ int cudssRefactor(int argc, char* argv[]) workspace_type workspace; workspace.initializeHandles(); - // Create a helper object (computing errors, printing summaries, etc.) + // Create a helper object (for computing errors, printing summaries, etc.) ExampleHelper helper(workspace); std::cout << "cudssRefactor with CUDA backend\n";