diff --git a/CHANGELOG.md b/CHANGELOG.md index 77086b150..ad188d2e5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## HyKKT Release changes +- Exported HyKKT libraries and fixed solver reuse by refreshing numerical data, reusing allocations, resizing GPU transpose workspaces, and handling zero residuals. - 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. diff --git a/resolve/CMakeLists.txt b/resolve/CMakeLists.txt index e64085db3..19eb51430 100644 --- a/resolve/CMakeLists.txt +++ b/resolve/CMakeLists.txt @@ -157,6 +157,16 @@ endif() # Add HyKKT solver if(RESOLVE_USE_KLU) add_subdirectory(hykkt) + list( + APPEND + ReSolve_Targets_List + resolve_hykkt + resolve_hykkt_ruiz + resolve_hykkt_chol + resolve_hykkt_spgemm + resolve_hykkt_sccg + resolve_hykkt_solver + ) endif() # Set installable targets diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 1a37685cc..d83a2fc22 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -7,10 +7,12 @@ #include "HyKKTSolver.hpp" #include +#include namespace ReSolve { using namespace constants; + using out = io::Logger; /** * @brief basic constructor @@ -62,6 +64,9 @@ namespace ReSolve * values. It will only set pointers to user provided data; it is user's * responsibility to supply and later delete that memory. * + * Reusing the solver requires the sparsity patterns of the matrix blocks + * to remain unchanged. Changing J_d between empty and nonempty is rejected. + * * @param[in] H_plus_D_x - Pointer to the Hessian matrix block (n_x x n_x), * corresponding to H + D_x in the HyKKT paper. * @param[in] D_s - Pointer to the slack variables derivatives matrix block @@ -70,18 +75,32 @@ namespace ReSolve * (m_c x n_x). * @param[in] J_d - Pointer to the inequality constraints Jacobian block * (m_d x n_x) + * + * @return 0 if successful, 1 if the supplied blocks are incompatible with + * the allocated solver state. On failure, the previously stored + * blocks are left unchanged. */ - void hykkt::HyKKTSolver::setMatrixBlocks(matrix::Csr* H_plus_D_x, matrix::Csr* D_s, matrix::Csr* J, matrix::Csr* J_d) + int hykkt::HyKKTSolver::setMatrixBlocks(matrix::Csr* H_plus_D_x, matrix::Csr* D_s, matrix::Csr* J, matrix::Csr* J_d) { + const bool J_d_flag = J_d->getNnz() > 0; + + if (allocated_ && J_d_flag_ != J_d_flag) + { + out::error() << "Changing J_d between empty and nonempty is not " + "supported when reusing HyKKT.\n"; + return 1; + } + H_ = H_plus_D_x; D_s_ = D_s; J_ = J; J_d_ = J_d; - bool J_d_flag = J_d->getNnz() > 0; - // status_ = (J_d_flag_ == J_d_flag); - status_ = true; // when using API, we can't check if sparsity pattern changed - J_d_flag_ = J_d_flag; + if (!allocated_) + { + J_d_flag_ = J_d_flag; + } + return 0; } /** @@ -156,16 +175,6 @@ namespace ReSolve */ real_type hykkt::HyKKTSolver::solve() { - // TODO: Review sparsity pattern checking in HyKKT - if (!status_ && allocated_) - { - printf("\n\nERROR: USING HYKKT WITH NEW NONZERO STRUCTURE\n\n"); - std::cout << "status = " << status_ - << ", allocated = " << allocated_ - << "\n"; - return 1; - } - setupParameters(); if (!allocated_) @@ -200,6 +209,16 @@ namespace ReSolve } computeHgammaFactorization(); + if (!allocated_) + { + sccg_ = new SchurComplementConjugateGradient(J_->getNumRows(), + J_->getNumColumns(), + cholesky_, + matrixHandler_, + vectorHandler_, + memspace_); + sccg_->setup(); + } setupConjugateGradient(); computeConjugateGradient(); @@ -210,8 +229,8 @@ namespace ReSolve /** * @brief allocates and initiates variables for KKT system * - * @pre jd_flag_ determines if variables used for Spgemm H_tilde - * should be initiated + * @pre J_d_flag_ determines whether variables used to form H_tilde with + * SpGEMM should be initialized. * * @post all variables used for hykkt are allocated for; J_d- * related variables are not initiated if J_d nnz == 0 @@ -242,7 +261,6 @@ namespace ReSolve J_tr_perm_ = new matrix::Csr(J_tr_->getNumRows(), J_tr_->getNumColumns(), J_tr_->getNnz()); J_d_scaled_ = new matrix::Csr(J_d_->getNumRows(), J_d_->getNumColumns(), J_d_->getNnz()); - D_s_vals_->setData(D_s_->getValues(memspace_), memspace_); r_yd_scaled_->allocate(memspace_); r_x_perm_->allocate(memspace_); omega_perm_->allocate(memspace_); @@ -258,14 +276,12 @@ namespace ReSolve J_d_scaled_->allocateWithExternalSparsityPattern(J_d_->getRowData(memspace_), J_d_->getColData(memspace_), J_d_->getNnz(), memspace_); // H_tilde_ does not need to be allocated because loadResultMatrix() does it later } - else if (memspace_ == memory::DEVICE) - { - J_d_->syncData(memory::DEVICE); // check if this is redundant - } + // D_s may be replaced between solves, so refresh the external value pointer. + D_s_vals_->setData(D_s_->getValues(memspace_), memspace_); r_y_copy_->copyFromExternal(r_y_, memspace_, memspace_); - // check if this is redundant in later iterations + // Matrix values may change between solves, so refresh the transpose. matrixHandler_->transpose(J_, J_tr_, memspace_); if (J_d_flag_) { @@ -314,7 +330,10 @@ namespace ReSolve else { H_tilde_->setNnz(H_->getNnz()); - H_tilde_->allocateMatrixData(memspace_); + if (!allocated_) + { + H_tilde_->allocateMatrixData(memspace_); + } H_tilde_->copyFromExternal(H_->getRowData(memspace_), H_->getColData(memspace_), H_->getValues(memspace_), @@ -388,9 +407,6 @@ namespace ReSolve void hykkt::HyKKTSolver::setupSpGEMMHgamma() { spgemm_hgamma_ = new SpGEMM(memspace_, gamma_, ONE); - spgemm_hgamma_->loadProductMatrices(J_tr_, J_); - spgemm_hgamma_->loadSumMatrix(H_tilde_); - spgemm_hgamma_->loadResultMatrix(&H_gamma_); // H_gamma_ will be created when calling SpGEMM->compute() } /* @@ -403,6 +419,14 @@ namespace ReSolve */ void hykkt::HyKKTSolver::computeSpGEMMHgamma() { + // Numerical values can change between solves while the sparsity pattern + // remains fixed, so refresh the SpGEMM inputs before recomputing H_gamma. + spgemm_hgamma_->setCoefficients(gamma_, ONE); + spgemm_hgamma_->loadProductMatrices(J_tr_, J_); + spgemm_hgamma_->loadSumMatrix(H_tilde_); + // HIP initializes the result descriptor using dimensions established by + // the product and sum inputs, so load the result matrix after both inputs. + spgemm_hgamma_->loadResultMatrix(&H_gamma_); spgemm_hgamma_->compute(); r_x_hat_->copyFromExternal(r_x_til_, memspace_, memspace_); matrixHandler_->matvec(J_tr_, r_y_, r_x_hat_, &gamma_, &ONE, memspace_); @@ -514,11 +538,9 @@ namespace ReSolve schur_->copyFromExternal(r_y_, memspace_, memspace_); matrixHandler_->matvec(J_perm_, omega_perm_, schur_, &ONE, &MINUS_ONE, memspace_); - sccg_ = new SchurComplementConjugateGradient(J_->getNumRows(), J_->getNumColumns(), cholesky_, matrixHandler_, vectorHandler_, memspace_); sccg_->addMatrixInfo(J_perm_, J_tr_perm_); y_->setToZero(memspace_); sccg_->addVectorInfo(y_, schur_); - sccg_->setup(); } /** @@ -553,10 +575,6 @@ namespace ReSolve // block-recovering the solution to the original system by parts // this part is to recover delta_x cholesky_->solve(z_, r_x_perm_); - if (memspace_ == memory::DEVICE) - { - x_->syncData(memory::DEVICE); - } permutation_->mapIndex(REV_PERM_V, z_->getData(memspace_), x_->getData(memspace_)); x_->setDataUpdated(memspace_); @@ -627,10 +645,17 @@ namespace ReSolve matrixHandler_->matvec(J_copy_, x_, r_y_copy_, &MINUS_ONE, &ONE, memspace_); norm_resy_sq = vectorHandler_->dot(r_y_copy_, r_y_copy_, memspace_); - // Calculate final relative norm norm_resx_sq += norm_resy_sq; - real_type norm_res = sqrt(norm_resx_sq) / sqrt(norm_r_x_sq); - printf("||Ax-b||/||b|| = %32.32g\n\n", norm_res); + real_type norm_res = sqrt(norm_resx_sq); + if (norm_r_x_sq > 0) + { + norm_res /= sqrt(norm_r_x_sq); + printf("||Ax-b||/||b|| = %32.32g\n\n", norm_res); + } + else + { + printf("||Ax-b|| = %32.32g\n\n", norm_res); + } allocated_ = true; diff --git a/resolve/hykkt/HyKKTSolver.hpp b/resolve/hykkt/HyKKTSolver.hpp index a3129d137..a1c93a07e 100644 --- a/resolve/hykkt/HyKKTSolver.hpp +++ b/resolve/hykkt/HyKKTSolver.hpp @@ -41,7 +41,7 @@ namespace ReSolve std::istream& r_y_file, std::istream& r_yd_file); - void setMatrixBlocks(matrix::Csr* H_plus_D_x, matrix::Csr* D_s, matrix::Csr* J, matrix::Csr* J_d); + int setMatrixBlocks(matrix::Csr* H_plus_D_x, matrix::Csr* D_s, matrix::Csr* J, matrix::Csr* J_d); void setRHSBlocks(vector::Vector* r_x, vector::Vector* r_s, vector::Vector* r_y, vector::Vector* r_yd); void setLHSPointers(vector::Vector* x, vector::Vector* s, vector::Vector* y, vector::Vector* y_d); @@ -78,10 +78,6 @@ namespace ReSolve bool allocated_ = false; bool J_d_flag_ = false; - // Whether the solver is correctly used with matrices of - // the same nonzero structure - bool status_ = true; - RuizScaling* ruiz_{nullptr}; SpGEMM* spgemm_htil_{nullptr}; SpGEMM* spgemm_hgamma_{nullptr}; diff --git a/resolve/hykkt/cholesky/CholeskySolverCpu.cpp b/resolve/hykkt/cholesky/CholeskySolverCpu.cpp index 296585341..1cae503c8 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCpu.cpp +++ b/resolve/hykkt/cholesky/CholeskySolverCpu.cpp @@ -39,11 +39,12 @@ namespace ReSolve void CholeskySolverCpu::addMatrixInfo(matrix::Csr* A) { + A_ = A; if (A_chol_) { cholmod_free_sparse(&A_chol_, &Common_); } - A_chol_ = convertToCholmod(A); + A_chol_ = convertToCholmod(A_); } /** @@ -72,6 +73,12 @@ namespace ReSolve { (void) tol; // Mark tol as unused + // CHOLMOD stores its own copy of the matrix values. Refresh that copy + // before numerical refactorization when the solver is reused. + mem_.copyArrayHostToHost(static_cast(A_chol_->x), + A_->getValues(memory::HOST), + A_->getNnz()); + cholmod_factorize(A_chol_, factorization_, &Common_); if (Common_.status < 0) { diff --git a/resolve/hykkt/cholesky/CholeskySolverCpu.hpp b/resolve/hykkt/cholesky/CholeskySolverCpu.hpp index fcfda9e31..60ea0b69f 100644 --- a/resolve/hykkt/cholesky/CholeskySolverCpu.hpp +++ b/resolve/hykkt/cholesky/CholeskySolverCpu.hpp @@ -27,6 +27,7 @@ namespace ReSolve MemoryHandler mem_; cholmod_common Common_; + matrix::Csr* A_ = nullptr; cholmod_sparse* A_chol_; // cholmod sparse matrix representation cholmod_factor* factorization_; diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index 8fff342da..791f8cd9f 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -99,6 +99,11 @@ namespace ReSolve p_->allocate(memspace_); s_->allocate(memspace_); w_->allocate(memspace_); + } + + int SchurComplementConjugateGradient::solve() + { + using namespace constants; y_->setToZero(memspace_); z_->setToZero(memspace_); @@ -110,17 +115,18 @@ namespace ReSolve x_0_->setToZero(memspace_); beta_ = 0; - } - - int SchurComplementConjugateGradient::solve() - { - using namespace constants; matrix_handler_->matvec(J_tr_, x_0_, y_, &ONE, &ZERO, memspace_); choleskySolver_->solve(z_, y_); matrix_handler_->matvec(J_, z_, r_, &MINUS_ONE, &ONE, memspace_); gamma_i_ = vector_handler_->dot(r_, r_, memspace_); + if (sqrt(gamma_i_) < tol_) + { + gamma_i1_ = gamma_i_; + return 0; + } + matrix_handler_->matvec(J_tr_, r_, y_, &ONE, &ZERO, memspace_); choleskySolver_->solve(z_, y_); matrix_handler_->matvec(J_, z_, w_, &ONE, &ZERO, memspace_); diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp index 8ba032bb7..bf0700983 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp @@ -49,6 +49,11 @@ namespace ReSolve void setSolverTolerance(double tol); void setSolverItmax(int itmax); + /** + * @brief Allocates internal work vectors. + * + * @pre Must be called exactly once before the first call to solve(). + */ void setup(); int solve(); diff --git a/resolve/hykkt/spgemm/SpGEMM.cpp b/resolve/hykkt/spgemm/SpGEMM.cpp index 1170f919f..2610dcdd1 100644 --- a/resolve/hykkt/spgemm/SpGEMM.cpp +++ b/resolve/hykkt/spgemm/SpGEMM.cpp @@ -53,6 +53,17 @@ namespace ReSolve delete impl_; } + /** + * Updates the coefficients for the SpGEMM operation. + * + * @param[in] alpha - Scalar multiplier for the matrix product. + * @param[in] beta - Scalar multiplier for the sum matrix. + */ + void SpGEMM::setCoefficients(real_type alpha, real_type beta) + { + impl_->setCoefficients(alpha, beta); + } + /** * Loads the two matrices for the product * @param A[in] - Pointer to CSR matrix diff --git a/resolve/hykkt/spgemm/SpGEMM.hpp b/resolve/hykkt/spgemm/SpGEMM.hpp index edc3a267a..21a3390e5 100644 --- a/resolve/hykkt/spgemm/SpGEMM.hpp +++ b/resolve/hykkt/spgemm/SpGEMM.hpp @@ -22,6 +22,8 @@ namespace ReSolve SpGEMM(memory::MemorySpace memspace, real_type alpha, real_type beta); ~SpGEMM(); + void setCoefficients(real_type alpha, real_type beta); + void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); void loadSumMatrix(matrix::Csr* D); void loadResultMatrix(matrix::Csr** E_ptr); diff --git a/resolve/hykkt/spgemm/SpGEMMCpu.cpp b/resolve/hykkt/spgemm/SpGEMMCpu.cpp index f634231d8..93ad651b3 100644 --- a/resolve/hykkt/spgemm/SpGEMMCpu.cpp +++ b/resolve/hykkt/spgemm/SpGEMMCpu.cpp @@ -47,6 +47,12 @@ namespace ReSolve cholmod_finish(&Common_); } + void SpGEMMCpu::setCoefficients(real_type alpha, real_type beta) + { + alpha_ = alpha; + beta_ = beta; + } + void SpGEMMCpu::loadProductMatrices(matrix::Csr* A, matrix::Csr* B) { if (!A_) diff --git a/resolve/hykkt/spgemm/SpGEMMCpu.hpp b/resolve/hykkt/spgemm/SpGEMMCpu.hpp index 7c61c9ede..3522be0ca 100644 --- a/resolve/hykkt/spgemm/SpGEMMCpu.hpp +++ b/resolve/hykkt/spgemm/SpGEMMCpu.hpp @@ -16,11 +16,13 @@ namespace ReSolve SpGEMMCpu(real_type alpha, real_type beta); ~SpGEMMCpu(); - void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); - void loadSumMatrix(matrix::Csr* D); - void loadResultMatrix(matrix::Csr** E_ptr); + void setCoefficients(real_type alpha, real_type beta) override; - void compute(); + void loadProductMatrices(matrix::Csr* A, matrix::Csr* B) override; + void loadSumMatrix(matrix::Csr* D) override; + void loadResultMatrix(matrix::Csr** E_ptr) override; + + void compute() override; private: real_type alpha_; diff --git a/resolve/hykkt/spgemm/SpGEMMCuda.cpp b/resolve/hykkt/spgemm/SpGEMMCuda.cpp index 734d2473c..0d4310a9b 100644 --- a/resolve/hykkt/spgemm/SpGEMMCuda.cpp +++ b/resolve/hykkt/spgemm/SpGEMMCuda.cpp @@ -29,6 +29,12 @@ namespace ReSolve cusparseDestroy(handle_); } + void SpGEMMCuda::setCoefficients(real_type alpha, real_type beta) + { + alpha_ = alpha; + beta_ = beta; + } + void SpGEMMCuda::loadProductMatrices(matrix::Csr* A, matrix::Csr* B) { A_descr_ = convertToCusparseType(A); @@ -149,7 +155,7 @@ namespace ReSolve mem_.deleteOnDevice(temp_buffer_2); int64_t C_num_cols = 0; - int64_t C_nnz_ = 0; + C_nnz_ = 0; cusparseSpMatGetSize(C_descr_, &n_, &C_num_cols, &C_nnz_); mem_.allocateArrayOnDevice(&C_col_ind_, (index_type) C_nnz_); diff --git a/resolve/hykkt/spgemm/SpGEMMCuda.hpp b/resolve/hykkt/spgemm/SpGEMMCuda.hpp index a39b2cdc9..8bd22fc99 100644 --- a/resolve/hykkt/spgemm/SpGEMMCuda.hpp +++ b/resolve/hykkt/spgemm/SpGEMMCuda.hpp @@ -22,6 +22,8 @@ namespace ReSolve SpGEMMCuda(real_type alpha, real_type beta); ~SpGEMMCuda(); + void setCoefficients(real_type alpha, real_type beta) override; + void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); void loadSumMatrix(matrix::Csr* D); void loadResultMatrix(matrix::Csr** E_ptr); diff --git a/resolve/hykkt/spgemm/SpGEMMHip.cpp b/resolve/hykkt/spgemm/SpGEMMHip.cpp index 38fafc928..fc2c62bc2 100644 --- a/resolve/hykkt/spgemm/SpGEMMHip.cpp +++ b/resolve/hykkt/spgemm/SpGEMMHip.cpp @@ -30,6 +30,12 @@ namespace ReSolve mem_.deleteOnDevice(buffer_); } + void SpGEMMHip::setCoefficients(real_type alpha, real_type beta) + { + alpha_ = alpha; + beta_ = beta; + } + void SpGEMMHip::loadProductMatrices(matrix::Csr* A, matrix::Csr* B) { if (A_descr_) diff --git a/resolve/hykkt/spgemm/SpGEMMHip.hpp b/resolve/hykkt/spgemm/SpGEMMHip.hpp index 00627cd32..51423db9d 100644 --- a/resolve/hykkt/spgemm/SpGEMMHip.hpp +++ b/resolve/hykkt/spgemm/SpGEMMHip.hpp @@ -22,11 +22,13 @@ namespace ReSolve SpGEMMHip(real_type alpha, real_type beta); ~SpGEMMHip(); - void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); - void loadSumMatrix(matrix::Csr* D); - void loadResultMatrix(matrix::Csr** E_ptr); + void setCoefficients(real_type alpha, real_type beta) override; - void compute(); + void loadProductMatrices(matrix::Csr* A, matrix::Csr* B) override; + void loadSumMatrix(matrix::Csr* D) override; + void loadResultMatrix(matrix::Csr** E_ptr) override; + + void compute() override; private: MemoryHandler mem_; diff --git a/resolve/hykkt/spgemm/SpGEMMImpl.hpp b/resolve/hykkt/spgemm/SpGEMMImpl.hpp index ee8575dbf..e70f0a0bd 100644 --- a/resolve/hykkt/spgemm/SpGEMMImpl.hpp +++ b/resolve/hykkt/spgemm/SpGEMMImpl.hpp @@ -23,6 +23,8 @@ namespace ReSolve SpGEMMImpl() = default; virtual ~SpGEMMImpl() = default; + virtual void setCoefficients(real_type alpha, real_type beta) = 0; + virtual void loadProductMatrices(matrix::Csr* A, matrix::Csr* B) = 0; virtual void loadSumMatrix(matrix::Csr* D) = 0; virtual void loadResultMatrix(matrix::Csr** E_ptr) = 0; diff --git a/resolve/matrix/MatrixHandlerCuda.cpp b/resolve/matrix/MatrixHandlerCuda.cpp index 92890dc4e..a98a46a4d 100644 --- a/resolve/matrix/MatrixHandlerCuda.cpp +++ b/resolve/matrix/MatrixHandlerCuda.cpp @@ -306,30 +306,25 @@ namespace ReSolve index_type n = A->getNumColumns(); index_type nnz = A->getNnz(); cusparseStatus_t status; - bool allocated = workspace_->isTransposeBufferAllocated(); - // check dimensions of A and At - if (!allocated) - { - // allocate transpose buffer - size_t bufferSize; - status = cusparseCsr2cscEx2_bufferSize(workspace_->getCusparseHandle(), - m, - n, - nnz, - A->getValues(memory::DEVICE), - A->getRowData(memory::DEVICE), - A->getColData(memory::DEVICE), - At->getValues(memory::DEVICE), - At->getRowData(memory::DEVICE), - At->getColData(memory::DEVICE), - CUDA_R_64F, - CUSPARSE_ACTION_NUMERIC, - CUSPARSE_INDEX_BASE_ZERO, - CUSPARSE_CSR2CSC_ALG1, - &bufferSize); - error_sum += status; - workspace_->setTransposeBufferWorkspace(bufferSize); - } + // Ensure the shared transpose workspace is large enough for this matrix. + size_t bufferSize; + status = cusparseCsr2cscEx2_bufferSize(workspace_->getCusparseHandle(), + m, + n, + nnz, + A->getValues(memory::DEVICE), + A->getRowData(memory::DEVICE), + A->getColData(memory::DEVICE), + At->getValues(memory::DEVICE), + At->getRowData(memory::DEVICE), + At->getColData(memory::DEVICE), + CUDA_R_64F, + CUSPARSE_ACTION_NUMERIC, + CUSPARSE_INDEX_BASE_ZERO, + CUSPARSE_CSR2CSC_ALG1, + &bufferSize); + error_sum += status; + error_sum += workspace_->setTransposeBufferWorkspace(bufferSize); status = cusparseCsr2cscEx2(workspace_->getCusparseHandle(), m, n, diff --git a/resolve/matrix/MatrixHandlerHip.cpp b/resolve/matrix/MatrixHandlerHip.cpp index eabe14817..56abbbf6d 100644 --- a/resolve/matrix/MatrixHandlerHip.cpp +++ b/resolve/matrix/MatrixHandlerHip.cpp @@ -281,22 +281,18 @@ namespace ReSolve index_type n = A->getNumColumns(); index_type nnz = A->getNnz(); rocsparse_status status; - bool allocated = workspace_->isTransposeBufferAllocated(); - if (!allocated) - { - // allocate transpose buffer - size_t bufferSize; - status = rocsparse_csr2csc_buffer_size(workspace_->getRocsparseHandle(), - m, - n, - nnz, - A->getRowData(memory::DEVICE), - A->getColData(memory::DEVICE), - rocsparse_action_numeric, - &bufferSize); - error_sum += status; - workspace_->setTransposeBufferWorkspace(bufferSize); - } + // Ensure the shared transpose workspace is large enough for this matrix. + size_t bufferSize; + status = rocsparse_csr2csc_buffer_size(workspace_->getRocsparseHandle(), + m, + n, + nnz, + A->getRowData(memory::DEVICE), + A->getColData(memory::DEVICE), + rocsparse_action_numeric, + &bufferSize); + error_sum += status; + error_sum += workspace_->setTransposeBufferWorkspace(bufferSize); status = rocsparse_dcsr2csc(workspace_->getRocsparseHandle(), m, n, diff --git a/resolve/workspace/LinAlgWorkspaceCUDA.cpp b/resolve/workspace/LinAlgWorkspaceCUDA.cpp index 98bb934ac..3ac6a7915 100644 --- a/resolve/workspace/LinAlgWorkspaceCUDA.cpp +++ b/resolve/workspace/LinAlgWorkspaceCUDA.cpp @@ -95,12 +95,19 @@ namespace ReSolve int LinAlgWorkspaceCUDA::setTransposeBufferWorkspace(size_t bufferSize) { + if (transpose_workspace_ready_ && bufferSize <= transpose_workspace_size_) + { + return 0; + } + if (transpose_workspace_ready_) { - out::error() << "Transpose workspace already set!\n"; - return 1; + mem_.deleteOnDevice(transpose_workspace_); + transpose_workspace_ = nullptr; } + mem_.allocateBufferOnDevice(&transpose_workspace_, bufferSize); + transpose_workspace_size_ = bufferSize; transpose_workspace_ready_ = true; return 0; } diff --git a/resolve/workspace/LinAlgWorkspaceCUDA.hpp b/resolve/workspace/LinAlgWorkspaceCUDA.hpp index 512f173b2..544348c3a 100644 --- a/resolve/workspace/LinAlgWorkspaceCUDA.hpp +++ b/resolve/workspace/LinAlgWorkspaceCUDA.hpp @@ -74,8 +74,9 @@ namespace ReSolve bool matvec_setup_done_{false}; // check if setup is done for matvec i.e. if buffer is allocated, csr structure is set etc. - void* transpose_workspace_{nullptr}; // needed for transpose - bool transpose_workspace_ready_{false}; // to track if allocated + void* transpose_workspace_{nullptr}; // needed for transpose + size_t transpose_workspace_size_{0}; // allocated size in bytes + bool transpose_workspace_ready_{false}; // to track if allocated real_type* d_r_{nullptr}; // needed for one-norm index_type d_r_size_{0}; diff --git a/resolve/workspace/LinAlgWorkspaceHIP.cpp b/resolve/workspace/LinAlgWorkspaceHIP.cpp index 4ef861374..5092fc720 100644 --- a/resolve/workspace/LinAlgWorkspaceHIP.cpp +++ b/resolve/workspace/LinAlgWorkspaceHIP.cpp @@ -193,12 +193,19 @@ namespace ReSolve int LinAlgWorkspaceHIP::setTransposeBufferWorkspace(size_t bufferSize) { + if (transpose_workspace_ready_ && bufferSize <= transpose_workspace_size_) + { + return 0; + } + if (transpose_workspace_ready_) { - out::error() << "Transpose workspace already set!\n"; - return 1; + mem_.deleteOnDevice(transpose_workspace_); + transpose_workspace_ = nullptr; } + mem_.allocateBufferOnDevice(&transpose_workspace_, bufferSize); + transpose_workspace_size_ = bufferSize; transpose_workspace_ready_ = true; return 0; } diff --git a/resolve/workspace/LinAlgWorkspaceHIP.hpp b/resolve/workspace/LinAlgWorkspaceHIP.hpp index db0815204..7f8eed986 100644 --- a/resolve/workspace/LinAlgWorkspaceHIP.hpp +++ b/resolve/workspace/LinAlgWorkspaceHIP.hpp @@ -71,6 +71,7 @@ namespace ReSolve real_type* d_r_{nullptr}; // needed for inf-norm real_type* norm_buffer_{nullptr}; // needed for inf-norm void* transpose_workspace_{nullptr}; // needed for transpose + size_t transpose_workspace_size_{0}; // allocated size in bytes bool transpose_workspace_ready_{false}; // to track if allocated index_type d_r_size_{0}; bool norm_buffer_ready_{false}; // to track if allocated diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index 4a3fe5331..f3958b310 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -104,6 +104,15 @@ namespace ReSolve testname += " n=" + std::to_string(n) + ", m=" + std::to_string(m) + ", nnz =" + std::to_string(nnz); status *= validateResult(x_0, converged_n); + // A zero initial residual is already converged and must not enter + // conjugate-gradient divisions with zero numerator and denominator. + x_0->setToZero(memspace_); + b->setToZero(memspace_); + sccg.addVectorInfo(x_0, b); + int zero_residual_converged_n = sccg.solve(); + status *= (zero_residual_converged_n == 0); + status *= (vector_handler_.dot(x_0, x_0, memspace_) <= sccg_tol); + delete H; delete J; delete J_tr; diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index 9d2a9d947..b71232a91 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -142,6 +142,86 @@ namespace ReSolve testname += " N=" + std::to_string(N) + ", nnz =" + std::to_string(nnz) + '\n'; status *= validateResult(error, tol); + // Replace D_s to exercise its pointer refresh; restore other inputs + // in place. + std::ifstream D_s_reuse_file(D_s_file_name); + std::ifstream J_update_file(J_file_name); + std::ifstream r_x_update_file(r_x_file_name); + std::ifstream r_s_update_file(r_s_file_name); + std::ifstream r_y_update_file(r_y_file_name); + std::ifstream r_yd_update_file(r_yd_file_name); + + matrix::Csr* D_s_reuse = io::createCsrFromFile(D_s_reuse_file, false); + + io::updateMatrixFromFile(J_update_file, J); + io::updateVectorFromFile(r_x_update_file, r_x); + io::updateVectorFromFile(r_s_update_file, r_s); + io::updateVectorFromFile(r_y_update_file, r_y); + io::updateVectorFromFile(r_yd_update_file, r_yd); + + real_type* D_s_values = D_s_reuse->getValues(memory::HOST); + for (index_type i = 0; i < D_s_reuse->getNnz(); ++i) + { + D_s_values[i] *= 1.1; + } + + D_s_reuse->setUpdated(memory::HOST); + J->setUpdated(memory::HOST); + + if (memspace_ == memory::DEVICE) + { + int restore_status = 0; + restore_status |= D_s_reuse->allocateMatrixData(memory::DEVICE); + restore_status |= D_s_reuse->syncData(memory::DEVICE); + restore_status |= J->syncData(memory::DEVICE); + restore_status |= r_x->syncData(memory::DEVICE); + restore_status |= r_s->syncData(memory::DEVICE); + restore_status |= r_y->syncData(memory::DEVICE); + restore_status |= r_yd->syncData(memory::DEVICE); + status *= (restore_status == 0); + } + + hykktSolver.setMatrixBlocks(H, D_s_reuse, J, J_d); + + // Change gamma to verify the cached SpGEMM coefficient is refreshed. + hykktSolver.setGamma(gamma * 1.1); + real_type second_error = hykktSolver.solve(); + status *= validateResult(second_error, tol); + + // Check that a zero RHS doesn't result in NaNs. + r_x->setToZero(memspace_); + r_s->setToZero(memspace_); + r_y->setToZero(memspace_); + r_yd->setToZero(memspace_); + real_type zero_rhs_error = hykktSolver.solve(); + status *= validateResult(zero_rhs_error, tol); + + // Check that the solver raises an error when trying to change J_d + // from nonempty to empty. + matrix::Csr* J_d_empty = new matrix::Csr(J_d->getNumRows(), + J_d->getNumColumns(), + 0); + + int structure_status = hykktSolver.setMatrixBlocks(H, D_s_reuse, J, J_d_empty); + status *= (structure_status != 0); + + // Exercise initialization and reuse with an empty J_d. + hykkt::HyKKTSolver no_jd_solver(n_x, m_d, m_c, memspace_); + no_jd_solver.setMatrixBlocks(H, D_s_reuse, J, J_d_empty); + no_jd_solver.setRHSBlocks(r_x, r_s, r_y, r_yd); + no_jd_solver.setLHSPointers(x, s, y, y_d); + no_jd_solver.setGamma(gamma); + no_jd_solver.addHandlers(&matrixHandler_, &vectorHandler_); + + real_type no_jd_error = no_jd_solver.solve(); + status *= validateResult(no_jd_error, tol); + + real_type no_jd_reuse_error = no_jd_solver.solve(); + status *= validateResult(no_jd_reuse_error, tol); + + delete D_s_reuse; + delete J_d_empty; + delete H; delete D_s; delete J; diff --git a/tests/unit/hykkt/HykktSpGEMMTests.hpp b/tests/unit/hykkt/HykktSpGEMMTests.hpp index ecaeb987a..24fd945c2 100644 --- a/tests/unit/hykkt/HykktSpGEMMTests.hpp +++ b/tests/unit/hykkt/HykktSpGEMMTests.hpp @@ -292,6 +292,58 @@ namespace ReSolve status *= verifyResult(E, 2.0); + // Recompute with updated product and sum coefficients. + spgemm.setCoefficients(3.0, 3.0); + spgemm.compute(); + + if (memspace_ == memory::DEVICE) + { + E->syncData(memory::HOST); + } + + status *= verifyResult(E, 3.0); + + // Change only beta and verify that the sum-matrix contribution changes. + real_type* equal_coefficient_values = new real_type[E->getNnz()]; + for (index_type i = 0; i < E->getNnz(); ++i) + { + equal_coefficient_values[i] = E->getValues(memory::HOST)[i]; + } + + spgemm.setCoefficients(3.0, 5.0); + spgemm.compute(); + + if (memspace_ == memory::DEVICE) + { + E->syncData(memory::HOST); + } + + bool beta_changed_result = false; + for (index_type i = 0; i < E->getNnz(); ++i) + { + if (fabs(E->getValues(memory::HOST)[i] - equal_coefficient_values[i]) > 1e-12) + { + beta_changed_result = true; + break; + } + } + + if (!beta_changed_result) + { + std::cerr << "Changing beta did not change the SpGEMM result.\n"; + } + status *= beta_changed_result; + delete[] equal_coefficient_values; + + // Restore equal coefficients for the matrix-value reuse check below. + spgemm.setCoefficients(3.0, 3.0); + spgemm.compute(); + + if (memspace_ == memory::DEVICE) + { + E->syncData(memory::HOST); + } + for (index_type j = 0; j < A->getNnz(); j++) { A->getValues(memory::HOST)[j] *= 2.0; @@ -318,7 +370,7 @@ namespace ReSolve E->syncData(memory::HOST); } - status *= verifyResult(E, 4.0); + status *= verifyResult(E, 6.0); delete A; delete B; diff --git a/tests/unit/matrix/runMatrixHandlerTests.cpp b/tests/unit/matrix/runMatrixHandlerTests.cpp index 84e9a6011..03e0f72fa 100644 --- a/tests/unit/matrix/runMatrixHandlerTests.cpp +++ b/tests/unit/matrix/runMatrixHandlerTests.cpp @@ -46,6 +46,9 @@ void runTests(const std::string& backend, ReSolve::tests::TestingResults& result result += test.csc2csr(1200, 1024); workspace.resetLinAlgWorkspace(); result += test.transpose(3, 3); + // Reuse the same workspace for a much larger transpose. CUDA and HIP + // should grow the buffer instead of reusing the first allocation. + result += test.transpose(2048, 1024); workspace.resetLinAlgWorkspace(); result += test.transpose(5, 3); workspace.resetLinAlgWorkspace(); @@ -59,8 +62,6 @@ void runTests(const std::string& backend, ReSolve::tests::TestingResults& result workspace.resetLinAlgWorkspace(); result += test.transpose(1024, 2048); workspace.resetLinAlgWorkspace(); - result += test.transpose(2048, 1024); - workspace.resetLinAlgWorkspace(); result += test.transpose(1024, 1200); workspace.resetLinAlgWorkspace(); result += test.transpose(1200, 1024);