From 189e6e02f723b438eb14d2298e3b3554502fd6c4 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 15:36:57 -0700 Subject: [PATCH 01/22] Export HyKKT libraries --- resolve/CMakeLists.txt | 10 ++++++++++ 1 file changed, 10 insertions(+) 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 From 5bc050b9a46cef4bcfeb82275e9795c2c1dd6955 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 15:46:51 -0700 Subject: [PATCH 02/22] Fix HyKKT numerical updates on solver reuse --- resolve/hykkt/HyKKTSolver.cpp | 4 ++ resolve/hykkt/cholesky/CholeskySolverCpu.cpp | 9 +++- resolve/hykkt/cholesky/CholeskySolverCpu.hpp | 1 + tests/unit/hykkt/HykktSolverTests.hpp | 53 ++++++++++++++++++++ 4 files changed, 66 insertions(+), 1 deletion(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 1a37685cc..40e1f83cc 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -403,6 +403,10 @@ 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_->loadProductMatrices(J_tr_, J_); + spgemm_hgamma_->loadSumMatrix(H_tilde_); spgemm_hgamma_->compute(); r_x_hat_->copyFromExternal(r_x_til_, memspace_, memspace_); matrixHandler_->matvec(J_tr_, r_y_, r_x_hat_, &gamma_, &ONE, memspace_); 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/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index 9d2a9d947..e197f96fa 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -142,6 +142,59 @@ namespace ReSolve testname += " N=" + std::to_string(N) + ", nnz =" + std::to_string(nnz) + '\n'; status *= validateResult(error, tol); + // Update D_s and restore data modified by the first solve. + real_type* D_s_values = D_s->getValues(memory::HOST); + for (index_type i = 0; i < D_s->getNnz(); ++i) + { + D_s_values[i] *= 1.1; + } + D_s->setUpdated(memory::HOST); + + std::ifstream J_reuse_file(J_file_name); + std::ifstream r_x_reuse_file(r_x_file_name); + std::ifstream r_s_reuse_file(r_s_file_name); + std::ifstream r_y_reuse_file(r_y_file_name); + std::ifstream r_yd_reuse_file(r_yd_file_name); + + matrix::Csr* J_reuse = io::createCsrFromFile(J_reuse_file, false); + + vector::Vector* r_x_reuse = io::createVectorFromFile(r_x_reuse_file); + vector::Vector* r_s_reuse = io::createVectorFromFile(r_s_reuse_file); + vector::Vector* r_y_reuse = io::createVectorFromFile(r_y_reuse_file); + vector::Vector* r_yd_reuse = io::createVectorFromFile(r_yd_reuse_file); + + int restore_status = 0; + restore_status |= J->copyFromExternal(J_reuse->getRowData(memory::HOST), + J_reuse->getColData(memory::HOST), + J_reuse->getValues(memory::HOST), + memory::HOST, + memory::HOST); + restore_status |= r_x->copyFromExternal(r_x_reuse, memory::HOST, memory::HOST); + restore_status |= r_s->copyFromExternal(r_s_reuse, memory::HOST, memory::HOST); + restore_status |= r_y->copyFromExternal(r_y_reuse, memory::HOST, memory::HOST); + restore_status |= r_yd->copyFromExternal(r_yd_reuse, memory::HOST, memory::HOST); + + if (memspace_ == memory::DEVICE) + { + restore_status |= D_s->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); + + delete J_reuse; + delete r_x_reuse; + delete r_s_reuse; + delete r_y_reuse; + delete r_yd_reuse; + + real_type second_error = hykktSolver.solve(); + status *= validateResult(second_error, tol); + delete H; delete D_s; delete J; From a4eccb257bec783c5b85fe5e8a2d6532f5b50109 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 15:49:21 -0700 Subject: [PATCH 03/22] Fix CUDA SpGEMM nonzero count state --- resolve/hykkt/spgemm/SpGEMMCuda.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/hykkt/spgemm/SpGEMMCuda.cpp b/resolve/hykkt/spgemm/SpGEMMCuda.cpp index 734d2473c..a9f93b60e 100644 --- a/resolve/hykkt/spgemm/SpGEMMCuda.cpp +++ b/resolve/hykkt/spgemm/SpGEMMCuda.cpp @@ -149,7 +149,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_); From 8bf63dc1f08df266fa273cd336baebf7e7602e40 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 15:50:30 -0700 Subject: [PATCH 04/22] Resize GPU transpose workspaces as needed --- resolve/matrix/MatrixHandlerCuda.cpp | 43 +++++++++------------ resolve/matrix/MatrixHandlerHip.cpp | 28 ++++++-------- resolve/workspace/LinAlgWorkspaceCUDA.cpp | 11 +++++- resolve/workspace/LinAlgWorkspaceCUDA.hpp | 5 ++- resolve/workspace/LinAlgWorkspaceHIP.cpp | 11 +++++- resolve/workspace/LinAlgWorkspaceHIP.hpp | 1 + tests/unit/matrix/runMatrixHandlerTests.cpp | 5 ++- 7 files changed, 56 insertions(+), 48 deletions(-) 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/matrix/runMatrixHandlerTests.cpp b/tests/unit/matrix/runMatrixHandlerTests.cpp index 84e9a6011..cba5b6ad4 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 GPU workspace for a much larger transpose to exercise + // transpose-buffer growth rather than first-use-only 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); From bc8e955a628068756f081e4f4b7b0517cd497833 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 16:00:12 -0700 Subject: [PATCH 05/22] Remove redundant HyKKT device synchronizations --- resolve/hykkt/HyKKTSolver.cpp | 9 --------- 1 file changed, 9 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 40e1f83cc..6bba28658 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -258,11 +258,6 @@ 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 - } - r_y_copy_->copyFromExternal(r_y_, memspace_, memspace_); // check if this is redundant in later iterations @@ -557,10 +552,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_); From fad8ee42747c13bd1025b62cbf9350e788fdbd01 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 16:34:40 -0700 Subject: [PATCH 06/22] Reuse HyKKT conjugate gradient solver --- resolve/hykkt/HyKKTSolver.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 6bba28658..cbd495005 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -513,7 +513,15 @@ 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_); + if (!allocated_) + { + 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_); From edb239d033ed42b726a79bf156ffe7541d4e62f5 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 16:55:47 -0700 Subject: [PATCH 07/22] Reject incompatible HyKKT solver reuse --- resolve/hykkt/HyKKTSolver.cpp | 10 +++++++--- tests/unit/hykkt/HykktSolverTests.hpp | 9 +++++++++ 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index cbd495005..97a06d31c 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -79,9 +79,13 @@ namespace ReSolve 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; + // Arbitrary sparsity changes remain the caller's responsibility, but + // switching between empty and nonempty J_d invalidates cached HyKKT data. + if (!allocated_) + { + J_d_flag_ = J_d_flag; + } + status_ = (J_d_flag_ == J_d_flag); } /** diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index e197f96fa..65b79a372 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -195,6 +195,15 @@ namespace ReSolve real_type second_error = hykktSolver.solve(); status *= validateResult(second_error, tol); + // Changing J_d between nonempty and empty invalidates cached solver data. + matrix::Csr* J_d_empty = new matrix::Csr(J_d->getNumRows(), + J_d->getNumColumns(), + 0); + hykktSolver.setMatrixBlocks(H, D_s, J, J_d_empty); + real_type structure_error = hykktSolver.solve(); + status *= (structure_error == 1); + delete J_d_empty; + delete H; delete D_s; delete J; From fe27691afb6ca8ae9220b3912e012709ccda3e3c Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 18:04:57 -0700 Subject: [PATCH 08/22] Reuse HyKKT matrix allocation without J_d --- resolve/hykkt/HyKKTSolver.cpp | 5 ++++- tests/unit/hykkt/HykktSolverTests.hpp | 15 +++++++++++++++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 97a06d31c..6d9b646b4 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -313,7 +313,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_), diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index 65b79a372..f36e5cad9 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -202,6 +202,21 @@ namespace ReSolve hykktSolver.setMatrixBlocks(H, D_s, J, J_d_empty); real_type structure_error = hykktSolver.solve(); status *= (structure_error == 1); + + // Starting without J_d is valid and must remain reusable. + hykkt::HyKKTSolver no_jd_solver(n_x, m_d, m_c, memspace_); + no_jd_solver.setMatrixBlocks(H, D_s, 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 J_d_empty; delete H; From 723046483a74970830bf65766d8bbd010699d4f0 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 18:06:23 -0700 Subject: [PATCH 09/22] Refresh HyKKT SpGEMM gamma on solver reuse --- resolve/hykkt/HyKKTSolver.cpp | 1 + resolve/hykkt/spgemm/SpGEMM.hpp | 5 +++++ resolve/hykkt/spgemm/SpGEMMCpu.hpp | 5 +++++ resolve/hykkt/spgemm/SpGEMMCuda.hpp | 5 +++++ resolve/hykkt/spgemm/SpGEMMHip.hpp | 5 +++++ resolve/hykkt/spgemm/SpGEMMImpl.hpp | 2 ++ tests/unit/hykkt/HykktSolverTests.hpp | 2 ++ 7 files changed, 25 insertions(+) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 6d9b646b4..cda84bff1 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -407,6 +407,7 @@ namespace ReSolve { // Numerical values can change between solves while the sparsity pattern // remains fixed, so refresh the SpGEMM inputs before recomputing H_gamma. + spgemm_hgamma_->setAlpha(gamma_); spgemm_hgamma_->loadProductMatrices(J_tr_, J_); spgemm_hgamma_->loadSumMatrix(H_tilde_); spgemm_hgamma_->compute(); diff --git a/resolve/hykkt/spgemm/SpGEMM.hpp b/resolve/hykkt/spgemm/SpGEMM.hpp index edc3a267a..0bd14a4ff 100644 --- a/resolve/hykkt/spgemm/SpGEMM.hpp +++ b/resolve/hykkt/spgemm/SpGEMM.hpp @@ -22,6 +22,11 @@ namespace ReSolve SpGEMM(memory::MemorySpace memspace, real_type alpha, real_type beta); ~SpGEMM(); + void setAlpha(real_type alpha) + { + impl_->setAlpha(alpha); + } + 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.hpp b/resolve/hykkt/spgemm/SpGEMMCpu.hpp index 7c61c9ede..36c50b3a4 100644 --- a/resolve/hykkt/spgemm/SpGEMMCpu.hpp +++ b/resolve/hykkt/spgemm/SpGEMMCpu.hpp @@ -16,6 +16,11 @@ namespace ReSolve SpGEMMCpu(real_type alpha, real_type beta); ~SpGEMMCpu(); + void setAlpha(real_type alpha) override + { + alpha_ = alpha; + } + 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/SpGEMMCuda.hpp b/resolve/hykkt/spgemm/SpGEMMCuda.hpp index a39b2cdc9..d7f6cb79e 100644 --- a/resolve/hykkt/spgemm/SpGEMMCuda.hpp +++ b/resolve/hykkt/spgemm/SpGEMMCuda.hpp @@ -22,6 +22,11 @@ namespace ReSolve SpGEMMCuda(real_type alpha, real_type beta); ~SpGEMMCuda(); + void setAlpha(real_type alpha) override + { + alpha_ = alpha; + } + 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.hpp b/resolve/hykkt/spgemm/SpGEMMHip.hpp index 00627cd32..fafc187f3 100644 --- a/resolve/hykkt/spgemm/SpGEMMHip.hpp +++ b/resolve/hykkt/spgemm/SpGEMMHip.hpp @@ -22,6 +22,11 @@ namespace ReSolve SpGEMMHip(real_type alpha, real_type beta); ~SpGEMMHip(); + void setAlpha(real_type alpha) override + { + alpha_ = alpha; + } + 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/SpGEMMImpl.hpp b/resolve/hykkt/spgemm/SpGEMMImpl.hpp index ee8575dbf..6ddad76ce 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 setAlpha(real_type alpha) = 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/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index f36e5cad9..46bef3559 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -192,6 +192,8 @@ namespace ReSolve delete r_y_reuse; delete r_yd_reuse; + // 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); From eaa2c802468d8945173f4412df7cc90c8ae01862 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 18:10:55 -0700 Subject: [PATCH 10/22] Refresh HyKKT D_s data on solver reuse --- resolve/hykkt/HyKKTSolver.cpp | 4 +++- tests/unit/hykkt/HykktSolverTests.hpp | 17 ++++++++++++----- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index cda84bff1..13b395092 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -246,7 +246,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_); @@ -262,6 +261,9 @@ 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 } + + // 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 diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index 46bef3559..6a0e01ee8 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -142,13 +142,16 @@ namespace ReSolve testname += " N=" + std::to_string(N) + ", nnz =" + std::to_string(nnz) + '\n'; status *= validateResult(error, tol); - // Update D_s and restore data modified by the first solve. - real_type* D_s_values = D_s->getValues(memory::HOST); - for (index_type i = 0; i < D_s->getNnz(); ++i) + // Replace D_s and restore data modified by the first solve. + std::ifstream D_s_reuse_file(D_s_file_name); + matrix::Csr* D_s_reuse = io::createCsrFromFile(D_s_reuse_file, false); + + 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->setUpdated(memory::HOST); + D_s_reuse->setUpdated(memory::HOST); std::ifstream J_reuse_file(J_file_name); std::ifstream r_x_reuse_file(r_x_file_name); @@ -176,7 +179,8 @@ namespace ReSolve if (memspace_ == memory::DEVICE) { - restore_status |= D_s->syncData(memory::DEVICE); + 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); @@ -192,6 +196,8 @@ namespace ReSolve delete r_y_reuse; delete r_yd_reuse; + 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(); @@ -202,6 +208,7 @@ namespace ReSolve J_d->getNumColumns(), 0); hykktSolver.setMatrixBlocks(H, D_s, J, J_d_empty); + delete D_s_reuse; real_type structure_error = hykktSolver.solve(); status *= (structure_error == 1); From 10ea6147330b3e8e76be92bb7e2100db3d422596 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 18:38:37 -0700 Subject: [PATCH 11/22] Reuse HyKKT SCCG work vectors --- .../sccg/SchurComplementConjugateGradient.cpp | 29 ++++++++++--------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index 8fff342da..ee543c9d1 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -86,19 +86,22 @@ namespace ReSolve void SchurComplementConjugateGradient::setup() { - y_ = new vector::Vector(m_); - z_ = new vector::Vector(m_); - r_ = new vector::Vector(n_); - p_ = new vector::Vector(n_); - s_ = new vector::Vector(n_); - w_ = new vector::Vector(n_); - - y_->allocate(memspace_); - z_->allocate(memspace_); - r_->allocate(memspace_); - p_->allocate(memspace_); - s_->allocate(memspace_); - w_->allocate(memspace_); + if (!y_) + { + y_ = new vector::Vector(m_); + z_ = new vector::Vector(m_); + r_ = new vector::Vector(n_); + p_ = new vector::Vector(n_); + s_ = new vector::Vector(n_); + w_ = new vector::Vector(n_); + + y_->allocate(memspace_); + z_->allocate(memspace_); + r_->allocate(memspace_); + p_->allocate(memspace_); + s_->allocate(memspace_); + w_->allocate(memspace_); + } y_->setToZero(memspace_); z_->setToZero(memspace_); From fc7c3c78739f783fe69ef0a2fd137074aaec8965 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 27 Jul 2026 18:39:02 -0700 Subject: [PATCH 12/22] Handle zero residuals in HyKKT solves --- resolve/hykkt/HyKKTSolver.cpp | 6 +++++- .../hykkt/sccg/SchurComplementConjugateGradient.cpp | 6 ++++++ tests/unit/hykkt/HykktSCCGTests.hpp | 10 ++++++++++ tests/unit/hykkt/HykktSolverTests.hpp | 8 ++++++++ 4 files changed, 29 insertions(+), 1 deletion(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 13b395092..54ec2fa03 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -642,7 +642,11 @@ namespace ReSolve // Calculate final relative norm norm_resx_sq += norm_resy_sq; - real_type norm_res = sqrt(norm_resx_sq) / sqrt(norm_r_x_sq); + 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); allocated_ = true; diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index ee543c9d1..ae5459c82 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -124,6 +124,12 @@ namespace ReSolve 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/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index 4a3fe5331..e1e38e00c 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -104,6 +104,16 @@ 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); + sccg.setup(); + 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 6a0e01ee8..b22f1d1ca 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -203,6 +203,14 @@ namespace ReSolve real_type second_error = hykktSolver.solve(); status *= validateResult(second_error, tol); + // Verify an exact zero RHS is handled without producing 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); + // Changing J_d between nonempty and empty invalidates cached solver data. matrix::Csr* J_d_empty = new matrix::Csr(J_d->getNumRows(), J_d->getNumColumns(), From b9916c2fd7de6083e129f06be55516dd68b06cf8 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 31 Jul 2026 15:13:12 -0700 Subject: [PATCH 13/22] Clarify HyKKT reuse checks and residual reporting --- resolve/hykkt/HyKKTSolver.cpp | 19 +++++++++++-------- tests/unit/matrix/runMatrixHandlerTests.cpp | 4 ++-- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 54ec2fa03..58ed3db34 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -62,6 +62,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 @@ -160,13 +163,10 @@ 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"; + std::cout << "ERROR: Changing J_d between empty and nonempty is not " + "supported when reusing HyKKT.\n"; return 1; } @@ -266,7 +266,7 @@ namespace ReSolve 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_) { @@ -640,14 +640,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); 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); } - printf("||Ax-b||/||b|| = %32.32g\n\n", norm_res); allocated_ = true; diff --git a/tests/unit/matrix/runMatrixHandlerTests.cpp b/tests/unit/matrix/runMatrixHandlerTests.cpp index cba5b6ad4..03e0f72fa 100644 --- a/tests/unit/matrix/runMatrixHandlerTests.cpp +++ b/tests/unit/matrix/runMatrixHandlerTests.cpp @@ -46,8 +46,8 @@ 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 GPU workspace for a much larger transpose to exercise - // transpose-buffer growth rather than first-use-only allocation. + // 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); From b5f773e4cd77cf8ccf03fad99ff4edf5725ee298 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 31 Jul 2026 19:31:04 -0400 Subject: [PATCH 14/22] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) 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. From 7de350ba28af36a9443751de4b3ead681455f228 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Sat, 1 Aug 2026 23:00:33 -0700 Subject: [PATCH 15/22] Clean up HyKKT SpGEMM setup --- resolve/hykkt/HyKKTSolver.cpp | 10 +++++----- resolve/hykkt/spgemm/SpGEMM.cpp | 10 ++++++++++ resolve/hykkt/spgemm/SpGEMM.hpp | 5 +---- resolve/hykkt/spgemm/SpGEMMCpu.cpp | 5 +++++ resolve/hykkt/spgemm/SpGEMMCpu.hpp | 5 +---- resolve/hykkt/spgemm/SpGEMMCuda.cpp | 5 +++++ resolve/hykkt/spgemm/SpGEMMCuda.hpp | 5 +---- resolve/hykkt/spgemm/SpGEMMHip.cpp | 5 +++++ resolve/hykkt/spgemm/SpGEMMHip.hpp | 5 +---- tests/unit/hykkt/HykktSolverTests.hpp | 5 +++-- 10 files changed, 37 insertions(+), 23 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 58ed3db34..48e447f89 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -214,8 +214,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 @@ -392,9 +392,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() } /* @@ -412,6 +409,9 @@ namespace ReSolve spgemm_hgamma_->setAlpha(gamma_); 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_); diff --git a/resolve/hykkt/spgemm/SpGEMM.cpp b/resolve/hykkt/spgemm/SpGEMM.cpp index 1170f919f..8167f4432 100644 --- a/resolve/hykkt/spgemm/SpGEMM.cpp +++ b/resolve/hykkt/spgemm/SpGEMM.cpp @@ -53,6 +53,16 @@ namespace ReSolve delete impl_; } + /** + * Updates the scalar multiplier for the matrix product. + * + * @param[in] alpha - Scalar multiplier for the product. + */ + void SpGEMM::setAlpha(real_type alpha) + { + impl_->setAlpha(alpha); + } + /** * 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 0bd14a4ff..c90292320 100644 --- a/resolve/hykkt/spgemm/SpGEMM.hpp +++ b/resolve/hykkt/spgemm/SpGEMM.hpp @@ -22,10 +22,7 @@ namespace ReSolve SpGEMM(memory::MemorySpace memspace, real_type alpha, real_type beta); ~SpGEMM(); - void setAlpha(real_type alpha) - { - impl_->setAlpha(alpha); - } + void setAlpha(real_type alpha); void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); void loadSumMatrix(matrix::Csr* D); diff --git a/resolve/hykkt/spgemm/SpGEMMCpu.cpp b/resolve/hykkt/spgemm/SpGEMMCpu.cpp index f634231d8..3391dffdb 100644 --- a/resolve/hykkt/spgemm/SpGEMMCpu.cpp +++ b/resolve/hykkt/spgemm/SpGEMMCpu.cpp @@ -47,6 +47,11 @@ namespace ReSolve cholmod_finish(&Common_); } + void SpGEMMCpu::setAlpha(real_type alpha) + { + alpha_ = alpha; + } + 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 36c50b3a4..d2c8e44a7 100644 --- a/resolve/hykkt/spgemm/SpGEMMCpu.hpp +++ b/resolve/hykkt/spgemm/SpGEMMCpu.hpp @@ -16,10 +16,7 @@ namespace ReSolve SpGEMMCpu(real_type alpha, real_type beta); ~SpGEMMCpu(); - void setAlpha(real_type alpha) override - { - alpha_ = alpha; - } + void setAlpha(real_type alpha) override; void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); void loadSumMatrix(matrix::Csr* D); diff --git a/resolve/hykkt/spgemm/SpGEMMCuda.cpp b/resolve/hykkt/spgemm/SpGEMMCuda.cpp index a9f93b60e..cf35fb075 100644 --- a/resolve/hykkt/spgemm/SpGEMMCuda.cpp +++ b/resolve/hykkt/spgemm/SpGEMMCuda.cpp @@ -29,6 +29,11 @@ namespace ReSolve cusparseDestroy(handle_); } + void SpGEMMCuda::setAlpha(real_type alpha) + { + alpha_ = alpha; + } + void SpGEMMCuda::loadProductMatrices(matrix::Csr* A, matrix::Csr* B) { A_descr_ = convertToCusparseType(A); diff --git a/resolve/hykkt/spgemm/SpGEMMCuda.hpp b/resolve/hykkt/spgemm/SpGEMMCuda.hpp index d7f6cb79e..0085e3600 100644 --- a/resolve/hykkt/spgemm/SpGEMMCuda.hpp +++ b/resolve/hykkt/spgemm/SpGEMMCuda.hpp @@ -22,10 +22,7 @@ namespace ReSolve SpGEMMCuda(real_type alpha, real_type beta); ~SpGEMMCuda(); - void setAlpha(real_type alpha) override - { - alpha_ = alpha; - } + void setAlpha(real_type alpha) override; void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); void loadSumMatrix(matrix::Csr* D); diff --git a/resolve/hykkt/spgemm/SpGEMMHip.cpp b/resolve/hykkt/spgemm/SpGEMMHip.cpp index 38fafc928..1be5f76b4 100644 --- a/resolve/hykkt/spgemm/SpGEMMHip.cpp +++ b/resolve/hykkt/spgemm/SpGEMMHip.cpp @@ -30,6 +30,11 @@ namespace ReSolve mem_.deleteOnDevice(buffer_); } + void SpGEMMHip::setAlpha(real_type alpha) + { + alpha_ = alpha; + } + 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 fafc187f3..72311fc12 100644 --- a/resolve/hykkt/spgemm/SpGEMMHip.hpp +++ b/resolve/hykkt/spgemm/SpGEMMHip.hpp @@ -22,10 +22,7 @@ namespace ReSolve SpGEMMHip(real_type alpha, real_type beta); ~SpGEMMHip(); - void setAlpha(real_type alpha) override - { - alpha_ = alpha; - } + void setAlpha(real_type alpha) override; void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); void loadSumMatrix(matrix::Csr* D); diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index b22f1d1ca..e9ca1f747 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -203,7 +203,7 @@ namespace ReSolve real_type second_error = hykktSolver.solve(); status *= validateResult(second_error, tol); - // Verify an exact zero RHS is handled without producing NaNs. + // Check that a zero RHS doesn't result in NaNs. r_x->setToZero(memspace_); r_s->setToZero(memspace_); r_y->setToZero(memspace_); @@ -211,7 +211,8 @@ namespace ReSolve real_type zero_rhs_error = hykktSolver.solve(); status *= validateResult(zero_rhs_error, tol); - // Changing J_d between nonempty and empty invalidates cached solver data. + // 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); From fe5fd93bc97f09ba105266d2e0e4841eaa486dfd Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Sun, 2 Aug 2026 14:26:10 -0700 Subject: [PATCH 16/22] Move HyKKT SCCG allocation into solve --- resolve/hykkt/HyKKTSolver.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 48e447f89..6454a1ca8 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -204,6 +204,15 @@ namespace ReSolve } computeHgammaFactorization(); + if (!allocated_) + { + sccg_ = new SchurComplementConjugateGradient(J_->getNumRows(), + J_->getNumColumns(), + cholesky_, + matrixHandler_, + vectorHandler_, + memspace_); + } setupConjugateGradient(); computeConjugateGradient(); @@ -523,15 +532,6 @@ namespace ReSolve schur_->copyFromExternal(r_y_, memspace_, memspace_); matrixHandler_->matvec(J_perm_, omega_perm_, schur_, &ONE, &MINUS_ONE, memspace_); - if (!allocated_) - { - 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_); From b5a0b049386210a5e3fca31138837f33916b6d41 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Sun, 2 Aug 2026 14:52:26 -0700 Subject: [PATCH 17/22] Reset HyKKT SCCG state in solve --- resolve/hykkt/HyKKTSolver.cpp | 2 +- .../sccg/SchurComplementConjugateGradient.cpp | 39 +++++++++---------- .../sccg/SchurComplementConjugateGradient.hpp | 5 +++ tests/unit/hykkt/HykktSCCGTests.hpp | 1 - 4 files changed, 24 insertions(+), 23 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 6454a1ca8..a454d037d 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -212,6 +212,7 @@ namespace ReSolve matrixHandler_, vectorHandler_, memspace_); + sccg_->setup(); } setupConjugateGradient(); computeConjugateGradient(); @@ -535,7 +536,6 @@ namespace ReSolve sccg_->addMatrixInfo(J_perm_, J_tr_perm_); y_->setToZero(memspace_); sccg_->addVectorInfo(y_, schur_); - sccg_->setup(); } /** diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index ae5459c82..791f8cd9f 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -86,22 +86,24 @@ namespace ReSolve void SchurComplementConjugateGradient::setup() { - if (!y_) - { - y_ = new vector::Vector(m_); - z_ = new vector::Vector(m_); - r_ = new vector::Vector(n_); - p_ = new vector::Vector(n_); - s_ = new vector::Vector(n_); - w_ = new vector::Vector(n_); - - y_->allocate(memspace_); - z_->allocate(memspace_); - r_->allocate(memspace_); - p_->allocate(memspace_); - s_->allocate(memspace_); - w_->allocate(memspace_); - } + y_ = new vector::Vector(m_); + z_ = new vector::Vector(m_); + r_ = new vector::Vector(n_); + p_ = new vector::Vector(n_); + s_ = new vector::Vector(n_); + w_ = new vector::Vector(n_); + + y_->allocate(memspace_); + z_->allocate(memspace_); + r_->allocate(memspace_); + p_->allocate(memspace_); + s_->allocate(memspace_); + w_->allocate(memspace_); + } + + int SchurComplementConjugateGradient::solve() + { + using namespace constants; y_->setToZero(memspace_); z_->setToZero(memspace_); @@ -113,11 +115,6 @@ 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_); 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/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index e1e38e00c..f3958b310 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -109,7 +109,6 @@ namespace ReSolve x_0->setToZero(memspace_); b->setToZero(memspace_); sccg.addVectorInfo(x_0, b); - sccg.setup(); int zero_residual_converged_n = sccg.solve(); status *= (zero_residual_converged_n == 0); status *= (vector_handler_.dot(x_0, x_0, memspace_) <= sccg_tol); From b1e5f22b7063736c3dd14b7663a32614640f461a Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Sun, 2 Aug 2026 15:11:31 -0700 Subject: [PATCH 18/22] Set HyKKT SpGEMM coefficients --- resolve/hykkt/HyKKTSolver.cpp | 2 +- resolve/hykkt/spgemm/SpGEMM.cpp | 9 +++-- resolve/hykkt/spgemm/SpGEMM.hpp | 2 +- resolve/hykkt/spgemm/SpGEMMCpu.cpp | 3 +- resolve/hykkt/spgemm/SpGEMMCpu.hpp | 2 +- resolve/hykkt/spgemm/SpGEMMCuda.cpp | 3 +- resolve/hykkt/spgemm/SpGEMMCuda.hpp | 2 +- resolve/hykkt/spgemm/SpGEMMHip.cpp | 3 +- resolve/hykkt/spgemm/SpGEMMHip.hpp | 2 +- resolve/hykkt/spgemm/SpGEMMImpl.hpp | 2 +- tests/unit/hykkt/HykktSpGEMMTests.hpp | 55 ++++++++++++++++++++++++++- 11 files changed, 71 insertions(+), 14 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index a454d037d..bb569e031 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -416,7 +416,7 @@ namespace ReSolve { // Numerical values can change between solves while the sparsity pattern // remains fixed, so refresh the SpGEMM inputs before recomputing H_gamma. - spgemm_hgamma_->setAlpha(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 diff --git a/resolve/hykkt/spgemm/SpGEMM.cpp b/resolve/hykkt/spgemm/SpGEMM.cpp index 8167f4432..2610dcdd1 100644 --- a/resolve/hykkt/spgemm/SpGEMM.cpp +++ b/resolve/hykkt/spgemm/SpGEMM.cpp @@ -54,13 +54,14 @@ namespace ReSolve } /** - * Updates the scalar multiplier for the matrix product. + * Updates the coefficients for the SpGEMM operation. * - * @param[in] alpha - Scalar multiplier for the product. + * @param[in] alpha - Scalar multiplier for the matrix product. + * @param[in] beta - Scalar multiplier for the sum matrix. */ - void SpGEMM::setAlpha(real_type alpha) + void SpGEMM::setCoefficients(real_type alpha, real_type beta) { - impl_->setAlpha(alpha); + impl_->setCoefficients(alpha, beta); } /** diff --git a/resolve/hykkt/spgemm/SpGEMM.hpp b/resolve/hykkt/spgemm/SpGEMM.hpp index c90292320..21a3390e5 100644 --- a/resolve/hykkt/spgemm/SpGEMM.hpp +++ b/resolve/hykkt/spgemm/SpGEMM.hpp @@ -22,7 +22,7 @@ namespace ReSolve SpGEMM(memory::MemorySpace memspace, real_type alpha, real_type beta); ~SpGEMM(); - void setAlpha(real_type alpha); + void setCoefficients(real_type alpha, real_type beta); void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); void loadSumMatrix(matrix::Csr* D); diff --git a/resolve/hykkt/spgemm/SpGEMMCpu.cpp b/resolve/hykkt/spgemm/SpGEMMCpu.cpp index 3391dffdb..93ad651b3 100644 --- a/resolve/hykkt/spgemm/SpGEMMCpu.cpp +++ b/resolve/hykkt/spgemm/SpGEMMCpu.cpp @@ -47,9 +47,10 @@ namespace ReSolve cholmod_finish(&Common_); } - void SpGEMMCpu::setAlpha(real_type alpha) + void SpGEMMCpu::setCoefficients(real_type alpha, real_type beta) { alpha_ = alpha; + beta_ = beta; } void SpGEMMCpu::loadProductMatrices(matrix::Csr* A, matrix::Csr* B) diff --git a/resolve/hykkt/spgemm/SpGEMMCpu.hpp b/resolve/hykkt/spgemm/SpGEMMCpu.hpp index d2c8e44a7..ccf018bc8 100644 --- a/resolve/hykkt/spgemm/SpGEMMCpu.hpp +++ b/resolve/hykkt/spgemm/SpGEMMCpu.hpp @@ -16,7 +16,7 @@ namespace ReSolve SpGEMMCpu(real_type alpha, real_type beta); ~SpGEMMCpu(); - void setAlpha(real_type alpha) override; + void setCoefficients(real_type alpha, real_type beta) override; void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); void loadSumMatrix(matrix::Csr* D); diff --git a/resolve/hykkt/spgemm/SpGEMMCuda.cpp b/resolve/hykkt/spgemm/SpGEMMCuda.cpp index cf35fb075..0d4310a9b 100644 --- a/resolve/hykkt/spgemm/SpGEMMCuda.cpp +++ b/resolve/hykkt/spgemm/SpGEMMCuda.cpp @@ -29,9 +29,10 @@ namespace ReSolve cusparseDestroy(handle_); } - void SpGEMMCuda::setAlpha(real_type alpha) + void SpGEMMCuda::setCoefficients(real_type alpha, real_type beta) { alpha_ = alpha; + beta_ = beta; } void SpGEMMCuda::loadProductMatrices(matrix::Csr* A, matrix::Csr* B) diff --git a/resolve/hykkt/spgemm/SpGEMMCuda.hpp b/resolve/hykkt/spgemm/SpGEMMCuda.hpp index 0085e3600..8bd22fc99 100644 --- a/resolve/hykkt/spgemm/SpGEMMCuda.hpp +++ b/resolve/hykkt/spgemm/SpGEMMCuda.hpp @@ -22,7 +22,7 @@ namespace ReSolve SpGEMMCuda(real_type alpha, real_type beta); ~SpGEMMCuda(); - void setAlpha(real_type alpha) override; + void setCoefficients(real_type alpha, real_type beta) override; void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); void loadSumMatrix(matrix::Csr* D); diff --git a/resolve/hykkt/spgemm/SpGEMMHip.cpp b/resolve/hykkt/spgemm/SpGEMMHip.cpp index 1be5f76b4..fc2c62bc2 100644 --- a/resolve/hykkt/spgemm/SpGEMMHip.cpp +++ b/resolve/hykkt/spgemm/SpGEMMHip.cpp @@ -30,9 +30,10 @@ namespace ReSolve mem_.deleteOnDevice(buffer_); } - void SpGEMMHip::setAlpha(real_type alpha) + void SpGEMMHip::setCoefficients(real_type alpha, real_type beta) { alpha_ = alpha; + beta_ = beta; } void SpGEMMHip::loadProductMatrices(matrix::Csr* A, matrix::Csr* B) diff --git a/resolve/hykkt/spgemm/SpGEMMHip.hpp b/resolve/hykkt/spgemm/SpGEMMHip.hpp index 72311fc12..d4b430289 100644 --- a/resolve/hykkt/spgemm/SpGEMMHip.hpp +++ b/resolve/hykkt/spgemm/SpGEMMHip.hpp @@ -22,7 +22,7 @@ namespace ReSolve SpGEMMHip(real_type alpha, real_type beta); ~SpGEMMHip(); - void setAlpha(real_type alpha) override; + void setCoefficients(real_type alpha, real_type beta) override; void loadProductMatrices(matrix::Csr* A, matrix::Csr* B); void loadSumMatrix(matrix::Csr* D); diff --git a/resolve/hykkt/spgemm/SpGEMMImpl.hpp b/resolve/hykkt/spgemm/SpGEMMImpl.hpp index 6ddad76ce..e70f0a0bd 100644 --- a/resolve/hykkt/spgemm/SpGEMMImpl.hpp +++ b/resolve/hykkt/spgemm/SpGEMMImpl.hpp @@ -23,7 +23,7 @@ namespace ReSolve SpGEMMImpl() = default; virtual ~SpGEMMImpl() = default; - virtual void setAlpha(real_type alpha) = 0; + 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; diff --git a/tests/unit/hykkt/HykktSpGEMMTests.hpp b/tests/unit/hykkt/HykktSpGEMMTests.hpp index ecaeb987a..b0ecc6284 100644 --- a/tests/unit/hykkt/HykktSpGEMMTests.hpp +++ b/tests/unit/hykkt/HykktSpGEMMTests.hpp @@ -292,6 +292,59 @@ 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 +371,7 @@ namespace ReSolve E->syncData(memory::HOST); } - status *= verifyResult(E, 4.0); + status *= verifyResult(E, 6.0); delete A; delete B; From e2806038a0a9bc89a323dcf558b8ea033c559d69 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Sun, 2 Aug 2026 15:20:42 -0700 Subject: [PATCH 19/22] Update HyKKT reuse inputs in place --- tests/unit/hykkt/HykktSolverTests.hpp | 61 ++++++++++----------------- 1 file changed, 23 insertions(+), 38 deletions(-) diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index e9ca1f747..173c50278 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -142,43 +142,35 @@ namespace ReSolve testname += " N=" + std::to_string(N) + ", nnz =" + std::to_string(nnz) + '\n'; status *= validateResult(error, tol); - // Replace D_s and restore data modified by the first solve. + // Replace D_s to exercise its pointer refresh; restore other inputs + // in place. std::ifstream D_s_reuse_file(D_s_file_name); - matrix::Csr* D_s_reuse = io::createCsrFromFile(D_s_reuse_file, false); + 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); - - std::ifstream J_reuse_file(J_file_name); - std::ifstream r_x_reuse_file(r_x_file_name); - std::ifstream r_s_reuse_file(r_s_file_name); - std::ifstream r_y_reuse_file(r_y_file_name); - std::ifstream r_yd_reuse_file(r_yd_file_name); - - matrix::Csr* J_reuse = io::createCsrFromFile(J_reuse_file, false); - - vector::Vector* r_x_reuse = io::createVectorFromFile(r_x_reuse_file); - vector::Vector* r_s_reuse = io::createVectorFromFile(r_s_reuse_file); - vector::Vector* r_y_reuse = io::createVectorFromFile(r_y_reuse_file); - vector::Vector* r_yd_reuse = io::createVectorFromFile(r_yd_reuse_file); - int restore_status = 0; - restore_status |= J->copyFromExternal(J_reuse->getRowData(memory::HOST), - J_reuse->getColData(memory::HOST), - J_reuse->getValues(memory::HOST), - memory::HOST, - memory::HOST); - restore_status |= r_x->copyFromExternal(r_x_reuse, memory::HOST, memory::HOST); - restore_status |= r_s->copyFromExternal(r_s_reuse, memory::HOST, memory::HOST); - restore_status |= r_y->copyFromExternal(r_y_reuse, memory::HOST, memory::HOST); - restore_status |= r_yd->copyFromExternal(r_yd_reuse, memory::HOST, memory::HOST); + 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); @@ -186,16 +178,9 @@ namespace ReSolve 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); } - status *= (restore_status == 0); - - delete J_reuse; - delete r_x_reuse; - delete r_s_reuse; - delete r_y_reuse; - delete r_yd_reuse; - hykktSolver.setMatrixBlocks(H, D_s_reuse, J, J_d); // Change gamma to verify the cached SpGEMM coefficient is refreshed. @@ -216,14 +201,13 @@ namespace ReSolve matrix::Csr* J_d_empty = new matrix::Csr(J_d->getNumRows(), J_d->getNumColumns(), 0); - hykktSolver.setMatrixBlocks(H, D_s, J, J_d_empty); - delete D_s_reuse; + hykktSolver.setMatrixBlocks(H, D_s_reuse, J, J_d_empty); real_type structure_error = hykktSolver.solve(); status *= (structure_error == 1); - // Starting without J_d is valid and must remain reusable. + // 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, J, J_d_empty); + 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); @@ -235,6 +219,7 @@ namespace ReSolve 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; From 1ff0ea77bbb10c1be8252dae8ee09829cef6dae8 Mon Sep 17 00:00:00 2001 From: tamar-dewilde Date: Mon, 3 Aug 2026 01:33:51 +0000 Subject: [PATCH 20/22] Apply pre-commmit fixes --- tests/unit/hykkt/HykktSpGEMMTests.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/unit/hykkt/HykktSpGEMMTests.hpp b/tests/unit/hykkt/HykktSpGEMMTests.hpp index b0ecc6284..24fd945c2 100644 --- a/tests/unit/hykkt/HykktSpGEMMTests.hpp +++ b/tests/unit/hykkt/HykktSpGEMMTests.hpp @@ -321,8 +321,7 @@ namespace ReSolve 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) + if (fabs(E->getValues(memory::HOST)[i] - equal_coefficient_values[i]) > 1e-12) { beta_changed_result = true; break; From b782d0b6f4ba7e0362bbf1b87bb62f29f093cb65 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 3 Aug 2026 09:26:55 -0700 Subject: [PATCH 21/22] Report HyKKT block errors from setMatrixBlocks --- resolve/hykkt/HyKKTSolver.cpp | 29 ++++++++++++++++----------- resolve/hykkt/HyKKTSolver.hpp | 6 +----- tests/unit/hykkt/HykktSolverTests.hpp | 6 +++--- 3 files changed, 21 insertions(+), 20 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index bb569e031..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 @@ -73,22 +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; - // Arbitrary sparsity changes remain the caller's responsibility, but - // switching between empty and nonempty J_d invalidates cached HyKKT data. if (!allocated_) { J_d_flag_ = J_d_flag; } - status_ = (J_d_flag_ == J_d_flag); + return 0; } /** @@ -163,13 +175,6 @@ namespace ReSolve */ real_type hykkt::HyKKTSolver::solve() { - if (!status_ && allocated_) - { - std::cout << "ERROR: Changing J_d between empty and nonempty is not " - "supported when reusing HyKKT.\n"; - return 1; - } - setupParameters(); if (!allocated_) 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/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index 173c50278..b71232a91 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -201,9 +201,9 @@ namespace ReSolve matrix::Csr* J_d_empty = new matrix::Csr(J_d->getNumRows(), J_d->getNumColumns(), 0); - hykktSolver.setMatrixBlocks(H, D_s_reuse, J, J_d_empty); - real_type structure_error = hykktSolver.solve(); - status *= (structure_error == 1); + + 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_); From 70234c0f091bb1f6318aba74065f6c6ebb2e2558 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 3 Aug 2026 10:13:40 -0700 Subject: [PATCH 22/22] Fix SpGEMM override warnings --- resolve/hykkt/spgemm/SpGEMMCpu.hpp | 8 ++++---- resolve/hykkt/spgemm/SpGEMMHip.hpp | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/resolve/hykkt/spgemm/SpGEMMCpu.hpp b/resolve/hykkt/spgemm/SpGEMMCpu.hpp index ccf018bc8..3522be0ca 100644 --- a/resolve/hykkt/spgemm/SpGEMMCpu.hpp +++ b/resolve/hykkt/spgemm/SpGEMMCpu.hpp @@ -18,11 +18,11 @@ namespace ReSolve 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); + void loadProductMatrices(matrix::Csr* A, matrix::Csr* B) override; + void loadSumMatrix(matrix::Csr* D) override; + void loadResultMatrix(matrix::Csr** E_ptr) override; - void compute(); + void compute() override; private: real_type alpha_; diff --git a/resolve/hykkt/spgemm/SpGEMMHip.hpp b/resolve/hykkt/spgemm/SpGEMMHip.hpp index d4b430289..51423db9d 100644 --- a/resolve/hykkt/spgemm/SpGEMMHip.hpp +++ b/resolve/hykkt/spgemm/SpGEMMHip.hpp @@ -24,11 +24,11 @@ namespace ReSolve 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); + void loadProductMatrices(matrix::Csr* A, matrix::Csr* B) override; + void loadSumMatrix(matrix::Csr* D) override; + void loadResultMatrix(matrix::Csr** E_ptr) override; - void compute(); + void compute() override; private: MemoryHandler mem_;