diff --git a/CHANGELOG.md b/CHANGELOG.md index 3c760b738..fdeb1ade5 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 af30adadb..5d5ae88aa 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 cc8ae4faa..358ebea9d 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/resolve/hykkt/sccg/CMakeLists.txt b/resolve/hykkt/sccg/CMakeLists.txt index a89fd53b6..74e3939e8 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_spgemm 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/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index 3e81d116e..41d5141d6 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -13,26 +13,28 @@ 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 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) + memory::MemorySpace memspace, + MatrixHandler& matrix_handler, + VectorHandler& vector_handler) : n_(n), m_(m), + matrix_handler_(matrix_handler), + vector_handler_(vector_handler), choleskySolver_(choleskySolver), - memspace_(memspace), - mem_(), - workspace_(), - matrixhandler_(&workspace_), - vectorhandler_(&workspace_), y_(m_), z_(m_), r_(n_), p_(n_), s_(n_), - w_(n_) + w_(n_), + memspace_(memspace) { ; } @@ -101,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 88664c6e7..9f2131e59 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,7 +23,19 @@ namespace ReSolve class SchurComplementConjugateGradient { public: - SchurComplementConjugateGradient(index_type n, index_type m, CholeskySolver* choleskySolver, memory::MemorySpace memspace); + /** + * @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] 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& matrix_handler, VectorHandler& vector_handler); void addMatrixInfo(matrix::Csr* jc, matrix::Csr* jc_tr); void addVectorInfo(vector::Vector* x0, vector::Vector* b); @@ -41,15 +52,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& matrix_handler_; ///< Backend-specific matrix handler. + VectorHandler& vector_handler_; ///< Backend-specific vector handler. CholeskySolver* choleskySolver_; // Cholesky factorization on 1,1 block @@ -75,7 +79,6 @@ namespace ReSolve vector::Vector s_; vector::Vector w_; - MemoryHandler mem_; memory::MemorySpace memspace_; }; // class SchurComplementConjugateGradient } // namespace hykkt diff --git a/resolve/matrix/MatrixHandlerCuda.cpp b/resolve/matrix/MatrixHandlerCuda.cpp index 7b3108438..db95dbe9b 100644 --- a/resolve/matrix/MatrixHandlerCuda.cpp +++ b/resolve/matrix/MatrixHandlerCuda.cpp @@ -76,17 +76,23 @@ 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_) + // 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(); + } + cusparseSpMatDescr_t matA = workspace_->getSpmvMatrixDescriptor(); + void* buffer_spmv = workspace_->getSpmvBuffer(); + if (!workspace_->matvecSetup()) { + // setup first, allocate, etc. status = cusparseCreateCsr(&matA, A->getNumRows(), A->getNumColumns(), @@ -99,11 +105,6 @@ namespace ReSolve CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F); error_sum += status; - values_changed_ = false; - } - if (!workspace_->matvecSetup()) - { - // setup first, allocate, etc. size_t bufferSize = 0; status = cusparseSpMV_bufferSize(handle_cusparse, @@ -123,6 +124,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 f0efbcee8..13b099903 100644 --- a/resolve/matrix/MatrixHandlerCuda.hpp +++ b/resolve/matrix/MatrixHandlerCuda.hpp @@ -56,7 +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 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/matrix/MatrixHandlerHip.cpp b/resolve/matrix/MatrixHandlerHip.cpp index 8586101e7..b73e7b6c3 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 7e69c3f1c..044362784 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 fdb19c303..1c6facee0 100644 --- a/resolve/workspace/LinAlgWorkspaceCUDA.cpp +++ b/resolve/workspace/LinAlgWorkspaceCUDA.cpp @@ -204,6 +204,23 @@ namespace ReSolve matvec_setup_done_ = true; } + /** + * @brief Reset cached SpMV resources. + */ + 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 ac996b9bb..9c7b555be 100644 --- a/resolve/workspace/LinAlgWorkspaceCUDA.hpp +++ b/resolve/workspace/LinAlgWorkspaceCUDA.hpp @@ -47,6 +47,13 @@ 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: // handles diff --git a/resolve/workspace/LinAlgWorkspaceHIP.cpp b/resolve/workspace/LinAlgWorkspaceHIP.cpp index fc13e6896..239dcb305 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 4db4b89fa..911712579 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); diff --git a/tests/unit/hykkt/HykktCholeskyTests.hpp b/tests/unit/hykkt/HykktCholeskyTests.hpp index ec214381c..cf1a9b8a1 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 c26cbc34a..035c82cd2 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -6,12 +6,14 @@ #pragma once #include +#include #include #include #include #include #include +#include #include namespace ReSolve @@ -26,8 +28,18 @@ namespace ReSolve class HykktSchurComplementConjugateGradientTests : public TestBase { public: - HykktSchurComplementConjugateGradientTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler) - : memspace_(memspace), matrixHandler_(matrixHandler) + /** + * @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] 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, std::mt19937& generator) + : memspace_(memspace), matrix_handler_(matrix_handler), vector_handler_(vector_handler), generator_(generator) { } @@ -44,17 +56,23 @@ 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(memspace_); - io::updateMatrixFromFile(hFile, h); + h->allocateMatrixData(memory::HOST); + io::updateMatrixFromFile(h_file, h); + if (memspace_ == memory::DEVICE) + { + h->syncData(memory::DEVICE); + } hykkt::CholeskySolver choleskySolver(memspace_); choleskySolver.addMatrixInfo(h); choleskySolver.symbolicAnalysis(); @@ -62,28 +80,34 @@ namespace ReSolve choleskySolver.numericalFactorization(); matrix::Csr* jc = new matrix::Csr(1386, 2278, 6784, false, false); - jc->allocateMatrixData(memspace_); - io::updateMatrixFromFile(jcFile, jc); + jc->allocateMatrixData(memory::HOST); + io::updateMatrixFromFile(jc_file, jc); + if (memspace_ == memory::DEVICE) + { + jc->syncData(memory::DEVICE); + } 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_, matrix_handler_, vector_handler_); 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_); + matrix_handler_.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_); - io::updateVectorFromFile(bFile, b); + b->allocate(memory::HOST); + io::updateVectorFromFile(b_file, b); + if (memspace_ == memory::DEVICE) + { + b->syncData(memory::DEVICE); + } sccg.addMatrixInfo(jc, jc_tr); sccg.addVectorInfo(x0, b); @@ -105,18 +129,21 @@ namespace ReSolve } private: - memory::MemorySpace memspace_; - MatrixHandler& matrixHandler_; + 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) @@ -126,11 +153,11 @@ 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) + bool validateResult(vector::Vector*, real_type) { bool test_passed = true; @@ -138,6 +165,6 @@ namespace ReSolve return test_passed; } - }; // class HykktPermutationTests + }; // class HykktSchurComplementConjugateGradientTests } // namespace tests } // namespace ReSolve diff --git a/tests/unit/hykkt/runHykktCholeskyTests.cpp b/tests/unit/hykkt/runHykktCholeskyTests.cpp index d05fb5603..55ebeb7d6 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 c4a052467..33f0e3856 100644 --- a/tests/unit/hykkt/runHykktSCCGTests.cpp +++ b/tests/unit/hykkt/runHykktSCCGTests.cpp @@ -5,9 +5,11 @@ */ #include #include +#include #include #include +#include #include #ifdef RESOLVE_USE_CUDA #include @@ -32,9 +34,10 @@ 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 matrix_handler(&workspace); + ReSolve::VectorHandler vector_handler(&workspace); + std::mt19937 generator(ReSolve::constants::SEED); + ReSolve::tests::HykktSchurComplementConjugateGradientTests test(memspace, matrix_handler, vector_handler, generator); result += test.SCCGTest();