From 2ec4c5d603640d7d0094edcbed70acb49fa9526e Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Thu, 21 May 2026 16:33:29 -0700 Subject: [PATCH 01/11] use caller-provided backend-specific handlers in SCCG --- .../hykkt/sccg/SchurComplementConjugateGradient.cpp | 9 +++++---- .../hykkt/sccg/SchurComplementConjugateGradient.hpp | 13 +++---------- tests/unit/hykkt/HykktSCCGTests.hpp | 12 ++++++------ tests/unit/hykkt/runHykktSCCGTests.cpp | 7 ++++--- 4 files changed, 18 insertions(+), 23 deletions(-) diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index 3e81d116..bdeacce3 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -18,15 +18,16 @@ namespace ReSolve index_type n, index_type m, CholeskySolver* choleskySolver, - memory::MemorySpace memspace) + memory::MemorySpace memspace, + MatrixHandler& matrixHandler, + VectorHandler& vectorHandler) : n_(n), m_(m), choleskySolver_(choleskySolver), memspace_(memspace), mem_(), - workspace_(), - matrixhandler_(&workspace_), - vectorhandler_(&workspace_), + matrixhandler_(matrixHandler), + vectorhandler_(vectorHandler), y_(m_), z_(m_), r_(n_), diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp index 88664c6e..a1b75221 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp @@ -24,7 +24,7 @@ namespace ReSolve class SchurComplementConjugateGradient { public: - SchurComplementConjugateGradient(index_type n, index_type m, CholeskySolver* choleskySolver, memory::MemorySpace memspace); + SchurComplementConjugateGradient(index_type n, index_type m, CholeskySolver* choleskySolver, memory::MemorySpace memspace, MatrixHandler& matrixHandler, VectorHandler& vectorHandler); void addMatrixInfo(matrix::Csr* jc, matrix::Csr* jc_tr); void addVectorInfo(vector::Vector* x0, vector::Vector* b); @@ -41,15 +41,8 @@ namespace ReSolve int itmax_ = 100; // Maximum iterations for conjugate gradient double tol_ = 1e-12; // Solver tolerance for Schur -#ifdef RESOLVE_USE_CUDA - LinAlgWorkspaceCUDA workspace_; -#elif defined(RESOLVE_USE_HIP) - LinAlgWorkspaceHIP workspace_; -#else - LinAlgWorkspaceCpu workspace_; -#endif - MatrixHandler matrixhandler_; - VectorHandler vectorhandler_; + MatrixHandler& matrixhandler_; + VectorHandler& vectorhandler_; CholeskySolver* choleskySolver_; // Cholesky factorization on 1,1 block diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index c26cbc34..18727207 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -26,8 +27,8 @@ namespace ReSolve class HykktSchurComplementConjugateGradientTests : public TestBase { public: - HykktSchurComplementConjugateGradientTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler) - : memspace_(memspace), matrixHandler_(matrixHandler) + HykktSchurComplementConjugateGradientTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler, VectorHandler& vectorHandler) + : memspace_(memspace), matrixHandler_(matrixHandler), vectorHandler_(vectorHandler) { } @@ -68,14 +69,12 @@ namespace ReSolve index_type n = jc->getNumRows(); index_type m = jc->getNumColumns(); index_type nnz = jc->getNnz(); - hykkt::SchurComplementConjugateGradient sccg(n, m, &choleskySolver, memspace_); + hykkt::SchurComplementConjugateGradient sccg(n, m, &choleskySolver, memspace_, matrixHandler_, vectorHandler_); sccg.setSolverTolerance(tol); - LinAlgWorkspaceCpu workspace; - MatrixHandler matrix_handler(&workspace); matrix::Csr* jc_tr = new matrix::Csr(m, n, nnz); jc_tr->allocateMatrixData(memspace_); - matrix_handler.transpose(jc, jc_tr, memspace_); + matrixHandler_.transpose(jc, jc_tr, memspace_); vector::Vector* x0 = new vector::Vector(n); x0->allocate(memspace_); @@ -107,6 +106,7 @@ namespace ReSolve private: memory::MemorySpace memspace_; MatrixHandler& matrixHandler_; + VectorHandler& vectorHandler_; /** * @brief Generate a random vector of doubles between 0 and RAND_MAX. Copied from HykktCholeskyTests.hpp. diff --git a/tests/unit/hykkt/runHykktSCCGTests.cpp b/tests/unit/hykkt/runHykktSCCGTests.cpp index c4a05246..9fc32f81 100644 --- a/tests/unit/hykkt/runHykktSCCGTests.cpp +++ b/tests/unit/hykkt/runHykktSCCGTests.cpp @@ -8,6 +8,7 @@ #include #include +#include #include #ifdef RESOLVE_USE_CUDA #include @@ -32,9 +33,9 @@ void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace, WorkspaceType workspace; workspace.initializeHandles(); - ReSolve::MatrixHandler handler(&workspace); - - ReSolve::tests::HykktSchurComplementConjugateGradientTests test(memspace, handler); + ReSolve::MatrixHandler matrixHandler(&workspace); + ReSolve::VectorHandler vectorHandler(&workspace); + ReSolve::tests::HykktSchurComplementConjugateGradientTests test(memspace, matrixHandler, vectorHandler); result += test.SCCGTest(); From 65c0903ba82e7daba97545e6499082e7d8d51159 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 22 May 2026 10:43:55 -0700 Subject: [PATCH 02/11] load SCCG data on host before device sync --- tests/unit/hykkt/HykktSCCGTests.hpp | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index 18727207..e1db0481 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -54,8 +54,12 @@ namespace ReSolve std::ifstream bFile(bFileName); matrix::Csr* h = new matrix::Csr(2278, 2278, 11304, true, false); - h->allocateMatrixData(memspace_); + h->allocateMatrixData(memory::HOST); io::updateMatrixFromFile(hFile, h); + if (memspace_ == memory::DEVICE) + { + h->syncData(memory::DEVICE); + } hykkt::CholeskySolver choleskySolver(memspace_); choleskySolver.addMatrixInfo(h); choleskySolver.symbolicAnalysis(); @@ -63,8 +67,12 @@ namespace ReSolve choleskySolver.numericalFactorization(); matrix::Csr* jc = new matrix::Csr(1386, 2278, 6784, false, false); - jc->allocateMatrixData(memspace_); + jc->allocateMatrixData(memory::HOST); io::updateMatrixFromFile(jcFile, jc); + if (memspace_ == memory::DEVICE) + { + jc->syncData(memory::DEVICE); + } index_type n = jc->getNumRows(); index_type m = jc->getNumColumns(); @@ -77,12 +85,16 @@ namespace ReSolve matrixHandler_.transpose(jc, jc_tr, memspace_); vector::Vector* x0 = new vector::Vector(n); - x0->allocate(memspace_); + x0->allocate(memory::HOST); randomVector(x0); vector::Vector* b = new vector::Vector(n); - b->allocate(memspace_); + b->allocate(memory::HOST); io::updateVectorFromFile(bFile, b); + if (memspace_ == memory::DEVICE) + { + b->syncData(memory::DEVICE); + } sccg.addMatrixInfo(jc, jc_tr); sccg.addVectorInfo(x0, b); From 63515edd8576840bf9a74a943698afc6f28a6af5 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 22 May 2026 15:04:24 -0700 Subject: [PATCH 03/11] refresh CUDA SpMV cache when matrix changes --- resolve/matrix/MatrixHandlerCuda.cpp | 43 +++++++++++++---------- resolve/matrix/MatrixHandlerCuda.hpp | 5 +++ resolve/workspace/LinAlgWorkspaceCUDA.cpp | 14 ++++++++ resolve/workspace/LinAlgWorkspaceCUDA.hpp | 1 + 4 files changed, 45 insertions(+), 18 deletions(-) diff --git a/resolve/matrix/MatrixHandlerCuda.cpp b/resolve/matrix/MatrixHandlerCuda.cpp index 7b310843..4b6b1cce 100644 --- a/resolve/matrix/MatrixHandlerCuda.cpp +++ b/resolve/matrix/MatrixHandlerCuda.cpp @@ -76,34 +76,34 @@ namespace ReSolve // result = alpha *A*x + beta * result cusparseStatus_t status; cusparseDnVecDescr_t vecx = workspace_->getVecX(); - cusparseCreateDnVec(&vecx, A->getNumRows(), vec_x->getData(memory::DEVICE), CUDA_R_64F); + cusparseCreateDnVec(&vecx, A->getNumColumns(), vec_x->getData(memory::DEVICE), CUDA_R_64F); cusparseDnVecDescr_t vecAx = workspace_->getVecY(); cusparseCreateDnVec(&vecAx, A->getNumRows(), vec_result->getData(memory::DEVICE), CUDA_R_64F); - cusparseSpMatDescr_t matA = workspace_->getSpmvMatrixDescriptor(); - - void* buffer_spmv = workspace_->getSpmvBuffer(); cusparseHandle_t handle_cusparse = workspace_->getCusparseHandle(); - if (values_changed_) + bool matrix_changed = (matrix_for_matvec_ != A) || (matvec_num_rows_ != A->getNumRows()) || (matvec_num_cols_ != A->getNumColumns()) || (matvec_nnz_ != A->getNnz()); + if (matrix_changed || values_changed_) { - status = cusparseCreateCsr(&matA, - A->getNumRows(), - A->getNumColumns(), - A->getNnz(), - A->getRowData(memory::DEVICE), - A->getColData(memory::DEVICE), - A->getValues(memory::DEVICE), - CUSPARSE_INDEX_32I, - CUSPARSE_INDEX_32I, - CUSPARSE_INDEX_BASE_ZERO, - CUDA_R_64F); - error_sum += status; - values_changed_ = false; + workspace_->resetMatvecSetup(); } + cusparseSpMatDescr_t matA = workspace_->getSpmvMatrixDescriptor(); + void* buffer_spmv = workspace_->getSpmvBuffer(); if (!workspace_->matvecSetup()) { // setup first, allocate, etc. + status = cusparseCreateCsr(&matA, + A->getNumRows(), + A->getNumColumns(), + A->getNnz(), + A->getRowData(memory::DEVICE), + A->getColData(memory::DEVICE), + A->getValues(memory::DEVICE), + CUSPARSE_INDEX_32I, + CUSPARSE_INDEX_32I, + CUSPARSE_INDEX_BASE_ZERO, + CUDA_R_64F); + error_sum += status; size_t bufferSize = 0; status = cusparseSpMV_bufferSize(handle_cusparse, @@ -123,6 +123,13 @@ namespace ReSolve workspace_->setSpmvBuffer(buffer_spmv); workspace_->matvecSetupDone(); + + matrix_for_matvec_ = A; + matvec_num_rows_ = A->getNumRows(); + matvec_num_cols_ = A->getNumColumns(); + matvec_nnz_ = A->getNnz(); + + values_changed_ = false; } status = cusparseSpMV(handle_cusparse, diff --git a/resolve/matrix/MatrixHandlerCuda.hpp b/resolve/matrix/MatrixHandlerCuda.hpp index f0efbcee..9c6d6d21 100644 --- a/resolve/matrix/MatrixHandlerCuda.hpp +++ b/resolve/matrix/MatrixHandlerCuda.hpp @@ -58,6 +58,11 @@ namespace ReSolve LinAlgWorkspaceCUDA* workspace_{nullptr}; bool values_changed_{true}; ///< needed for matvec + matrix::Sparse* matrix_for_matvec_{nullptr}; ///< matrix for cached matvec setup + index_type matvec_num_rows_{0}; ///< number of rows in matrix used for matvec + index_type matvec_num_cols_{0}; ///< number of columns in matrix used for matvec + index_type matvec_nnz_{0}; ///< number of nonzeros in matrix used for matvec + MemoryHandler mem_; ///< Device memory manager object }; diff --git a/resolve/workspace/LinAlgWorkspaceCUDA.cpp b/resolve/workspace/LinAlgWorkspaceCUDA.cpp index fdb19c30..ad492f78 100644 --- a/resolve/workspace/LinAlgWorkspaceCUDA.cpp +++ b/resolve/workspace/LinAlgWorkspaceCUDA.cpp @@ -204,6 +204,20 @@ namespace ReSolve matvec_setup_done_ = true; } + void LinAlgWorkspaceCUDA::resetMatvecSetup() + { + if (matvec_setup_done_) + { + cusparseDestroySpMat(mat_A_); + matvec_setup_done_ = false; + } + if (buffer_spmv_ != nullptr) + { + mem_.deleteOnDevice(buffer_spmv_); + buffer_spmv_ = nullptr; + } + } + void LinAlgWorkspaceCUDA::initializeHandles() { cusparseCreate(&handle_cusparse_); diff --git a/resolve/workspace/LinAlgWorkspaceCUDA.hpp b/resolve/workspace/LinAlgWorkspaceCUDA.hpp index ac996b9b..957e23fb 100644 --- a/resolve/workspace/LinAlgWorkspaceCUDA.hpp +++ b/resolve/workspace/LinAlgWorkspaceCUDA.hpp @@ -47,6 +47,7 @@ namespace ReSolve bool matvecSetup(); void matvecSetupDone(); + void resetMatvecSetup(); private: // handles From e6ad7e93b6eb3d0c5a82af73e5c960f478bfe201 Mon Sep 17 00:00:00 2001 From: tamar-dewilde Date: Mon, 25 May 2026 20:10:20 +0000 Subject: [PATCH 04/11] Apply pre-commmit fixes --- resolve/matrix/MatrixHandlerCuda.cpp | 26 +++++++++++++------------- resolve/matrix/MatrixHandlerCuda.hpp | 6 +++--- tests/unit/hykkt/HykktSCCGTests.hpp | 2 +- tests/unit/hykkt/runHykktSCCGTests.cpp | 4 ++-- 4 files changed, 19 insertions(+), 19 deletions(-) diff --git a/resolve/matrix/MatrixHandlerCuda.cpp b/resolve/matrix/MatrixHandlerCuda.cpp index 4b6b1cce..a27f2aa0 100644 --- a/resolve/matrix/MatrixHandlerCuda.cpp +++ b/resolve/matrix/MatrixHandlerCuda.cpp @@ -82,27 +82,27 @@ namespace ReSolve cusparseCreateDnVec(&vecAx, A->getNumRows(), vec_result->getData(memory::DEVICE), CUDA_R_64F); cusparseHandle_t handle_cusparse = workspace_->getCusparseHandle(); - bool matrix_changed = (matrix_for_matvec_ != A) || (matvec_num_rows_ != A->getNumRows()) || (matvec_num_cols_ != A->getNumColumns()) || (matvec_nnz_ != A->getNnz()); + bool matrix_changed = (matrix_for_matvec_ != A) || (matvec_num_rows_ != A->getNumRows()) || (matvec_num_cols_ != A->getNumColumns()) || (matvec_nnz_ != A->getNnz()); if (matrix_changed || values_changed_) { workspace_->resetMatvecSetup(); } - cusparseSpMatDescr_t matA = workspace_->getSpmvMatrixDescriptor(); - void* buffer_spmv = workspace_->getSpmvBuffer(); + cusparseSpMatDescr_t matA = workspace_->getSpmvMatrixDescriptor(); + void* buffer_spmv = workspace_->getSpmvBuffer(); if (!workspace_->matvecSetup()) { // setup first, allocate, etc. status = cusparseCreateCsr(&matA, - A->getNumRows(), - A->getNumColumns(), - A->getNnz(), - A->getRowData(memory::DEVICE), - A->getColData(memory::DEVICE), - A->getValues(memory::DEVICE), - CUSPARSE_INDEX_32I, - CUSPARSE_INDEX_32I, - CUSPARSE_INDEX_BASE_ZERO, - CUDA_R_64F); + A->getNumRows(), + A->getNumColumns(), + A->getNnz(), + A->getRowData(memory::DEVICE), + A->getColData(memory::DEVICE), + A->getValues(memory::DEVICE), + CUSPARSE_INDEX_32I, + CUSPARSE_INDEX_32I, + CUSPARSE_INDEX_BASE_ZERO, + CUDA_R_64F); error_sum += status; size_t bufferSize = 0; diff --git a/resolve/matrix/MatrixHandlerCuda.hpp b/resolve/matrix/MatrixHandlerCuda.hpp index 9c6d6d21..0ba32957 100644 --- a/resolve/matrix/MatrixHandlerCuda.hpp +++ b/resolve/matrix/MatrixHandlerCuda.hpp @@ -59,9 +59,9 @@ namespace ReSolve bool values_changed_{true}; ///< needed for matvec matrix::Sparse* matrix_for_matvec_{nullptr}; ///< matrix for cached matvec setup - index_type matvec_num_rows_{0}; ///< number of rows in matrix used for matvec - index_type matvec_num_cols_{0}; ///< number of columns in matrix used for matvec - index_type matvec_nnz_{0}; ///< number of nonzeros in matrix used for matvec + index_type matvec_num_rows_{0}; ///< number of rows in matrix used for matvec + index_type matvec_num_cols_{0}; ///< number of columns in matrix used for matvec + index_type matvec_nnz_{0}; ///< number of nonzeros in matrix used for matvec MemoryHandler mem_; ///< Device memory manager object }; diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index e1db0481..7304f9fc 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -11,8 +11,8 @@ #include #include #include -#include #include +#include #include namespace ReSolve diff --git a/tests/unit/hykkt/runHykktSCCGTests.cpp b/tests/unit/hykkt/runHykktSCCGTests.cpp index 9fc32f81..45b51518 100644 --- a/tests/unit/hykkt/runHykktSCCGTests.cpp +++ b/tests/unit/hykkt/runHykktSCCGTests.cpp @@ -33,8 +33,8 @@ void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace, WorkspaceType workspace; workspace.initializeHandles(); - ReSolve::MatrixHandler matrixHandler(&workspace); - ReSolve::VectorHandler vectorHandler(&workspace); + ReSolve::MatrixHandler matrixHandler(&workspace); + ReSolve::VectorHandler vectorHandler(&workspace); ReSolve::tests::HykktSchurComplementConjugateGradientTests test(memspace, matrixHandler, vectorHandler); result += test.SCCGTest(); From 2aa19c95dbbe3b6f3079c9b42ebc02da14850b0b Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Tue, 26 May 2026 10:08:33 -0700 Subject: [PATCH 05/11] fix HyKKT SCCG HIP target setup --- resolve/hykkt/sccg/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/hykkt/sccg/CMakeLists.txt b/resolve/hykkt/sccg/CMakeLists.txt index a89fd53b..0c28c55c 100644 --- a/resolve/hykkt/sccg/CMakeLists.txt +++ b/resolve/hykkt/sccg/CMakeLists.txt @@ -14,7 +14,7 @@ target_link_libraries( # Link to CUDA ReSolve backend if CUDA is support enabled target_include_directories( - resolve_hykkt_spgemm PUBLIC ${SUITESPARSE_INCLUDE_DIR} + resolve_hykkt_sccg PUBLIC ${SUITESPARSE_INCLUDE_DIR} ) if(RESOLVE_USE_CUDA) From 1ceb4d9ef2d169e0d81cfb8e29e0535b7183e44f Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Tue, 26 May 2026 11:51:25 -0700 Subject: [PATCH 06/11] document SCCG backend setup and CUDA SpMV cache handling --- .../sccg/SchurComplementConjugateGradient.cpp | 2 ++ .../sccg/SchurComplementConjugateGradient.hpp | 19 +++++++++++--- resolve/matrix/MatrixHandlerCuda.cpp | 1 + resolve/matrix/MatrixHandlerCuda.hpp | 10 ++++---- resolve/workspace/LinAlgWorkspaceCUDA.cpp | 4 ++- resolve/workspace/LinAlgWorkspaceCUDA.hpp | 6 +++++ tests/unit/hykkt/HykktSCCGTests.hpp | 25 +++++++++++++------ 7 files changed, 50 insertions(+), 17 deletions(-) diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index bdeacce3..2d2ac24c 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -13,6 +13,8 @@ namespace ReSolve * @param m[in] - Dimension of inner system. * @param choleskySolver[in] - Factorization of Hgamma to use for direct solve. * @param memspace[in] - Memory space of incoming data and for computation. + * @param matrixHandler[in] - Matrix handler for the selected backend. + * @param vectorHandler[in] - Vector handler for the selected backend. */ SchurComplementConjugateGradient::SchurComplementConjugateGradient( index_type n, diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp index a1b75221..94045df7 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp @@ -1,6 +1,6 @@ /** * @file SchurComplementConjugateGradient.hpp - * @brief Currently CPU-only. + * @brief Schur complement conjugate gradient solver for HyKKT. */ #pragma once @@ -12,7 +12,6 @@ #include #include #include -#include namespace ReSolve { @@ -24,6 +23,18 @@ namespace ReSolve class SchurComplementConjugateGradient { public: + /** + * @brief Constructor for SchurComplementConjugateGradient. + * + * The solver uses caller-provided matrix and vector handlers so the same solver can be run with CPU, CUDA, or HIP backends. + * + * @param[in] n Dimension of outer system. + * @param[in] m Dimension of inner system. + * @param[in] choleskySolver Factorization of Hgamma to use for direct solves. + * @param[in] memspace Memory space of incoming data and for computation. + * @param[in] matrixHandler Matrix handler for the selected backend. + * @param[in] vectorHandler Vector handler for the selected backend. + */ SchurComplementConjugateGradient(index_type n, index_type m, CholeskySolver* choleskySolver, memory::MemorySpace memspace, MatrixHandler& matrixHandler, VectorHandler& vectorHandler); void addMatrixInfo(matrix::Csr* jc, matrix::Csr* jc_tr); @@ -41,8 +52,8 @@ namespace ReSolve int itmax_ = 100; // Maximum iterations for conjugate gradient double tol_ = 1e-12; // Solver tolerance for Schur - MatrixHandler& matrixhandler_; - VectorHandler& vectorhandler_; + MatrixHandler& matrixhandler_; ///< Backend-specific matrix handler. + VectorHandler& vectorhandler_; ///< Backend-specific vector handler. CholeskySolver* choleskySolver_; // Cholesky factorization on 1,1 block diff --git a/resolve/matrix/MatrixHandlerCuda.cpp b/resolve/matrix/MatrixHandlerCuda.cpp index a27f2aa0..2427d126 100644 --- a/resolve/matrix/MatrixHandlerCuda.cpp +++ b/resolve/matrix/MatrixHandlerCuda.cpp @@ -82,6 +82,7 @@ namespace ReSolve cusparseCreateDnVec(&vecAx, A->getNumRows(), vec_result->getData(memory::DEVICE), CUDA_R_64F); cusparseHandle_t handle_cusparse = workspace_->getCusparseHandle(); + // Rebuild cached SpMV setup if the matrix object or dimensions changed. bool matrix_changed = (matrix_for_matvec_ != A) || (matvec_num_rows_ != A->getNumRows()) || (matvec_num_cols_ != A->getNumColumns()) || (matvec_nnz_ != A->getNnz()); if (matrix_changed || values_changed_) { diff --git a/resolve/matrix/MatrixHandlerCuda.hpp b/resolve/matrix/MatrixHandlerCuda.hpp index 0ba32957..13b09990 100644 --- a/resolve/matrix/MatrixHandlerCuda.hpp +++ b/resolve/matrix/MatrixHandlerCuda.hpp @@ -56,12 +56,12 @@ namespace ReSolve private: LinAlgWorkspaceCUDA* workspace_{nullptr}; - bool values_changed_{true}; ///< needed for matvec + bool values_changed_{true}; ///< Flag to indicate if matrix values changed since the cached SpMV setup. - matrix::Sparse* matrix_for_matvec_{nullptr}; ///< matrix for cached matvec setup - index_type matvec_num_rows_{0}; ///< number of rows in matrix used for matvec - index_type matvec_num_cols_{0}; ///< number of columns in matrix used for matvec - index_type matvec_nnz_{0}; ///< number of nonzeros in matrix used for matvec + matrix::Sparse* matrix_for_matvec_{nullptr}; ///< Matrix used for cached SpMV setup. + index_type matvec_num_rows_{0}; ///< Number of rows in cached SpMV matrix. + index_type matvec_num_cols_{0}; ///< Number of columns in cached SpMV matrix. + index_type matvec_nnz_{0}; ///< Number of nonzeros in cached SpMV matrix. MemoryHandler mem_; ///< Device memory manager object }; diff --git a/resolve/workspace/LinAlgWorkspaceCUDA.cpp b/resolve/workspace/LinAlgWorkspaceCUDA.cpp index ad492f78..89c7b9fb 100644 --- a/resolve/workspace/LinAlgWorkspaceCUDA.cpp +++ b/resolve/workspace/LinAlgWorkspaceCUDA.cpp @@ -203,7 +203,9 @@ namespace ReSolve { matvec_setup_done_ = true; } - + /** + * @brief Reset cached SpMV resources. + */ void LinAlgWorkspaceCUDA::resetMatvecSetup() { if (matvec_setup_done_) diff --git a/resolve/workspace/LinAlgWorkspaceCUDA.hpp b/resolve/workspace/LinAlgWorkspaceCUDA.hpp index 957e23fb..9c7b555b 100644 --- a/resolve/workspace/LinAlgWorkspaceCUDA.hpp +++ b/resolve/workspace/LinAlgWorkspaceCUDA.hpp @@ -47,6 +47,12 @@ namespace ReSolve bool matvecSetup(); void matvecSetupDone(); + /** + * @brief Reset the cached CUDA SpMV setup. + * + * Destroys the cached sparse matrix descriptor and frees the SpMV buffer so + * the next matvec call can rebuild the SpMV setup if the matrix or its dimensions have changed. + */ void resetMatvecSetup(); private: diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index 7304f9fc..44bc75f9 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -27,6 +27,15 @@ namespace ReSolve class HykktSchurComplementConjugateGradientTests : public TestBase { public: + /** + * @brief Constructs the SCCG test fixture with the specified memory space and handlers. + * + * The test fixture uses caller-provided matrix and vector handlers so the same test can be run with CPU, CUDA, or HIP backends. + * + * @param[in] memspace Memory space for the test (HOST or DEVICE). + * @param[in] matrixHandler Reference to a matrix handler for the selected backend. + * @param[in] vectorHandler Reference to a vector handler for the selected backend. + */ HykktSchurComplementConjugateGradientTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler, VectorHandler& vectorHandler) : memspace_(memspace), matrixHandler_(matrixHandler), vectorHandler_(vectorHandler) { @@ -53,6 +62,8 @@ namespace ReSolve std::ifstream hFile(hFileName); std::ifstream bFile(bFileName); + // The .mtx file readers write into host accessible memory. + // Load test data into HOST first, then sync to DEVICE for CUDA and HIP backends. matrix::Csr* h = new matrix::Csr(2278, 2278, 11304, true, false); h->allocateMatrixData(memory::HOST); io::updateMatrixFromFile(hFile, h); @@ -116,9 +127,9 @@ namespace ReSolve } private: - memory::MemorySpace memspace_; - MatrixHandler& matrixHandler_; - VectorHandler& vectorHandler_; + memory::MemorySpace memspace_; ///< Memory space used by the test. + MatrixHandler& matrixHandler_; ///< Backend-specific matrix handler. + VectorHandler& vectorHandler_; ///< Backend-specific vector handler. /** * @brief Generate a random vector of doubles between 0 and RAND_MAX. Copied from HykktCholeskyTests.hpp. @@ -138,9 +149,9 @@ namespace ReSolve } /** - * @brief Validate the results of the scaling - * @param[in] x0 Pointer to the output x0 vector - * @param[in] tol Solver tolerance + * @brief Validate the SCCG result. + * @param[in] x0 Pointer to the output x0 vector. + * @param[in] tol Solver tolerance. */ bool validateResult(vector::Vector* x0, real_type tol) { @@ -150,6 +161,6 @@ namespace ReSolve return test_passed; } - }; // class HykktPermutationTests + }; // class HykktSchurComplementConjugateGradientTests } // namespace tests } // namespace ReSolve From f8fca77cd49f80b9ba4b027720cb3e4da4075c3f Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Tue, 26 May 2026 13:06:54 -0700 Subject: [PATCH 07/11] mirror SpMV cache reset for HIP --- resolve/hykkt/sccg/CMakeLists.txt | 4 +-- resolve/matrix/MatrixHandlerCuda.cpp | 2 +- resolve/matrix/MatrixHandlerHip.cpp | 13 ++++++++++ resolve/matrix/MatrixHandlerHip.hpp | 7 +++++- resolve/workspace/LinAlgWorkspaceCUDA.cpp | 1 + resolve/workspace/LinAlgWorkspaceHIP.cpp | 30 ++++++++++++++++------- resolve/workspace/LinAlgWorkspaceHIP.hpp | 7 ++++++ 7 files changed, 50 insertions(+), 14 deletions(-) diff --git a/resolve/hykkt/sccg/CMakeLists.txt b/resolve/hykkt/sccg/CMakeLists.txt index 0c28c55c..74e3939e 100644 --- a/resolve/hykkt/sccg/CMakeLists.txt +++ b/resolve/hykkt/sccg/CMakeLists.txt @@ -13,9 +13,7 @@ target_link_libraries( # Link to CUDA ReSolve backend if CUDA is support enabled -target_include_directories( - resolve_hykkt_sccg PUBLIC ${SUITESPARSE_INCLUDE_DIR} -) +target_include_directories(resolve_hykkt_sccg PUBLIC ${SUITESPARSE_INCLUDE_DIR}) if(RESOLVE_USE_CUDA) target_sources(resolve_hykkt_sccg PRIVATE ${Matrix_CUDASDK_SRC}) diff --git a/resolve/matrix/MatrixHandlerCuda.cpp b/resolve/matrix/MatrixHandlerCuda.cpp index 2427d126..db95dbe9 100644 --- a/resolve/matrix/MatrixHandlerCuda.cpp +++ b/resolve/matrix/MatrixHandlerCuda.cpp @@ -83,7 +83,7 @@ namespace ReSolve cusparseHandle_t handle_cusparse = workspace_->getCusparseHandle(); // Rebuild cached SpMV setup if the matrix object or dimensions changed. - bool matrix_changed = (matrix_for_matvec_ != A) || (matvec_num_rows_ != A->getNumRows()) || (matvec_num_cols_ != A->getNumColumns()) || (matvec_nnz_ != A->getNnz()); + bool matrix_changed = (matrix_for_matvec_ != A) || (matvec_num_rows_ != A->getNumRows()) || (matvec_num_cols_ != A->getNumColumns()) || (matvec_nnz_ != A->getNnz()); if (matrix_changed || values_changed_) { workspace_->resetMatvecSetup(); diff --git a/resolve/matrix/MatrixHandlerHip.cpp b/resolve/matrix/MatrixHandlerHip.cpp index 8586101e..b73e7b6c 100644 --- a/resolve/matrix/MatrixHandlerHip.cpp +++ b/resolve/matrix/MatrixHandlerHip.cpp @@ -78,6 +78,13 @@ namespace ReSolve rocsparse_handle handle_rocsparse = workspace_->getRocsparseHandle(); + // Rebuild cached SpMV setup if the matrix object or dimensions changed. + bool matrix_changed = (matrix_for_matvec_ != A) || (matvec_num_rows_ != A->getNumRows()) || (matvec_num_cols_ != A->getNumColumns()) || (matvec_nnz_ != A->getNnz()); + if (matrix_changed || values_changed_) + { + workspace_->resetMatvecSetup(); + } + rocsparse_mat_info infoA = workspace_->getSpmvMatrixInfo(); rocsparse_mat_descr descrA = workspace_->getSpmvMatrixDescriptor(); @@ -106,6 +113,12 @@ namespace ReSolve workspace_->setSpmvMatrixDescriptor(descrA); workspace_->setSpmvMatrixInfo(infoA); workspace_->matvecSetupDone(); + + matrix_for_matvec_ = A; + matvec_num_rows_ = A->getNumRows(); + matvec_num_cols_ = A->getNumColumns(); + matvec_nnz_ = A->getNnz(); + values_changed_ = false; } status = rocsparse_dcsrmv(handle_rocsparse, diff --git a/resolve/matrix/MatrixHandlerHip.hpp b/resolve/matrix/MatrixHandlerHip.hpp index 7e69c3f1..04436278 100644 --- a/resolve/matrix/MatrixHandlerHip.hpp +++ b/resolve/matrix/MatrixHandlerHip.hpp @@ -56,7 +56,12 @@ namespace ReSolve private: LinAlgWorkspaceHIP* workspace_{nullptr}; - bool values_changed_{true}; ///< needed for matvec + bool values_changed_{true}; ///< Flag to indicate if matrix values changed since the cached SpMV setup. + + matrix::Sparse* matrix_for_matvec_{nullptr}; ///< Matrix used for cached SpMV setup. + index_type matvec_num_rows_{0}; ///< Number of rows in cached SpMV matrix. + index_type matvec_num_cols_{0}; ///< Number of columns in cached SpMV matrix. + index_type matvec_nnz_{0}; ///< Number of nonzeros in cached SpMV matrix. MemoryHandler mem_; ///< Device memory manager object }; diff --git a/resolve/workspace/LinAlgWorkspaceCUDA.cpp b/resolve/workspace/LinAlgWorkspaceCUDA.cpp index 89c7b9fb..1c6facee 100644 --- a/resolve/workspace/LinAlgWorkspaceCUDA.cpp +++ b/resolve/workspace/LinAlgWorkspaceCUDA.cpp @@ -203,6 +203,7 @@ namespace ReSolve { matvec_setup_done_ = true; } + /** * @brief Reset cached SpMV resources. */ diff --git a/resolve/workspace/LinAlgWorkspaceHIP.cpp b/resolve/workspace/LinAlgWorkspaceHIP.cpp index fc13e689..239dcb30 100644 --- a/resolve/workspace/LinAlgWorkspaceHIP.cpp +++ b/resolve/workspace/LinAlgWorkspaceHIP.cpp @@ -6,6 +6,8 @@ namespace ReSolve { handle_rocsparse_ = nullptr; handle_rocblas_ = nullptr; + mat_A_ = nullptr; + info_A_ = nullptr; matvec_setup_done_ = false; d_r_ = nullptr; @@ -18,12 +20,11 @@ namespace ReSolve LinAlgWorkspaceHIP::~LinAlgWorkspaceHIP() { + resetMatvecSetup(); + rocsparse_destroy_handle(handle_rocsparse_); rocblas_destroy_handle(handle_rocblas_); - if (matvec_setup_done_) - { - rocsparse_destroy_mat_descr(mat_A_); - } + if (d_r_size_ != 0) { mem_.deleteOnDevice(d_r_); @@ -47,11 +48,7 @@ namespace ReSolve */ void LinAlgWorkspaceHIP::resetLinAlgWorkspace() { - if (matvec_setup_done_) - { - rocsparse_destroy_mat_descr(mat_A_); - matvec_setup_done_ = false; - } + resetMatvecSetup(); if (d_r_size_ != 0) { mem_.deleteOnDevice(d_r_); @@ -143,6 +140,21 @@ namespace ReSolve matvec_setup_done_ = true; } + void LinAlgWorkspaceHIP::resetMatvecSetup() + { + if (mat_A_ != nullptr) + { + rocsparse_destroy_mat_descr(mat_A_); + mat_A_ = nullptr; + } + if (info_A_ != nullptr) + { + rocsparse_destroy_mat_info(info_A_); + info_A_ = nullptr; + } + matvec_setup_done_ = false; + } + void LinAlgWorkspaceHIP::initializeHandles() { rocsparse_create_handle(&handle_rocsparse_); diff --git a/resolve/workspace/LinAlgWorkspaceHIP.hpp b/resolve/workspace/LinAlgWorkspaceHIP.hpp index 4db4b89f..91171257 100644 --- a/resolve/workspace/LinAlgWorkspaceHIP.hpp +++ b/resolve/workspace/LinAlgWorkspaceHIP.hpp @@ -38,6 +38,13 @@ namespace ReSolve bool matvecSetup(); void matvecSetupDone(); + /** + * @brief Reset the cached HIP SpMV setup. + * + * Destroys the cached rocSPARSE matrix descriptor and matrix info so the + * next matvec call can rebuild the setup if the matrix or its dimensions have changed. + */ + void resetMatvecSetup(); void setDrSize(index_type new_sz); void setDr(real_type* new_dr); From 16cefc1fbfb1b14e6e1ea39f0e1dd4d3b6bb0160 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Wed, 27 May 2026 11:18:22 -0700 Subject: [PATCH 08/11] clean up build warnings --- resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp | 7 +++---- resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp | 1 - tests/unit/hykkt/HykktSCCGTests.hpp | 2 +- 3 files changed, 4 insertions(+), 6 deletions(-) diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index 2d2ac24c..bb5a23dc 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -25,17 +25,16 @@ namespace ReSolve VectorHandler& vectorHandler) : n_(n), m_(m), - choleskySolver_(choleskySolver), - memspace_(memspace), - mem_(), matrixhandler_(matrixHandler), vectorhandler_(vectorHandler), + choleskySolver_(choleskySolver), y_(m_), z_(m_), r_(n_), p_(n_), s_(n_), - w_(n_) + w_(n_), + memspace_(memspace) { ; } diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp index 94045df7..a42292ba 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp @@ -79,7 +79,6 @@ namespace ReSolve vector::Vector s_; vector::Vector w_; - MemoryHandler mem_; memory::MemorySpace memspace_; }; // class SchurComplementConjugateGradient } // namespace hykkt diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index 44bc75f9..b6d0a271 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -153,7 +153,7 @@ namespace ReSolve * @param[in] x0 Pointer to the output x0 vector. * @param[in] tol Solver tolerance. */ - bool validateResult(vector::Vector* x0, real_type tol) + bool validateResult(vector::Vector*, real_type) { bool test_passed = true; From b884811cad27bd7c25d01f8a0ffc4aa03c7291c0 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Wed, 27 May 2026 12:45:59 -0700 Subject: [PATCH 09/11] update variable names to follow conventions --- .../sccg/SchurComplementConjugateGradient.cpp | 44 +++++++++---------- .../sccg/SchurComplementConjugateGradient.hpp | 10 ++--- tests/unit/hykkt/HykktSCCGTests.hpp | 36 +++++++-------- tests/unit/hykkt/runHykktSCCGTests.cpp | 6 +-- 4 files changed, 48 insertions(+), 48 deletions(-) diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index bb5a23dc..41d5141d 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -13,20 +13,20 @@ namespace ReSolve * @param m[in] - Dimension of inner system. * @param choleskySolver[in] - Factorization of Hgamma to use for direct solve. * @param memspace[in] - Memory space of incoming data and for computation. - * @param matrixHandler[in] - Matrix handler for the selected backend. - * @param vectorHandler[in] - Vector handler for the selected backend. + * @param matrix_handler[in] - Matrix handler for the selected backend. + * @param vector_handler[in] - Vector handler for the selected backend. */ SchurComplementConjugateGradient::SchurComplementConjugateGradient( index_type n, index_type m, CholeskySolver* choleskySolver, memory::MemorySpace memspace, - MatrixHandler& matrixHandler, - VectorHandler& vectorHandler) + MatrixHandler& matrix_handler, + VectorHandler& vector_handler) : n_(n), m_(m), - matrixhandler_(matrixHandler), - vectorhandler_(vectorHandler), + matrix_handler_(matrix_handler), + vector_handler_(vector_handler), choleskySolver_(choleskySolver), y_(m_), z_(m_), @@ -103,37 +103,37 @@ namespace ReSolve { using namespace constants; - matrixhandler_.matvec(jc_tr_, x0_, &y_, &ONE, &ZERO, memspace_); + matrix_handler_.matvec(jc_tr_, x0_, &y_, &ONE, &ZERO, memspace_); choleskySolver_->solve(&z_, &y_); - matrixhandler_.matvec(jc_, &z_, &r_, &MINUS_ONE, &ONE, memspace_); - gam_i_ = vectorhandler_.dot(&r_, &r_, memspace_); + matrix_handler_.matvec(jc_, &z_, &r_, &MINUS_ONE, &ONE, memspace_); + gam_i_ = vector_handler_.dot(&r_, &r_, memspace_); - matrixhandler_.matvec(jc_tr_, &r_, &y_, &ONE, &ZERO, memspace_); + matrix_handler_.matvec(jc_tr_, &r_, &y_, &ONE, &ZERO, memspace_); choleskySolver_->solve(&z_, &y_); - matrixhandler_.matvec(jc_, &z_, &w_, &ONE, &ZERO, memspace_); - delta_ = vectorhandler_.dot(&w_, &r_, memspace_); + matrix_handler_.matvec(jc_, &z_, &w_, &ONE, &ZERO, memspace_); + delta_ = vector_handler_.dot(&w_, &r_, memspace_); alpha_ = gam_i_ / delta_; int i; for (i = 0; i < itmax_; i++) { - vectorhandler_.scal(beta_, &p_, memspace_); - vectorhandler_.axpy(ONE, &r_, &p_, memspace_); - vectorhandler_.scal(beta_, &s_, memspace_); - vectorhandler_.axpy(ONE, &w_, &s_, memspace_); - vectorhandler_.axpy(alpha_, &p_, x0_, memspace_); + vector_handler_.scal(beta_, &p_, memspace_); + vector_handler_.axpy(ONE, &r_, &p_, memspace_); + vector_handler_.scal(beta_, &s_, memspace_); + vector_handler_.axpy(ONE, &w_, &s_, memspace_); + vector_handler_.axpy(alpha_, &p_, x0_, memspace_); minalpha_ = -alpha_; - vectorhandler_.axpy(minalpha_, &s_, &r_, memspace_); - gam_i1_ = vectorhandler_.dot(&r_, &r_, memspace_); + vector_handler_.axpy(minalpha_, &s_, &r_, memspace_); + gam_i1_ = vector_handler_.dot(&r_, &r_, memspace_); if (sqrt(gam_i1_) < tol_) { printf("Convergence occured at iteration %d\n", i); break; } - matrixhandler_.matvec(jc_tr_, &r_, &y_, &ONE, &ZERO, memspace_); + matrix_handler_.matvec(jc_tr_, &r_, &y_, &ONE, &ZERO, memspace_); choleskySolver_->solve(&z_, &y_); - matrixhandler_.matvec(jc_, &z_, &w_, &ONE, &ZERO, memspace_); - delta_ = vectorhandler_.dot(&w_, &r_, memspace_); + matrix_handler_.matvec(jc_, &z_, &w_, &ONE, &ZERO, memspace_); + delta_ = vector_handler_.dot(&w_, &r_, memspace_); beta_ = gam_i1_ / gam_i_; gam_i_ = gam_i1_; alpha_ = gam_i_ / (delta_ - beta_ * gam_i_ / alpha_); diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp index a42292ba..9f2131e5 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp @@ -32,10 +32,10 @@ namespace ReSolve * @param[in] m Dimension of inner system. * @param[in] choleskySolver Factorization of Hgamma to use for direct solves. * @param[in] memspace Memory space of incoming data and for computation. - * @param[in] matrixHandler Matrix handler for the selected backend. - * @param[in] vectorHandler Vector handler for the selected backend. + * @param[in] matrix_handler Matrix handler for the selected backend. + * @param[in] vector_handler Vector handler for the selected backend. */ - SchurComplementConjugateGradient(index_type n, index_type m, CholeskySolver* choleskySolver, memory::MemorySpace memspace, MatrixHandler& matrixHandler, VectorHandler& vectorHandler); + SchurComplementConjugateGradient(index_type n, index_type m, CholeskySolver* choleskySolver, memory::MemorySpace memspace, MatrixHandler& matrix_handler, VectorHandler& vector_handler); void addMatrixInfo(matrix::Csr* jc, matrix::Csr* jc_tr); void addVectorInfo(vector::Vector* x0, vector::Vector* b); @@ -52,8 +52,8 @@ namespace ReSolve int itmax_ = 100; // Maximum iterations for conjugate gradient double tol_ = 1e-12; // Solver tolerance for Schur - MatrixHandler& matrixhandler_; ///< Backend-specific matrix handler. - VectorHandler& vectorhandler_; ///< Backend-specific vector handler. + MatrixHandler& matrix_handler_; ///< Backend-specific matrix handler. + VectorHandler& vector_handler_; ///< Backend-specific vector handler. CholeskySolver* choleskySolver_; // Cholesky factorization on 1,1 block diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index b6d0a271..32eb6ab2 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -33,11 +33,11 @@ namespace ReSolve * The test fixture uses caller-provided matrix and vector handlers so the same test can be run with CPU, CUDA, or HIP backends. * * @param[in] memspace Memory space for the test (HOST or DEVICE). - * @param[in] matrixHandler Reference to a matrix handler for the selected backend. - * @param[in] vectorHandler Reference to a vector handler for the selected backend. + * @param[in] matrix_handler Reference to a matrix handler for the selected backend. + * @param[in] vector_handler Reference to a vector handler for the selected backend. */ - HykktSchurComplementConjugateGradientTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler, VectorHandler& vectorHandler) - : memspace_(memspace), matrixHandler_(matrixHandler), vectorHandler_(vectorHandler) + HykktSchurComplementConjugateGradientTests(memory::MemorySpace memspace, MatrixHandler& matrix_handler, VectorHandler& vector_handler) + : memspace_(memspace), matrix_handler_(matrix_handler), vector_handler_(vector_handler) { } @@ -54,19 +54,19 @@ namespace ReSolve { constexpr double tol = 1e-12; - std::string sourceDir = std::string(SOURCE_DIR); - std::string jcFileName = sourceDir + "/SCCGTestMatrices/JC_matrix_ACTIVSg200_AC_00.mtx"; - std::string hFileName = sourceDir + "/SCCGTestMatrices/H_matrix_ACTIVSg200_AC_00.mtx"; - std::string bFileName = sourceDir + "/SCCGTestMatrices/CG_rhs_ACTIVSg200_AC_00.mtx"; // rhs - std::ifstream jcFile(jcFileName); - std::ifstream hFile(hFileName); - std::ifstream bFile(bFileName); + std::string source_dir = std::string(SOURCE_DIR); + std::string jc_file_name = source_dir + "/SCCGTestMatrices/JC_matrix_ACTIVSg200_AC_00.mtx"; + std::string h_file_name = source_dir + "/SCCGTestMatrices/H_matrix_ACTIVSg200_AC_00.mtx"; + std::string b_file_name = source_dir + "/SCCGTestMatrices/CG_rhs_ACTIVSg200_AC_00.mtx"; // rhs + std::ifstream jc_file(jc_file_name); + std::ifstream h_file(h_file_name); + std::ifstream b_file(b_file_name); // The .mtx file readers write into host accessible memory. // Load test data into HOST first, then sync to DEVICE for CUDA and HIP backends. matrix::Csr* h = new matrix::Csr(2278, 2278, 11304, true, false); h->allocateMatrixData(memory::HOST); - io::updateMatrixFromFile(hFile, h); + io::updateMatrixFromFile(h_file, h); if (memspace_ == memory::DEVICE) { h->syncData(memory::DEVICE); @@ -79,7 +79,7 @@ namespace ReSolve matrix::Csr* jc = new matrix::Csr(1386, 2278, 6784, false, false); jc->allocateMatrixData(memory::HOST); - io::updateMatrixFromFile(jcFile, jc); + io::updateMatrixFromFile(jc_file, jc); if (memspace_ == memory::DEVICE) { jc->syncData(memory::DEVICE); @@ -88,12 +88,12 @@ namespace ReSolve index_type n = jc->getNumRows(); index_type m = jc->getNumColumns(); index_type nnz = jc->getNnz(); - hykkt::SchurComplementConjugateGradient sccg(n, m, &choleskySolver, memspace_, matrixHandler_, vectorHandler_); + hykkt::SchurComplementConjugateGradient sccg(n, m, &choleskySolver, memspace_, matrix_handler_, vector_handler_); sccg.setSolverTolerance(tol); matrix::Csr* jc_tr = new matrix::Csr(m, n, nnz); jc_tr->allocateMatrixData(memspace_); - matrixHandler_.transpose(jc, jc_tr, memspace_); + matrix_handler_.transpose(jc, jc_tr, memspace_); vector::Vector* x0 = new vector::Vector(n); x0->allocate(memory::HOST); @@ -101,7 +101,7 @@ namespace ReSolve vector::Vector* b = new vector::Vector(n); b->allocate(memory::HOST); - io::updateVectorFromFile(bFile, b); + io::updateVectorFromFile(b_file, b); if (memspace_ == memory::DEVICE) { b->syncData(memory::DEVICE); @@ -128,8 +128,8 @@ namespace ReSolve private: memory::MemorySpace memspace_; ///< Memory space used by the test. - MatrixHandler& matrixHandler_; ///< Backend-specific matrix handler. - VectorHandler& vectorHandler_; ///< Backend-specific vector handler. + MatrixHandler& matrix_handler_; ///< Backend-specific matrix handler. + VectorHandler& vector_handler_; ///< Backend-specific vector handler. /** * @brief Generate a random vector of doubles between 0 and RAND_MAX. Copied from HykktCholeskyTests.hpp. diff --git a/tests/unit/hykkt/runHykktSCCGTests.cpp b/tests/unit/hykkt/runHykktSCCGTests.cpp index 45b51518..b3462676 100644 --- a/tests/unit/hykkt/runHykktSCCGTests.cpp +++ b/tests/unit/hykkt/runHykktSCCGTests.cpp @@ -33,9 +33,9 @@ void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace, WorkspaceType workspace; workspace.initializeHandles(); - ReSolve::MatrixHandler matrixHandler(&workspace); - ReSolve::VectorHandler vectorHandler(&workspace); - ReSolve::tests::HykktSchurComplementConjugateGradientTests test(memspace, matrixHandler, vectorHandler); + ReSolve::MatrixHandler matrix_handler(&workspace); + ReSolve::VectorHandler vector_handler(&workspace); + ReSolve::tests::HykktSchurComplementConjugateGradientTests test(memspace, matrix_handler, vector_handler); result += test.SCCGTest(); From 6a54f81047517f0c383da78d852df167503a63ff Mon Sep 17 00:00:00 2001 From: tamar-dewilde Date: Wed, 27 May 2026 19:46:46 +0000 Subject: [PATCH 10/11] Apply pre-commmit fixes --- tests/unit/hykkt/HykktSCCGTests.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index 32eb6ab2..2d853474 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -54,7 +54,7 @@ namespace ReSolve { constexpr double tol = 1e-12; - std::string source_dir = std::string(SOURCE_DIR); + std::string source_dir = std::string(SOURCE_DIR); std::string jc_file_name = source_dir + "/SCCGTestMatrices/JC_matrix_ACTIVSg200_AC_00.mtx"; std::string h_file_name = source_dir + "/SCCGTestMatrices/H_matrix_ACTIVSg200_AC_00.mtx"; std::string b_file_name = source_dir + "/SCCGTestMatrices/CG_rhs_ACTIVSg200_AC_00.mtx"; // rhs @@ -127,7 +127,7 @@ namespace ReSolve } private: - memory::MemorySpace memspace_; ///< Memory space used by the test. + memory::MemorySpace memspace_; ///< Memory space used by the test. MatrixHandler& matrix_handler_; ///< Backend-specific matrix handler. VectorHandler& vector_handler_; ///< Backend-specific vector handler. From e76a84560f5644ddd9ecfd976f2bd7b2b13d967e Mon Sep 17 00:00:00 2001 From: Shaked Regev <35384901+shakedregev@users.noreply.github.com> Date: Thu, 28 May 2026 09:32:05 -0400 Subject: [PATCH 11/11] Shaked/random error fix (#424) Fixed stochastic failures due to a faulty random number generator and asynchronous operations. Co-authored-by: shakedregev --- CHANGELOG.md | 4 ++++ resolve/Common.hpp | 11 ++++++----- resolve/hykkt/cholesky/CholeskySolverHip.cpp | 3 +++ tests/unit/hykkt/HykktCholeskyTests.hpp | 19 ++++++++++++------- tests/unit/hykkt/HykktSCCGTests.hpp | 12 ++++++++---- tests/unit/hykkt/runHykktCholeskyTests.cpp | 8 +++++--- tests/unit/hykkt/runHykktSCCGTests.cpp | 4 +++- 7 files changed, 41 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3c760b73..fdeb1ade 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,10 @@ - Added classes and tests for permutation, Ruiz scaling, Cholesky factorization, Schur complement conjugate gradient and matrix multiplication and addition. +- Changed random number generation int tests to be C++ style and fixed-seed, to avoid random failures. + +- Added Schur Complement Conjugate Gradient class. + ## Changes to Re::Solve since release 0.99.2 - Added cmake-format. diff --git a/resolve/Common.hpp b/resolve/Common.hpp index af30adad..5d5ae88a 100644 --- a/resolve/Common.hpp +++ b/resolve/Common.hpp @@ -14,11 +14,12 @@ namespace ReSolve namespace constants { - constexpr real_type ZERO = 0.0; - constexpr real_type ONE = 1.0; - constexpr real_type TWO = 2.0; - constexpr real_type HALF = 0.5; - constexpr real_type MINUS_ONE = -1.0; + constexpr real_type ZERO = 0.0; + constexpr real_type ONE = 1.0; + constexpr real_type TWO = 2.0; + constexpr real_type HALF = 0.5; + constexpr real_type MINUS_ONE = -1.0; + constexpr index_type SEED = 12345; constexpr real_type MACHINE_EPSILON = std::numeric_limits::epsilon(); } // namespace constants diff --git a/resolve/hykkt/cholesky/CholeskySolverHip.cpp b/resolve/hykkt/cholesky/CholeskySolverHip.cpp index cc8ae4fa..358ebea9 100644 --- a/resolve/hykkt/cholesky/CholeskySolverHip.cpp +++ b/resolve/hykkt/cholesky/CholeskySolverHip.cpp @@ -60,6 +60,7 @@ namespace ReSolve } A_chol_ = convertToCholmod(A); A_ = A; + mem_.deviceSynchronize(); } void CholeskySolverHip::symbolicAnalysis() @@ -70,6 +71,7 @@ namespace ReSolve { out::error() << "Cholesky symbolic analysis failed with status: " << Common_.status << "\n"; } + mem_.deviceSynchronize(); } /** @@ -156,6 +158,7 @@ namespace ReSolve out::error() << "Refactorization step failed with status: " << status << "\n"; } L_->setUpdated(memory::DEVICE); + mem_.deviceSynchronize(); } } diff --git a/tests/unit/hykkt/HykktCholeskyTests.hpp b/tests/unit/hykkt/HykktCholeskyTests.hpp index ec214381..cf1a9b8a 100644 --- a/tests/unit/hykkt/HykktCholeskyTests.hpp +++ b/tests/unit/hykkt/HykktCholeskyTests.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -26,8 +27,8 @@ namespace ReSolve class HykktCholeskyTests : public TestBase { public: - HykktCholeskyTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler) - : memspace_(memspace), matrixHandler_(matrixHandler) + HykktCholeskyTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler, std::mt19937& generator) + : memspace_(memspace), matrixHandler_(matrixHandler), generator_(generator) { cholmod_start(&Common); } @@ -127,7 +128,6 @@ namespace ReSolve 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); @@ -242,7 +242,8 @@ namespace ReSolve // reset values for (size_t j = 0; j < L->nzmax; j++) { - static_cast(L->x)[j] = 2.0 * static_cast(rand()) / RAND_MAX - 1.0; + 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); @@ -266,6 +267,7 @@ namespace ReSolve private: ReSolve::memory::MemorySpace memspace_; MatrixHandler& matrixHandler_; + std::mt19937& generator_; cholmod_common Common; @@ -281,14 +283,16 @@ namespace ReSolve L_p[i + 1] = L_p[i]; for (size_t j = i; j < n; ++j) { - // with probability 'density', add a non-zero entry - if (i == j || static_cast(rand()) / RAND_MAX < density) + 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); @@ -311,9 +315,10 @@ namespace ReSolve { vector::Vector* v = new vector::Vector(n); v->allocate(memory::HOST); + std::uniform_real_distribution distribution(0.0, 1.0); for (index_type i = 0; i < n; ++i) { - v->getData(memory::HOST)[i] = static_cast(rand()) / RAND_MAX; + v->getData(memory::HOST)[i] = distribution(generator_); } v->setDataUpdated(memory::HOST); if (memspace_ == memory::DEVICE) diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index 2d853474..035c82cd 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -6,6 +6,7 @@ #pragma once #include +#include #include #include @@ -35,9 +36,10 @@ namespace ReSolve * @param[in] memspace Memory space for the test (HOST or DEVICE). * @param[in] matrix_handler Reference to a matrix handler for the selected backend. * @param[in] vector_handler Reference to a vector handler for the selected backend. + * @param[in] generator Reference to a C++ random number generator. */ - HykktSchurComplementConjugateGradientTests(memory::MemorySpace memspace, MatrixHandler& matrix_handler, VectorHandler& vector_handler) - : memspace_(memspace), matrix_handler_(matrix_handler), vector_handler_(vector_handler) + HykktSchurComplementConjugateGradientTests(memory::MemorySpace memspace, MatrixHandler& matrix_handler, VectorHandler& vector_handler, std::mt19937& generator) + : memspace_(memspace), matrix_handler_(matrix_handler), vector_handler_(vector_handler), generator_(generator) { } @@ -130,16 +132,18 @@ namespace ReSolve memory::MemorySpace memspace_; ///< Memory space used by the test. MatrixHandler& matrix_handler_; ///< Backend-specific matrix handler. VectorHandler& vector_handler_; ///< Backend-specific vector handler. + std::mt19937& generator_; ///< C++ random number generator. /** - * @brief Generate a random vector of doubles between 0 and RAND_MAX. Copied from HykktCholeskyTests.hpp. + * @brief Generate a random vector of doubles between 0 and 1. Copied from HykktCholeskyTests.hpp. * @param[in] vec Target vector to write to. */ void randomVector(vector::Vector* vec) { + std::uniform_real_distribution distribution(0.0, 1.0); for (index_type i = 0; i < vec->getSize(); ++i) { - vec->getData(memory::HOST)[i] = static_cast(rand()) / RAND_MAX; + vec->getData(memory::HOST)[i] = distribution(generator_); } vec->setDataUpdated(memory::HOST); if (memspace_ == memory::DEVICE) diff --git a/tests/unit/hykkt/runHykktCholeskyTests.cpp b/tests/unit/hykkt/runHykktCholeskyTests.cpp index d05fb560..55ebeb7d 100644 --- a/tests/unit/hykkt/runHykktCholeskyTests.cpp +++ b/tests/unit/hykkt/runHykktCholeskyTests.cpp @@ -7,8 +7,10 @@ */ #include #include +#include #include +#include "resolve/Common.hpp" #include "tests/unit/hykkt/HykktCholeskyTests.hpp" /** @@ -24,9 +26,9 @@ void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace, WorkspaceType workspace; workspace.initializeHandles(); - ReSolve::MatrixHandler handler(&workspace); - - ReSolve::tests::HykktCholeskyTests test(memspace, handler); + ReSolve::MatrixHandler handler(&workspace); + std::mt19937 generator(ReSolve::constants::SEED); // set random seed for reproducibility + ReSolve::tests::HykktCholeskyTests test(memspace, handler, generator); result += test.minimalCorrectness(); handler.setValuesChanged(true, memspace); diff --git a/tests/unit/hykkt/runHykktSCCGTests.cpp b/tests/unit/hykkt/runHykktSCCGTests.cpp index b3462676..33f0e385 100644 --- a/tests/unit/hykkt/runHykktSCCGTests.cpp +++ b/tests/unit/hykkt/runHykktSCCGTests.cpp @@ -5,6 +5,7 @@ */ #include #include +#include #include #include @@ -35,7 +36,8 @@ void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace, workspace.initializeHandles(); ReSolve::MatrixHandler matrix_handler(&workspace); ReSolve::VectorHandler vector_handler(&workspace); - ReSolve::tests::HykktSchurComplementConjugateGradientTests test(memspace, matrix_handler, vector_handler); + std::mt19937 generator(ReSolve::constants::SEED); + ReSolve::tests::HykktSchurComplementConjugateGradientTests test(memspace, matrix_handler, vector_handler, generator); result += test.SCCGTest();