From a5eac9ff4ceb2cc567b390469bf3d2844628784f Mon Sep 17 00:00:00 2001 From: Ishika Roy Date: Thu, 2 Jul 2026 05:26:01 +0000 Subject: [PATCH] add solver persistence pipeline --- .../pdlp/solver_settings.hpp | 8 + .../utilities/cython_solve.hpp | 4 +- .../utilities/cython_types.hpp | 4 + .../utilities/lp_solve_session.hpp | 65 ++ .../utilities/solver_cache_profiler.hpp | 152 +++ cpp/src/barrier/CMakeLists.txt | 1 + cpp/src/barrier/barrier.cu | 716 ++++++++++++-- cpp/src/barrier/barrier.hpp | 9 +- .../barrier_factorization_sparsity_hash.cu | 25 + .../barrier_factorization_sparsity_hash.hpp | 126 +++ cpp/src/barrier/barrier_symbolic_cache.hpp | 82 ++ cpp/src/barrier/cusparse_view.cu | 10 + cpp/src/barrier/cusparse_view.hpp | 2 + cpp/src/barrier/device_sparse_matrix.cuh | 8 + cpp/src/barrier/sparse_cholesky.cuh | 172 +++- cpp/src/dual_simplex/solve.cpp | 16 +- cpp/src/dual_simplex/solve.hpp | 22 +- cpp/src/pdlp/CMakeLists.txt | 1 + cpp/src/pdlp/solve.cu | 76 +- cpp/src/pdlp/utilities/cython_solve.cu | 80 +- cpp/src/pdlp/utilities/lp_solve_session.cu | 83 ++ .../cuopt/cuopt/linear_programming/problem.py | 160 ++- .../linear_programming/solution/solution.py | 7 + .../linear_programming/solver/solver.pxd | 8 + .../cuopt/linear_programming/solver/solver.py | 6 +- .../solver/solver_wrapper.pyx | 53 +- .../solver_settings/solver_settings.pxd | 11 + .../solver_settings/solver_settings.pyx | 11 + .../session_cache_helpers.py | 350 +++++++ .../test_lp_solve_session.py | 221 +++++ script_perf_eval.py | 927 ++++++++++++++++++ script_session_cache_tests.py | 150 +++ 32 files changed, 3362 insertions(+), 204 deletions(-) create mode 100644 cpp/include/cuopt/linear_programming/utilities/lp_solve_session.hpp create mode 100644 cpp/include/cuopt/linear_programming/utilities/solver_cache_profiler.hpp create mode 100644 cpp/src/barrier/barrier_factorization_sparsity_hash.cu create mode 100644 cpp/src/barrier/barrier_factorization_sparsity_hash.hpp create mode 100644 cpp/src/barrier/barrier_symbolic_cache.hpp create mode 100644 cpp/src/pdlp/utilities/lp_solve_session.cu create mode 100644 python/cuopt/cuopt/tests/linear_programming/session_cache_helpers.py create mode 100644 python/cuopt/cuopt/tests/linear_programming/test_lp_solve_session.py create mode 100644 script_perf_eval.py create mode 100644 script_session_cache_tests.py diff --git a/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp b/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp index b30286f9ce..500984f6f9 100644 --- a/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp @@ -21,6 +21,10 @@ #include +namespace cuopt::cython { +class lp_solve_session_t; +} + namespace cuopt::linear_programming { // Forward declare solver_settings_t for friend class @@ -326,6 +330,10 @@ class pdlp_solver_settings_t { // Used to force batch PDLP to solve a subbatch of the problems at a time // The 0 default value will make the solver use its heuristic to determine the subbatch size i_t fixed_batch_size{0}; + /** When true, first GPU barrier/QCQP solve returns an ``lp_solve_session_t`` capsule. */ + bool session_enabled{false}; + /** Non-owning session pointer set by ``call_solve`` for barrier symbolic reuse. */ + cuopt::cython::lp_solve_session_t* lp_solve_session{nullptr}; private: /** Initial primal solution */ diff --git a/cpp/include/cuopt/linear_programming/utilities/cython_solve.hpp b/cpp/include/cuopt/linear_programming/utilities/cython_solve.hpp index 9f4fc93923..a3092f81b2 100644 --- a/cpp/include/cuopt/linear_programming/utilities/cython_solve.hpp +++ b/cpp/include/cuopt/linear_programming/utilities/cython_solve.hpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -55,7 +56,8 @@ std::unique_ptr call_solve( cuopt::linear_programming::io::data_model_view_t*, linear_programming::solver_settings_t*, unsigned int flags = cudaStreamNonBlocking, - bool is_batch_mode = false); + bool is_batch_mode = false, + lp_solve_session_t* session_in = nullptr); std::pair>, double> solve_batch_remote( std::vector*>, diff --git a/cpp/include/cuopt/linear_programming/utilities/cython_types.hpp b/cpp/include/cuopt/linear_programming/utilities/cython_types.hpp index 20db133512..8385b688f6 100644 --- a/cpp/include/cuopt/linear_programming/utilities/cython_types.hpp +++ b/cpp/include/cuopt/linear_programming/utilities/cython_types.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include @@ -85,6 +86,9 @@ struct linear_programming_ret_t { double solve_time_{}; linear_programming::method_t solved_by_{}; + /** GPU barrier session (stream + handle + symbolic cache); moved to Python capsule when set. */ + std::unique_ptr lp_solve_session; + bool is_gpu() const { return std::holds_alternative(solutions_); } }; diff --git a/cpp/include/cuopt/linear_programming/utilities/lp_solve_session.hpp b/cpp/include/cuopt/linear_programming/utilities/lp_solve_session.hpp new file mode 100644 index 0000000000..98b7b3e0a9 --- /dev/null +++ b/cpp/include/cuopt/linear_programming/utilities/lp_solve_session.hpp @@ -0,0 +1,65 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include + +#include +#include + +namespace cuopt::linear_programming::dual_simplex { +template +class iteration_data_t; +template +struct barrier_symbolic_cache_t; + +template +void barrier_store_symbolic_cache_from_iteration_data(iteration_data_t& data, + barrier_symbolic_cache_t& cache); +} // namespace cuopt::linear_programming::dual_simplex + +namespace cuopt::cython { + +/** + * @brief Lean GPU solve session: owns RAFT handle + stream and optional barrier symbolic cache. + * + * Created on first solve when session_enabled; reused on subsequent solves with the same capsule. + * Per-solve state (optimization_problem_t, presolve, barrier_lp) remains stack-local. + */ +class lp_solve_session_t { + public: + static std::unique_ptr create(unsigned stream_flags); + + lp_solve_session_t(lp_solve_session_t&&) noexcept; + lp_solve_session_t& operator=(lp_solve_session_t&&) noexcept; + ~lp_solve_session_t(); + + [[nodiscard]] raft::handle_t* handle_ptr(); + [[nodiscard]] raft::handle_t const* handle_ptr() const; + [[nodiscard]] rmm::cuda_stream_view stream_view() const; + + /** + * @brief Returns cached symbolic state when valid and @p handle matches the stored handle. + */ + [[nodiscard]] linear_programming::dual_simplex::barrier_symbolic_cache_t* + symbolic_cache_for_reuse(raft::handle_t const* handle); + + void clear_symbolic_cache(); + + void store_symbolic_cache( + linear_programming::dual_simplex::iteration_data_t& data); + + private: + lp_solve_session_t(std::unique_ptr stream, + std::unique_ptr handle); + + struct impl; + std::unique_ptr impl_; +}; + +} // namespace cuopt::cython diff --git a/cpp/include/cuopt/linear_programming/utilities/solver_cache_profiler.hpp b/cpp/include/cuopt/linear_programming/utilities/solver_cache_profiler.hpp new file mode 100644 index 0000000000..0b72d45666 --- /dev/null +++ b/cpp/include/cuopt/linear_programming/utilities/solver_cache_profiler.hpp @@ -0,0 +1,152 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace cuopt::linear_programming::cache_profile { + +enum class cache_id : int { + C01 = 0, // raft::handle_t + stream + C02, // cuBLAS / cuSparse warmup (init_handler) + C03, // problem fingerprint (structural hash) + C04, // augmented vs ADAT choice + C05, // ADAT / augmented sparsity pattern + C06, // cuDSS handle + config + C07, // cuDSS symbolic factorization + C08, // dense-column / SOC layout metadata + C09, // device buffer allocation (iteration_data setup) + COUNT +}; + +inline constexpr int num_cache_ids = static_cast(cache_id::COUNT); + +inline const char* cache_id_label(cache_id id) +{ + switch (id) { + case cache_id::C01: return "C01 raft handle+stream"; + case cache_id::C02: return "C02 cuBLAS/cuSparse warmup"; + case cache_id::C03: return "C03 problem fingerprint"; + case cache_id::C04: return "C04 augmented vs ADAT choice"; + case cache_id::C05: return "C05 KKT sparsity pattern"; + case cache_id::C06: return "C06 cuDSS handle+config"; + case cache_id::C07: return "C07 cuDSS symbolic factorization"; + case cache_id::C08: return "C08 dense-column/SOC layout"; + case cache_id::C09: return "C09 device buffer allocation"; + default: return "C?? unknown"; + } +} + +class profiler_t { + public: + static profiler_t& instance() + { + static profiler_t prof; + return prof; + } + + bool enabled() const { return enabled_; } + + void reset() + { + times_.fill(0.0); + } + + void add(cache_id id, double seconds) + { + if (!enabled_) { return; } + times_[static_cast(id)] += seconds; + } + + double get(cache_id id) const { return times_[static_cast(id)]; } + + double total_measured() const + { + double sum = 0.0; + for (double t : times_) { + sum += t; + } + return sum; + } + + void log_summary() const + { + if (!enabled_) { return; } + auto emit = [](const char* fmt, ...) { + va_list args; + va_start(args, fmt); + char buf[512]; + vsnprintf(buf, sizeof(buf), fmt, args); + va_end(args); + CUOPT_LOG_INFO("%s", buf); + fprintf(stderr, "%s\n", buf); + }; + emit("=== Solver cache profile (ms) ==="); + for (int i = 0; i < num_cache_ids; ++i) { + const double ms = times_[i] * 1000.0; + emit("Cache profile: %s %.3f", cache_id_label(static_cast(i)), ms); + } + emit("Cache profile: TOTAL measured %.3f", total_measured() * 1000.0); + emit("=== End solver cache profile ==="); + } + + private: + profiler_t() + { + const char* env = std::getenv("CUOPT_CACHE_PROFILE"); + enabled_ = env != nullptr && env[0] != '\0' && std::strcmp(env, "0") != 0; + } + + bool enabled_{false}; + std::array times_{}; +}; + +inline bool enabled() { return profiler_t::instance().enabled(); } + +inline void reset() { profiler_t::instance().reset(); } + +inline void add(cache_id id, double seconds) { profiler_t::instance().add(id, seconds); } + +inline void log_summary() { profiler_t::instance().log_summary(); } + +class scoped_timer_t { + public: + explicit scoped_timer_t(cache_id id) : id_(id), start_(clock_::now()), active_(enabled()) {} + + ~scoped_timer_t() + { + if (!active_) { return; } + const double elapsed = + std::chrono::duration(clock_::now() - start_).count(); + add(id_, elapsed); + } + + private: + using clock_ = std::chrono::steady_clock; + cache_id id_; + clock_::time_point start_; + bool active_; +}; + +} // namespace cuopt::linear_programming::cache_profile + +#define CUOPT_CACHE_PROFILE_SCOPE(id) \ + ::cuopt::linear_programming::cache_profile::scoped_timer_t CUOPT_CACHE_PROFILE_CONCAT( \ + _cuopt_cache_scope_, __LINE__)(id) + +#define CUOPT_CACHE_PROFILE_CONCAT(a, b) CUOPT_CACHE_PROFILE_CONCAT_IMPL(a, b) +#define CUOPT_CACHE_PROFILE_CONCAT_IMPL(a, b) a##b diff --git a/cpp/src/barrier/CMakeLists.txt b/cpp/src/barrier/CMakeLists.txt index 650bc733e9..2d5fb27dc7 100644 --- a/cpp/src/barrier/CMakeLists.txt +++ b/cpp/src/barrier/CMakeLists.txt @@ -6,6 +6,7 @@ set(BARRIER_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/cusparse_view.cu ${CMAKE_CURRENT_SOURCE_DIR}/barrier.cu + ${CMAKE_CURRENT_SOURCE_DIR}/barrier_factorization_sparsity_hash.cu ${CMAKE_CURRENT_SOURCE_DIR}/device_sparse_matrix.cu ${CMAKE_CURRENT_SOURCE_DIR}/pinned_host_allocator.cu ) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index 1cbeed5cf3..830f37face 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -7,6 +7,8 @@ #include +#include +#include #include #include #include @@ -28,6 +30,10 @@ #include +#include + +#include + #include #include @@ -37,6 +43,7 @@ #include #include +#include #include #include @@ -211,7 +218,8 @@ class iteration_data_t { i_t num_upper_bounds, const std::vector& direct_free_variables, const csc_matrix_t& Qin, - const simplex_solver_settings_t& settings) + const simplex_solver_settings_t& settings, + barrier_symbolic_cache_t* adopt_symbolic = nullptr) : upper_bounds(num_upper_bounds), c(lp.objective), b(lp.rhs), @@ -246,7 +254,7 @@ class iteration_data_t { Q(Qin), cusparse_Q_view_(lp.handle_ptr, Q), cusparse_view_(lp.handle_ptr, lp.A), - cusparse_info(lp.handle_ptr), + cusparse_info_(nullptr), device_AD(lp.num_cols, lp.num_rows, 0, lp.handle_ptr->get_stream()), device_A(lp.num_cols, lp.num_rows, 0, lp.handle_ptr->get_stream()), device_ADAT(lp.num_rows, lp.num_rows, 0, lp.handle_ptr->get_stream()), @@ -324,12 +332,14 @@ class iteration_data_t { indefinite_Q(false), Q_diagonal(false), symbolic_status(0), + adopted_symbolic_(false), cone_combined_step_(false), cone_sigma_mu_(f_t(0)) { raft::common::nvtx::range fun_scope("Barrier: LP Data Creation"); { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C09); raft::common::nvtx::range scope("Barrier: LP Data: direct free linear"); // Setup tracking of direct free variables (linear columns only j < cone_start) n_direct_free_linear = direct_free_variables.size(); @@ -352,6 +362,7 @@ class iteration_data_t { bool has_Q = Q.x.size() > 0; indefinite_Q = false; { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C09); raft::common::nvtx::range scope("Barrier: LP Data: Q setup"); if (has_Q) { Qdiag.resize(lp.num_cols, 0.0); @@ -393,6 +404,7 @@ class iteration_data_t { } if (!lp.second_order_cone_dims.empty()) { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C08); raft::common::nvtx::range scope("Barrier: LP Data: SOC setup"); cone_var_start_ = lp.cone_var_start; i_t total_cone_dim = @@ -412,6 +424,7 @@ class iteration_data_t { } { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C09); raft::common::nvtx::range scope("Barrier: LP Data: complementarity buffers"); const i_t linear_xz_rhs_size = linear_xz_size(lp.num_cols); d_complementarity_xz_rhs_.resize(linear_xz_rhs_size, stream_view_); @@ -431,6 +444,7 @@ class iteration_data_t { } { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C09); raft::common::nvtx::range scope("Barrier: LP Data: upper bounds"); // Create the upper bounds vector n_upper_bounds = 0; @@ -444,6 +458,7 @@ class iteration_data_t { std::vector dense_columns_unordered; { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C04); raft::common::nvtx::range scope("Barrier: LP Data: dense columns and augmented"); // Decide if we are going to use the augmented system or not n_dense_columns = 0; @@ -515,6 +530,7 @@ class iteration_data_t { } { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C09); raft::common::nvtx::range scope("Barrier: LP Data: diag and inv_diag"); // D = I + EET diag.set_scalar(1.0); @@ -545,6 +561,7 @@ class iteration_data_t { if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C09); raft::common::nvtx::range scope("Barrier: LP Data: AD matrix setup"); // Copy A into AD AD = lp.A; @@ -580,6 +597,7 @@ class iteration_data_t { // device_AD / device_A / ADAT path is only used when forming ADAT (!use_augmented). if (!use_augmented) { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C09); raft::common::nvtx::range scope("Barrier: LP Data: device AD path"); device_AD.copy(AD, handle_ptr->get_stream()); d_original_A_values.resize(device_AD.x.size(), handle_ptr->get_stream()); @@ -598,36 +616,420 @@ class iteration_data_t { if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C06); raft::common::nvtx::range scope("Barrier: LP Data: Cholesky init"); i_t factorization_size = use_augmented ? lp.num_rows + lp.num_cols : lp.num_rows; - chol = std::make_unique>( - handle_ptr, settings, factorization_size); - chol->set_positive_definite(false); - } - if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } - { - raft::common::nvtx::range scope("Barrier: LP Data: symbolic analysis"); - // Perform symbolic analysis - symbolic_status = 0; - if (use_augmented) { - { - raft::common::nvtx::range form_scope("Barrier: LP Data: form augmented"); - // Build the sparsity pattern of the augmented system - form_augmented(true); + + auto adopt_augmented_symbolic = [&]() -> bool { + if (has_cones() || !use_augmented || adopt_symbolic == nullptr) { return false; } + + const barrier_sparsity_hash_t host_hash = hash_augmented_kkt_sparsity(A, AT, Q); + const bool matched = + adopt_symbolic->matches_reuse(host_hash, true, handle_ptr); + if (!matched) { return false; } + + chol = adopt_symbolic->chol; + static_cast*>(chol.get())->rebind_settings(settings); + pinned_device_augmented_ = &adopt_symbolic->device_augmented; + d_augmented_diagonal_indices_.resize(adopt_symbolic->d_augmented_diagonal_indices_.size(), + stream_view_); + raft::copy(d_augmented_diagonal_indices_.data(), + adopt_symbolic->d_augmented_diagonal_indices_.data(), + adopt_symbolic->d_augmented_diagonal_indices_.size(), + stream_view_); + handle_ptr->sync_stream(); + static_cast*>(chol.get())->rebind_csr_matrix(aug_mat()); + adopted_symbolic_ = true; + symbolic_status = 0; + return true; + }; + + auto unpin_adat_workspace = [&]() { + pinned_device_ADAT_ = nullptr; + pinned_device_A_ = nullptr; + pinned_device_AD_ = nullptr; + pinned_d_original_A_values_ = nullptr; + pinned_device_A_x_values_ = nullptr; + pinned_cusparse_info_ = nullptr; + }; + + auto pin_adat_from_cache = [&](barrier_symbolic_cache_t& cache) { + pinned_device_ADAT_ = &cache.device_ADAT; + pinned_device_A_ = &cache.device_A; + pinned_device_AD_ = &cache.device_AD; + pinned_d_original_A_values_ = &cache.d_original_A_values; + pinned_device_A_x_values_ = &cache.device_A_x_values; + pinned_cusparse_info_ = cache.cusparse_info.get(); + }; + + auto adopt_adat_symbolic = [&]() -> bool { + if (has_cones() || use_augmented || adopt_symbolic == nullptr || n_dense_columns > 0) { + return false; + } + if (!adopt_symbolic->valid || adopt_symbolic->use_augmented) { return false; } + + pin_adat_from_cache(*adopt_symbolic); + form_adat(true); + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + unpin_adat_workspace(); + return false; + } + + const barrier_sparsity_hash_t adat_hash = + hash_device_csr_sparsity_pattern(adat_mat(), stream_view_); + const bool matched = + adopt_symbolic->matches_reuse(adat_hash, false, handle_ptr); + if (!matched) { + unpin_adat_workspace(); + adopt_symbolic->clear(); + return false; + } + + chol = adopt_symbolic->chol; + static_cast*>(chol.get())->rebind_settings(settings); + handle_ptr->sync_stream(); + static_cast*>(chol.get())->rebind_csr_matrix(adat_mat()); + adopted_symbolic_ = true; + symbolic_status = 0; + return true; + }; + + if (!adopt_augmented_symbolic() && !adopt_adat_symbolic()) { + if (use_augmented) { + { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C05); + raft::common::nvtx::range form_scope("Barrier: LP Data: form augmented"); + form_augmented(true); + } + } else { + { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C05); + raft::common::nvtx::range form_scope("Barrier: LP Data: form ADAT"); + form_adat(true); + } } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } - symbolic_status = chol->analyze(device_augmented); - } else { + + chol = std::make_shared>( + handle_ptr, settings, factorization_size); + chol->set_positive_definite(false); + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } + symbolic_status = 0; { - raft::common::nvtx::range form_scope("Barrier: LP Data: form ADAT"); - form_adat(true); + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C07); + raft::common::nvtx::range analyze_scope("Barrier: LP Data: symbolic analysis"); + if (use_augmented) { + symbolic_status = chol->analyze(aug_mat()); + } else { + symbolic_status = chol->analyze(adat_mat()); + } } - if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } - symbolic_status = chol->analyze(device_ADAT); } } } + [[nodiscard]] bool adopted_symbolic() const { return adopted_symbolic_; } + + device_csr_matrix_t& augmented_system() { return aug_mat(); } + const device_csr_matrix_t& augmented_system() const { return aug_mat(); } + + void store_symbolic_cache(barrier_symbolic_cache_t& cache) + { + if (symbolic_status != 0 || has_cones()) { return; } + auto* cudss_chol = dynamic_cast*>(chol.get()); + if (cudss_chol == nullptr) { return; } + + cache.chol = std::static_pointer_cast>(chol); + cache.handle_ptr = handle_ptr; + + if (use_augmented) { + cache.cusparse_info.reset(); + cache.use_augmented = true; + + if (pinned_device_augmented_ != nullptr) { + // Warm reuse: sparsity_hash unchanged since adopt (values-only refresh); unpin only. + pinned_device_augmented_ = nullptr; + } else { + cache.sparsity_hash = + hash_device_csr_sparsity_pattern(device_augmented, handle_ptr->get_stream()); + cache.device_augmented = std::move(device_augmented); + cache.d_augmented_diagonal_indices_ = std::move(d_augmented_diagonal_indices_); + } + +#ifndef NDEBUG + const barrier_sparsity_hash_t host_hash = hash_augmented_kkt_sparsity(A, AT, Q); + cuopt_assert(cache.sparsity_hash == host_hash, + "store_symbolic_cache: device/host augmented sparsity hash mismatch"); +#endif + + cache.valid = true; + settings_.log.printf( + "Barrier: stored augmented symbolic cache hash=0x%016llx\n", + static_cast(cache.sparsity_hash)); + return; + } + + if (n_dense_columns > 0) { return; } + + cache.use_augmented = false; + if (pinned_device_ADAT_ != nullptr) { + // Warm reuse: sparsity_hash unchanged since adopt (values-only refresh); unpin only. + pinned_device_ADAT_ = nullptr; + pinned_device_A_ = nullptr; + pinned_device_AD_ = nullptr; + pinned_d_original_A_values_ = nullptr; + pinned_device_A_x_values_ = nullptr; + pinned_cusparse_info_ = nullptr; + } else { + cache.sparsity_hash = + hash_device_csr_sparsity_pattern(device_ADAT, handle_ptr->get_stream()); + cache.device_ADAT = std::move(device_ADAT); + cache.device_A = std::move(device_A); + cache.device_AD = std::move(device_AD); + cache.d_original_A_values = std::move(d_original_A_values); + cache.device_A_x_values = std::move(device_A_x_values); + cache.cusparse_info = std::move(cusparse_info_); + } + + cache.valid = true; + settings_.log.printf( + "Barrier: stored ADAT symbolic cache hash=0x%016llx\n", + static_cast(cache.sparsity_hash)); + } + + bool refresh_augmented_values() + { + i_t n = A.n; + i_t m = A.m; + i_t nnzA = A.col_start[n]; + i_t nnzQ = Q.n > 0 ? Q.col_start[n] : 0; + + i_t new_nnz = 2 * nnzA + n + m + nnzQ; + csr_matrix_t augmented_CSR(n + m, n + m, new_nnz); + i_t q = 0; + i_t off_diag_Qnz = 0; + + for (i_t i = 0; i < n; i++) { + augmented_CSR.row_start[i] = q; + if (nnzQ == 0) { + augmented_CSR.j[q] = i; + augmented_CSR.x[q++] = -diag[i] - dual_perturb; + } else { + const i_t q_col_beg = Q.col_start[i]; + const i_t q_col_end = Q.col_start[i + 1]; + bool has_diagonal = false; + for (i_t p = q_col_beg; p < q_col_end; ++p) { + augmented_CSR.j[q] = Q.i[p]; + if (Q.i[p] == i) { + has_diagonal = true; + augmented_CSR.x[q++] = -Q.x[p] - diag[i] - dual_perturb; + } else { + off_diag_Qnz++; + augmented_CSR.x[q++] = -Q.x[p]; + } + } + if (!has_diagonal) { + augmented_CSR.j[q] = i; + augmented_CSR.x[q++] = -diag[i] - dual_perturb; + } + } + const i_t col_beg = A.col_start[i]; + const i_t col_end = A.col_start[i + 1]; + for (i_t p = col_beg; p < col_end; ++p) { + augmented_CSR.j[q] = A.i[p] + n; + augmented_CSR.x[q++] = A.x[p]; + } + } + + for (i_t k = n; k < n + m; ++k) { + augmented_CSR.row_start[k] = q; + const i_t l = k - n; + const i_t col_beg = AT.col_start[l]; + const i_t col_end = AT.col_start[l + 1]; + for (i_t p = col_beg; p < col_end; ++p) { + augmented_CSR.j[q] = AT.i[p]; + augmented_CSR.x[q++] = AT.x[p]; + } + augmented_CSR.j[q] = k; + augmented_CSR.x[q++] = primal_perturb; + } + augmented_CSR.row_start[n + m] = q; + if (q != static_cast(aug_mat().x.size()) || q != 2 * nnzA + n + m + off_diag_Qnz) { + return false; + } + + augmented_CSR.j.resize(q); + augmented_CSR.x.resize(q); + raft::copy(aug_mat().x.data(), augmented_CSR.x.data(), q, handle_ptr->get_stream()); + RAFT_CHECK_CUDA(handle_ptr->get_stream()); + return true; + } + + bool rebuild_augmented_symbolic() + { + if (!use_augmented) { return false; } + + settings_.log.printf( + "Barrier: augmented nnz mismatch on cached symbolic; rebuilding symbolic analysis\n"); + + adopted_symbolic_ = false; + pinned_device_augmented_ = nullptr; + + const i_t factorization_size = A.n + A.m; + chol = std::make_shared>( + handle_ptr, settings_, factorization_size); + chol->set_positive_definite(false); + if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { return false; } + + form_augmented(true); + if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { return false; } + + symbolic_status = chol->analyze(aug_mat()); + if (symbolic_status != 0) { return false; } + + reset_for_new_solve(); + return true; + } + + bool refresh_adat_values() + { + if (use_augmented || n_dense_columns > 0) { return false; } + + const i_t expected_nnz = static_cast(adat_mat().x.size()); + form_adat(false); + handle_ptr->sync_stream(); + return static_cast(adat_mat().x.size()) == expected_nnz; + } + + bool rebuild_adat_symbolic() + { + if (use_augmented || n_dense_columns > 0) { return false; } + + settings_.log.printf( + "Barrier: ADAT nnz mismatch on cached symbolic; rebuilding symbolic analysis\n"); + + adopted_symbolic_ = false; + pinned_device_ADAT_ = nullptr; + pinned_device_A_ = nullptr; + pinned_device_AD_ = nullptr; + pinned_d_original_A_values_ = nullptr; + pinned_device_A_x_values_ = nullptr; + pinned_cusparse_info_ = nullptr; + cusparse_info_.reset(); + + const i_t factorization_size = A.m; + chol = std::make_shared>( + handle_ptr, settings_, factorization_size); + chol->set_positive_definite(false); + if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { return false; } + + form_adat(true); + if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { return false; } + + symbolic_status = chol->analyze(adat_mat()); + if (symbolic_status != 0) { return false; } + + reset_for_new_solve(); + return true; + } + + bool refresh_lp_numerics(const lp_problem_t& lp) + { + raft::common::nvtx::range fun_scope("Barrier: refresh LP numerics"); + + c = lp.objective; + b = lp.rhs; + + AD = A; + if (!use_augmented && n_dense_columns > 0) { AD.remove_columns(cols_to_remove); } + AT.transpose(AD); + + const bool has_Q = Q.n > 0; + if (has_Q) { + for (i_t j = 0; j < Q.n; j++) { + Qdiag[j] = 0.0; + const i_t col_start = Q.col_start[j]; + const i_t col_end = Q.col_start[j + 1]; + for (i_t p = col_start; p < col_end; p++) { + const i_t row = Q.i[p]; + if (j == row) { + Qdiag[j] = Q.x[p]; + break; + } + } + } + if (d_Q_diag_.size() > 0) { + raft::copy(d_Q_diag_.data(), Qdiag.data(), Qdiag.size(), stream_view_); + } + } + + diag.set_scalar(1.0); + if (n_upper_bounds > 0) { + for (i_t k = 0; k < n_upper_bounds; k++) { + const i_t j = upper_bounds[k]; + diag[j] = 2.0; + } + } + if (has_Q && !use_augmented) { + for (i_t j = 0; j < Q.n; j++) { + diag[j] += Qdiag[j]; + } + } + + inv_diag.set_scalar(1.0); + if (use_augmented) { diag.multiply_scalar(-1.0); } + if (n_upper_bounds > 0 || (has_Q && !use_augmented)) { diag.inverse(inv_diag); } + raft::copy(d_inv_diag.data(), inv_diag.data(), inv_diag.size(), stream_view_); + inv_sqrt_diag.set_scalar(1.0); + if (n_upper_bounds > 0 || (has_Q && !use_augmented)) { inv_diag.sqrt(inv_sqrt_diag); } + + if (!use_augmented) { + ad_mat().copy(AD, handle_ptr->get_stream()); + raft::copy(original_a_values().data(), + ad_mat().x.data(), + ad_mat().x.size(), + handle_ptr->get_stream()); + raft::copy(a_x_values().data(), ad_mat().x.data(), ad_mat().x.size(), handle_ptr->get_stream()); + ad_mat().to_compressed_row(a_mat(), handle_ptr->get_stream()); + RAFT_CHECK_CUDA(handle_ptr->get_stream()); + + if (adopted_symbolic_) { + if (!refresh_adat_values()) { + if (!rebuild_adat_symbolic()) { return false; } + } else { + handle_ptr->sync_stream(); + if (chol != nullptr) { chol->rebind_csr_matrix(adat_mat()); } + } + } + } + + if (use_augmented) { + if (!refresh_augmented_values()) { + if (!rebuild_augmented_symbolic()) { return false; } + } + } + + cusparse_view_.update_matrix_values(A); + if (Q.n > 0) { cusparse_Q_view_.update_matrix_values(Q); } + + reset_for_new_solve(); + return true; + } + + void reset_for_new_solve() + { + has_factorization = false; + has_solve_info = false; + relative_primal_residual_save = inf; + relative_dual_residual_save = inf; + relative_complementarity_residual_save = inf; + primal_residual_norm_save = inf; + dual_residual_norm_save = inf; + complementarity_residual_norm_save = inf; + if (chol != nullptr) { chol->invalidate_numeric_factor(); } + handle_ptr->sync_stream(); + } + bool has_cones() const { return cones_.has_value(); } cone_data_t& cones() @@ -887,7 +1289,7 @@ class iteration_data_t { thrust::for_each_n(rmm::exec_policy(handle_ptr->get_stream()), thrust::make_counting_iterator(0), linear_n, - [span_x = cuopt::make_span(device_augmented.x), + [span_x = cuopt::make_span(aug_mat().x), span_diag_indices = cuopt::make_span(d_augmented_diagonal_indices_), span_q_diag = cuopt::make_span(d_Q_diag_), span_diag = cuopt::make_span(d_diag_), @@ -901,7 +1303,7 @@ class iteration_data_t { thrust::for_each_n(rmm::exec_policy(handle_ptr->get_stream()), thrust::make_counting_iterator(n), i_t(m), - [span_x = cuopt::make_span(device_augmented.x), + [span_x = cuopt::make_span(aug_mat().x), span_diag_indices = cuopt::make_span(d_augmented_diagonal_indices_), primal_perturb_value = primal_perturb] __device__(i_t j) { span_x[span_diag_indices[j]] = primal_perturb_value; @@ -910,7 +1312,7 @@ class iteration_data_t { if (has_soc) { scatter_hessian_into_augmented(cones(), - device_augmented.x, + aug_mat().x, d_cone_csr_indices_, d_cone_Q_values_, handle_ptr->get_stream(), @@ -929,9 +1331,9 @@ class iteration_data_t { { raft::common::nvtx::range scope("Barrier: Form ADAT: restore A"); - raft::copy(device_AD.x.data(), - d_original_A_values.data(), - d_original_A_values.size(), + raft::copy(ad_mat().x.data(), + original_a_values().data(), + original_a_values().size(), handle_ptr->get_stream()); } { @@ -963,20 +1365,23 @@ class iteration_data_t { raft::common::nvtx::range scope("Barrier: Form ADAT: scale AD"); thrust::for_each_n(rmm::exec_policy(stream_view_), thrust::make_counting_iterator(0), - i_t(device_AD.x.size()), - [span_x = cuopt::make_span(device_AD.x), + i_t(ad_mat().x.size()), + [span_x = cuopt::make_span(ad_mat().x), span_scale = cuopt::make_span(d_inv_diag_prime), - span_col_ind = cuopt::make_span(device_AD.col_index)] __device__(i_t i) { + span_col_ind = cuopt::make_span(ad_mat().col_index)] __device__(i_t i) { span_x[i] *= span_scale[span_col_ind[i]]; }); RAFT_CHECK_CUDA(stream_view_); } if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { return; } - if (first_call) { + if (first_call && pinned_cusparse_info_ == nullptr) { raft::common::nvtx::range scope("Barrier: Form ADAT: cusparse init"); try { + if (!cusparse_info_) { + cusparse_info_ = std::make_unique>(handle_ptr); + } initialize_cusparse_data( - handle_ptr, device_A, device_AD, device_ADAT, cusparse_info); + handle_ptr, a_mat(), ad_mat(), adat_mat(), spgemm_info()); } catch (const raft::cuda_error& e) { settings_.log.printf("Error in initialize_cusparse_data: %s\n", e.what()); return; @@ -986,11 +1391,11 @@ class iteration_data_t { { raft::common::nvtx::range scope("Barrier: Form ADAT: ADAT multiply"); - multiply_kernels(handle_ptr, device_A, device_AD, device_ADAT, cusparse_info); + multiply_kernels(handle_ptr, a_mat(), ad_mat(), adat_mat(), spgemm_info()); handle_ptr->sync_stream(); } - auto adat_nnz = device_ADAT.row_start.element(device_ADAT.m, handle_ptr->get_stream()); + auto adat_nnz = adat_mat().row_start.element(adat_mat().m, handle_ptr->get_stream()); float64_t adat_time = toc(start_form_adat); if (num_factorizations == 0) { @@ -1000,7 +1405,7 @@ class iteration_data_t { settings_.log.printf( "ADAT density : %.2f\n", static_cast(adat_nnz) / - (static_cast(device_ADAT.m) * static_cast(device_ADAT.m))); + (static_cast(adat_mat().m) * static_cast(adat_mat().m))); } } @@ -1944,6 +2349,53 @@ class iteration_data_t { rmm::device_uvector d_augmented_diagonal_indices_; rmm::device_uvector d_cone_csr_indices_; rmm::device_uvector d_cone_Q_values_; + device_csr_matrix_t* pinned_device_augmented_{nullptr}; + device_csr_matrix_t* pinned_device_ADAT_{nullptr}; + device_csr_matrix_t* pinned_device_A_{nullptr}; + device_csc_matrix_t* pinned_device_AD_{nullptr}; + rmm::device_uvector* pinned_d_original_A_values_{nullptr}; + rmm::device_uvector* pinned_device_A_x_values_{nullptr}; + cusparse_info_t* pinned_cusparse_info_{nullptr}; + + device_csr_matrix_t& aug_mat() + { + return pinned_device_augmented_ != nullptr ? *pinned_device_augmented_ : device_augmented; + } + const device_csr_matrix_t& aug_mat() const + { + return pinned_device_augmented_ != nullptr ? *pinned_device_augmented_ : device_augmented; + } + device_csr_matrix_t& adat_mat() + { + return pinned_device_ADAT_ != nullptr ? *pinned_device_ADAT_ : device_ADAT; + } + const device_csr_matrix_t& adat_mat() const + { + return pinned_device_ADAT_ != nullptr ? *pinned_device_ADAT_ : device_ADAT; + } + device_csr_matrix_t& a_mat() + { + return pinned_device_A_ != nullptr ? *pinned_device_A_ : device_A; + } + device_csc_matrix_t& ad_mat() + { + return pinned_device_AD_ != nullptr ? *pinned_device_AD_ : device_AD; + } + rmm::device_uvector& original_a_values() + { + return pinned_d_original_A_values_ != nullptr ? *pinned_d_original_A_values_ + : d_original_A_values; + } + rmm::device_uvector& a_x_values() + { + return pinned_device_A_x_values_ != nullptr ? *pinned_device_A_x_values_ : device_A_x_values; + } + cusparse_info_t& spgemm_info() + { + cuopt_assert(pinned_cusparse_info_ != nullptr || cusparse_info_ != nullptr, + "spgemm_info: cusparse workspace unset"); + return pinned_cusparse_info_ != nullptr ? *pinned_cusparse_info_ : *cusparse_info_; + } bool indefinite_Q; cusparse_view_t cusparse_Q_view_; @@ -1952,6 +2404,7 @@ class iteration_data_t { bool use_augmented; i_t symbolic_status; + bool adopted_symbolic_; i_t n_direct_free_linear{0}; rmm::device_uvector d_is_direct_free_linear_; // 1 if variable is free in the linear block, else 0 @@ -1960,13 +2413,13 @@ class iteration_data_t { f_t dual_perturb{1e-8}; f_t primal_perturb{1e-8}; - std::unique_ptr> chol; + std::shared_ptr> chol; bool has_factorization; bool has_solve_info; i_t num_factorizations; - cusparse_info_t cusparse_info; + std::unique_ptr> cusparse_info_; cusparse_view_t cusparse_view_; detail::cusparse_dn_vec_descr_wrapper_t cusparse_tmp4_; detail::cusparse_dn_vec_descr_wrapper_t cusparse_h_; @@ -2198,13 +2651,13 @@ int barrier_solver_t::initial_point(iteration_data_t& data) // Perform a numerical factorization i_t status; if (use_augmented) { - status = data.chol->factorize(data.device_augmented); + status = data.chol->factorize(data.aug_mat()); #ifdef CHOLESKY_DEBUG_CHECK cholesky_debug_check(data, lp, use_augmented); #endif } else { - status = data.chol->factorize(data.device_ADAT); + status = data.chol->factorize(data.adat_mat()); } if (status == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; } if (status != 0) { @@ -2761,7 +3214,7 @@ i_t barrier_solver_t::gpu_compute_search_direction(iteration_data_tfactorize(data.device_augmented); + status = data.chol->factorize(data.aug_mat()); #ifdef CHOLESKY_DEBUG_CHECK cholesky_debug_check(data, lp, use_augmented); @@ -2776,7 +3229,7 @@ i_t barrier_solver_t::gpu_compute_search_direction(iteration_data_tfactorize(data.device_ADAT); + status = data.chol->factorize(data.adat_mat()); } data.has_factorization = true; data.num_factorizations++; @@ -3995,7 +4448,9 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( } template -lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t& solution) +lp_status_t barrier_solver_t::solve(f_t start_time, + lp_solution_t& solution, + cuopt::cython::lp_solve_session_t* session) { settings.log.printf("Barrier solver started at %.2f seconds\n", toc(start_time)); try { @@ -4031,18 +4486,57 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t Q(lp.num_cols, 0, 0); + std::unique_ptr> owned_data; + + auto finish_session = [&](lp_status_t status) -> lp_status_t { + if (session == nullptr) { return status; } + if (owned_data) { + if (status == lp_status_t::OPTIMAL) { + session->store_symbolic_cache(*owned_data); + } else { + session->clear_symbolic_cache(); + } + } + return status; + }; + if (lp.Q.n > 0) { create_Q(lp, Q); } + barrier_symbolic_cache_t* adopt_cache = nullptr; + if (session != nullptr) { adopt_cache = session->symbolic_cache_for_reuse(lp.handle_ptr); } + owned_data = std::make_unique>( + lp, num_upper_bounds, presolve_info.direct_free_variables, Q, settings, adopt_cache); + iteration_data_t& data = *owned_data; - iteration_data_t data( - lp, num_upper_bounds, presolve_info.direct_free_variables, Q, settings); + if (data.adopted_symbolic()) { + try { + if (!data.refresh_lp_numerics(lp)) { + settings.log.printf( + "Barrier: hash match but numeric refresh and symbolic rebuild failed\n"); + if (session != nullptr) { session->clear_symbolic_cache(); } + return finish_session(lp_status_t::NUMERICAL_ISSUES); + } + if (data.adopted_symbolic()) { + settings.log.printf("Barrier: reusing cuDSS symbolic analysis (sparsity hash match)\n"); + } else { + settings.log.printf( + "Barrier: rebuilt cuDSS symbolic analysis (%s nnz mismatch)\n", + data.use_augmented ? "augmented" : "adat"); + } + } catch (const raft::cuda_error&) { + settings.log.printf( + "Barrier: hash match but numeric refresh failed (CUDA); clearing symbolic cache\n"); + if (session != nullptr) { session->clear_symbolic_cache(); } + return finish_session(lp_status_t::NUMERICAL_ISSUES); + } + } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Barrier solver halted\n"); - return lp_status_t::CONCURRENT_LIMIT; + return finish_session(lp_status_t::CONCURRENT_LIMIT); } - if (data.indefinite_Q) { return lp_status_t::NUMERICAL_ISSUES; } + if (data.indefinite_Q) { return finish_session(lp_status_t::NUMERICAL_ISSUES); } if (data.symbolic_status != 0) { settings.log.printf("Error in symbolic analysis\n"); - return lp_status_t::NUMERICAL_ISSUES; + return finish_session(lp_status_t::NUMERICAL_ISSUES); } data.cusparse_dual_residual_ = data.cusparse_view_.create_vector(data.d_dual_residual_); @@ -4056,21 +4550,21 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t settings.time_limit) { settings.log.printf("Barrier time limit exceeded\n"); - return lp_status_t::TIME_LIMIT; + return finish_session(lp_status_t::TIME_LIMIT); } i_t initial_status = initial_point(data); if (toc(start_time) > settings.time_limit) { settings.log.printf("Barrier time limit exceeded\n"); - return lp_status_t::TIME_LIMIT; + return finish_session(lp_status_t::TIME_LIMIT); } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Barrier solver halted\n"); - return lp_status_t::CONCURRENT_LIMIT; + return finish_session(lp_status_t::CONCURRENT_LIMIT); } if (initial_status != 0) { settings.log.printf("Unable to compute initial point\n"); - return lp_status_t::NUMERICAL_ISSUES; + return finish_session(lp_status_t::NUMERICAL_ISSUES); } // Upload initial point to device and compute initial residuals/norms on GPU data.d_complementarity_wv_residual_.resize(data.n_upper_bounds, stream_view_); @@ -4175,11 +4669,11 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t settings.time_limit) { settings.log.printf("Barrier time limit exceeded\n"); - return lp_status_t::TIME_LIMIT; + return finish_session(lp_status_t::TIME_LIMIT); } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Barrier solver halted\n"); - return lp_status_t::CONCURRENT_LIMIT; + return finish_session(lp_status_t::CONCURRENT_LIMIT); } // Compute the affine step. This is the call that (re)factorizes the @@ -4192,7 +4686,7 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t settings.time_limit) { settings.log.printf("Barrier time limit exceeded\n"); - return lp_status_t::TIME_LIMIT; + return finish_session(lp_status_t::TIME_LIMIT); } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Barrier solver halted\n"); - return lp_status_t::CONCURRENT_LIMIT; + return finish_session(lp_status_t::CONCURRENT_LIMIT); } f_t mu_aff, sigma, new_mu; @@ -4241,30 +4735,30 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t settings.time_limit) { settings.log.printf("Barrier time limit exceeded\n"); - return lp_status_t::TIME_LIMIT; + return finish_session(lp_status_t::TIME_LIMIT); } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Barrier solver halted\n"); - return lp_status_t::CONCURRENT_LIMIT; + return finish_session(lp_status_t::CONCURRENT_LIMIT); } compute_final_direction(data); @@ -4331,17 +4825,17 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t::solve(f_t start_time, lp_solution_t +void barrier_store_symbolic_cache_from_iteration_data(iteration_data_t& data, + barrier_symbolic_cache_t& cache) +{ + data.store_symbolic_cache(cache); +} + #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE template bool validate_barrier_cone_layout( const lp_problem_t& problem, const simplex_solver_settings_t& settings); @@ -4446,6 +4947,9 @@ template class barrier_solver_t; template class sparse_cholesky_base_t; template class sparse_cholesky_cudss_t; template class iteration_data_t; + +template void barrier_store_symbolic_cache_from_iteration_data( + iteration_data_t& data, barrier_symbolic_cache_t& cache); #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/barrier/barrier.hpp b/cpp/src/barrier/barrier.hpp index b17cce3979..626c05dd36 100644 --- a/cpp/src/barrier/barrier.hpp +++ b/cpp/src/barrier/barrier.hpp @@ -16,6 +16,11 @@ #include #include + +namespace cuopt::cython { +class lp_solve_session_t; +} // namespace cuopt::cython + namespace cuopt::linear_programming::dual_simplex { /** Validates SOC layout on an lp_problem_t before barrier presolve/solve. */ @@ -32,7 +37,9 @@ class barrier_solver_t { barrier_solver_t(const lp_problem_t& lp, const presolve_info_t& presolve, const simplex_solver_settings_t& settings); - lp_status_t solve(f_t start_time, lp_solution_t& solution); + lp_status_t solve(f_t start_time, + lp_solution_t& solution, + cuopt::cython::lp_solve_session_t* session = nullptr); private: void my_pop_range(bool debug) const; diff --git a/cpp/src/barrier/barrier_factorization_sparsity_hash.cu b/cpp/src/barrier/barrier_factorization_sparsity_hash.cu new file mode 100644 index 0000000000..44cb4c89bc --- /dev/null +++ b/cpp/src/barrier/barrier_factorization_sparsity_hash.cu @@ -0,0 +1,25 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include + +namespace cuopt::linear_programming::dual_simplex { + +template +barrier_sparsity_hash_t hash_device_csr_sparsity_pattern( + device_csr_matrix_t& mat, rmm::cuda_stream_view stream) +{ + const csr_matrix_t host = mat.to_host(stream); + return hash_host_csr_sparsity_pattern(host.m, host.row_start, host.j); +} + +template barrier_sparsity_hash_t hash_device_csr_sparsity_pattern( + device_csr_matrix_t&, rmm::cuda_stream_view); +template barrier_sparsity_hash_t hash_device_csr_sparsity_pattern( + device_csr_matrix_t&, rmm::cuda_stream_view); + +} // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/barrier/barrier_factorization_sparsity_hash.hpp b/cpp/src/barrier/barrier_factorization_sparsity_hash.hpp new file mode 100644 index 0000000000..5451c9a215 --- /dev/null +++ b/cpp/src/barrier/barrier_factorization_sparsity_hash.hpp @@ -0,0 +1,126 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include + +#include + +#include +#include + +namespace cuopt::linear_programming::dual_simplex { + +using barrier_sparsity_hash_t = std::uint64_t; + +/// FNV-1a style mix for incremental hashing. +inline barrier_sparsity_hash_t barrier_hash_combine(barrier_sparsity_hash_t h, std::uint64_t value) +{ + constexpr barrier_sparsity_hash_t kPrime = 1099511628211ULL; + h ^= value; + h *= kPrime; + return h; +} + +inline barrier_sparsity_hash_t barrier_hash_u64(std::uint64_t value) +{ + return barrier_hash_combine(1469598103934665603ULL, value); +} + +/** + * @brief Hash CSR sparsity (row_start + col indices); numeric values are ignored. + */ +template +barrier_sparsity_hash_t hash_host_csr_sparsity_pattern(i_t num_rows, + const std::vector& row_start, + const std::vector& col_indices) +{ + barrier_sparsity_hash_t h = barrier_hash_u64(static_cast(num_rows)); + h = barrier_hash_combine(h, static_cast(col_indices.size())); + for (i_t k = 0; k <= num_rows; ++k) { + h = barrier_hash_combine(h, static_cast(row_start[static_cast(k)])); + } + for (i_t col : col_indices) { + h = barrier_hash_combine(h, static_cast(col)); + } + return h; +} + +/** + * @brief Hash the sparsity pattern of the augmented KKT matrix passed to cuDSS (host CSR). + * + * Must match the index layout produced by iteration_data_t::form_augmented(true). + */ +template +barrier_sparsity_hash_t hash_augmented_kkt_sparsity(const csc_matrix_t& A, + const csc_matrix_t& AT, + const csc_matrix_t& Q) +{ + const i_t n = A.n; + const i_t m = A.m; + const i_t size = n + m; + + std::vector row_start(static_cast(size + 1), 0); + std::vector col_indices; + col_indices.reserve(static_cast(2) * static_cast(A.col_start[n]) + + static_cast(n + m) + + (Q.n > 0 ? static_cast(Q.col_start[n]) : 0)); + + i_t q = 0; + for (i_t i = 0; i < n; ++i) { + row_start[static_cast(i)] = q; + if (Q.n == 0) { + col_indices.push_back(i); + ++q; + } else { + const i_t q_col_beg = Q.col_start[i]; + const i_t q_col_end = Q.col_start[i + 1]; + bool has_diagonal = false; + for (i_t p = q_col_beg; p < q_col_end; ++p) { + col_indices.push_back(Q.i[p]); + ++q; + if (Q.i[p] == i) { has_diagonal = true; } + } + if (!has_diagonal) { + col_indices.push_back(i); + ++q; + } + } + const i_t col_beg = A.col_start[i]; + const i_t col_end = A.col_start[i + 1]; + for (i_t p = col_beg; p < col_end; ++p) { + col_indices.push_back(A.i[p] + n); + ++q; + } + } + + for (i_t k = n; k < n + m; ++k) { + row_start[static_cast(k)] = q; + const i_t l = k - n; + const i_t col_beg = AT.col_start[l]; + const i_t col_end = AT.col_start[l + 1]; + for (i_t p = col_beg; p < col_end; ++p) { + col_indices.push_back(AT.i[p]); + ++q; + } + col_indices.push_back(k); + ++q; + } + row_start[static_cast(size)] = q; + + return hash_host_csr_sparsity_pattern(size, row_start, col_indices); +} + +/** + * @brief Hash CSR sparsity from a device matrix (copies row/col indices to host). + */ +template +barrier_sparsity_hash_t hash_device_csr_sparsity_pattern( + device_csr_matrix_t& mat, rmm::cuda_stream_view stream); + +} // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/barrier/barrier_symbolic_cache.hpp b/cpp/src/barrier/barrier_symbolic_cache.hpp new file mode 100644 index 0000000000..1ba7bde2b6 --- /dev/null +++ b/cpp/src/barrier/barrier_symbolic_cache.hpp @@ -0,0 +1,82 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include +#include +#include + +#include + +#include +#include + +namespace cuopt::linear_programming::dual_simplex { + +using barrier_sparsity_hash_t = std::uint64_t; + +template +class sparse_cholesky_cudss_t; + +/** + * @brief Cached cuDSS symbolic state and GPU buffers for hash-gated barrier reuse. + * + * Holds reordering + symbolic factorization in @p chol, the sparsity hash of the matrix + * that was analyzed, and path-specific GPU workspace (augmented KKT or ADAT + cuSPARSE). + */ +template +struct barrier_symbolic_cache_t { + std::shared_ptr> chol; + barrier_sparsity_hash_t sparsity_hash{0}; + raft::handle_t const* handle_ptr{nullptr}; + bool use_augmented{false}; + bool valid{false}; + + // --- Augmented KKT (use_augmented == true) --- + device_csr_matrix_t device_augmented; + rmm::device_uvector d_augmented_diagonal_indices_; + + // --- ADAT (use_augmented == false) --- + device_csr_matrix_t device_ADAT; + device_csc_matrix_t device_AD; + device_csr_matrix_t device_A; + rmm::device_uvector d_original_A_values; + rmm::device_uvector device_A_x_values; + std::unique_ptr> cusparse_info; + + explicit barrier_symbolic_cache_t(rmm::cuda_stream_view stream) + : device_augmented(stream), + d_augmented_diagonal_indices_(0, stream), + device_ADAT(stream), + device_AD(stream), + device_A(stream), + d_original_A_values(0, stream), + device_A_x_values(0, stream) + { + } + + void clear() + { + chol.reset(); + sparsity_hash = 0; + handle_ptr = nullptr; + use_augmented = false; + valid = false; + cusparse_info.reset(); + } + + [[nodiscard]] bool matches_reuse(barrier_sparsity_hash_t hash, + bool augmented, + raft::handle_t const* handle) const + { + return valid && handle != nullptr && handle_ptr == handle && use_augmented == augmented && + sparsity_hash == hash; + } +}; + +} // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/barrier/cusparse_view.cu b/cpp/src/barrier/cusparse_view.cu index b7673eacd5..587597dea9 100644 --- a/cpp/src/barrier/cusparse_view.cu +++ b/cpp/src/barrier/cusparse_view.cu @@ -254,6 +254,16 @@ cusparse_view_t::~cusparse_view_t() CUOPT_CUSPARSE_TRY_NO_THROW(cusparseDestroySpMat(A_T_)); } +template +void cusparse_view_t::update_matrix_values(const csc_matrix_t& A) +{ + const auto stream = handle_ptr_->get_stream(); + raft::copy(A_T_data_.data(), A.x.data(), A.x.size(), stream); + csr_matrix_t A_csr(A.m, A.n, 1); + A.to_compressed_row(A_csr); + raft::copy(A_data_.data(), A_csr.x.data(), A_csr.x.size(), stream); +} + template detail::cusparse_dn_vec_descr_wrapper_t cusparse_view_t::create_vector( rmm::device_uvector const& vec) diff --git a/cpp/src/barrier/cusparse_view.hpp b/cpp/src/barrier/cusparse_view.hpp index 802fc90f8b..36ebd3a851 100644 --- a/cpp/src/barrier/cusparse_view.hpp +++ b/cpp/src/barrier/cusparse_view.hpp @@ -55,6 +55,8 @@ class cusparse_view_t { f_t beta, detail::cusparse_dn_vec_descr_wrapper_t const& y); + void update_matrix_values(const csc_matrix_t& A); + raft::handle_t const* handle_ptr_{nullptr}; private: diff --git a/cpp/src/barrier/device_sparse_matrix.cuh b/cpp/src/barrier/device_sparse_matrix.cuh index 29927c81b3..085be6c713 100644 --- a/cpp/src/barrier/device_sparse_matrix.cuh +++ b/cpp/src/barrier/device_sparse_matrix.cuh @@ -125,6 +125,10 @@ class device_csc_matrix_t { { } + device_csc_matrix_t(device_csc_matrix_t&&) = default; + device_csc_matrix_t& operator=(device_csc_matrix_t&&) = default; + device_csc_matrix_t& operator=(const device_csc_matrix_t&) = delete; + device_csc_matrix_t(const csc_matrix_t& A, rmm::cuda_stream_view stream) : m(A.m), n(A.n), @@ -252,6 +256,10 @@ class device_csr_matrix_t { { } + device_csr_matrix_t(device_csr_matrix_t&&) = default; + device_csr_matrix_t& operator=(device_csr_matrix_t&&) = default; + device_csr_matrix_t& operator=(const device_csr_matrix_t&) = delete; + device_csr_matrix_t(const csr_matrix_t& A, rmm::cuda_stream_view stream) : m(A.m), n(A.n), diff --git a/cpp/src/barrier/sparse_cholesky.cuh b/cpp/src/barrier/sparse_cholesky.cuh index c2223e5080..d056a55a33 100644 --- a/cpp/src/barrier/sparse_cholesky.cuh +++ b/cpp/src/barrier/sparse_cholesky.cuh @@ -33,6 +33,9 @@ class sparse_cholesky_base_t { virtual i_t solve(const dense_vector_t& b, dense_vector_t& x) = 0; virtual i_t solve(rmm::device_uvector& b, rmm::device_uvector& x) = 0; virtual void set_positive_definite(bool positive_definite) = 0; + virtual void invalidate_numeric_factor() {} + virtual void rebind_csr_matrix(device_csr_matrix_t& Arow) {} + virtual void rebind_settings(const simplex_solver_settings_t& settings) {} }; #define CUDSS_EXAMPLE_FREE \ @@ -143,7 +146,9 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { first_factor(true), positive_definite(true), A_created(false), - settings_(settings), + settings_(&settings), + symbolic_done_(false), + numeric_factor_valid_(false), stream(handle_ptr->get_stream()) { int major, minor, patch; @@ -155,8 +160,8 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { cuda_error = cudaSuccess; status = CUDSS_STATUS_SUCCESS; - if (CUDART_VERSION >= 13000 && settings_.concurrent_halt != nullptr && - settings_.num_gpus == 1) { + if (CUDART_VERSION >= 13000 && settings_->concurrent_halt != nullptr && + settings_->num_gpus == 1) { cuGetErrorString_func = cuopt::detail::get_driver_entry_point("cuGetErrorString"); // 1. Set up the GPU resources CUdevResource initial_device_GPU_resources = {}; @@ -279,18 +284,18 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { "cudssConfigSet for device count"); #if CUDSS_VERSION_MAJOR >= 0 && CUDSS_VERSION_MINOR >= 7 - if (settings_.concurrent_halt != nullptr) { + if (settings_->concurrent_halt != nullptr) { CUDSS_CALL_AND_CHECK_EXIT(cudssDataSet(handle, solverData, CUDSS_DATA_USER_HOST_INTERRUPT, - (void*)settings_.concurrent_halt, + (void*)settings_->concurrent_halt, sizeof(int)), status, "cudssDataSet for interrupt"); } - if (settings_.cudss_deterministic) { - settings_.log.printf("cuDSS solve mode : deterministic\n"); + if (settings_->cudss_deterministic) { + settings_->log.printf("cuDSS solve mode : deterministic\n"); int32_t deterministic = 1; CUDSS_CALL_AND_CHECK_EXIT( cudssConfigSet( @@ -309,7 +314,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { #endif #if USE_MATCHING - settings_.log.printf("Using matching\n"); + settings_->log.printf("Using matching\n"); int32_t use_matching = 1; CUDSS_CALL_AND_CHECK_EXIT( cudssConfigSet(solverConfig, CUDSS_CONFIG_USE_MATCHING, &use_matching, sizeof(int32_t)), @@ -373,7 +378,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { CUDSS_CALL_AND_CHECK_EXIT(cudssDestroy(handle), status, "cudssDestroy"); CUDA_CALL_AND_CHECK_EXIT(cudaStreamSynchronize(stream), "cudaStreamSynchronize"); #if CUDART_VERSION >= 13000 - if (settings_.concurrent_halt != nullptr && settings_.num_gpus == 1) { + if (settings_->concurrent_halt != nullptr && settings_->num_gpus == 1) { auto cuStreamDestroy_func = cuopt::detail::get_driver_entry_point("cuStreamDestroy"); CU_CHECK(reinterpret_cast(cuStreamDestroy_func)(stream), reinterpret_cast(cuGetErrorString_func)); @@ -396,9 +401,9 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { csc_matrix_t A_col(Arow_host.m, Arow_host.n, 1); Arow_host.to_compressed_col(A_col); FILE* fid = fopen("A_to_factorize.mtx", "w"); - settings_.log.printf("writing matrix matrix\n"); + settings_->log.printf("writing matrix matrix\n"); A_col.write_matrix_market(fid); - settings_.log.printf("finished\n"); + settings_->log.printf("finished\n"); fclose(fid); } #endif @@ -407,9 +412,9 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { const f_t density = static_cast(nnz) / (static_cast(n) * static_cast(n)); if (first_factor && - ((settings_.ordering == -1 && density >= 0.05 && nnz > n) || settings_.ordering == 1) && + ((settings_->ordering == -1 && density >= 0.05 && nnz > n) || settings_->ordering == 1) && n > 1) { - settings_.log.printf("Reordering algorithm : AMD\n"); + settings_->log.printf("Reordering algorithm : AMD\n"); // Tell cuDSS to use AMD #if CUDSS_VERSION_MAJOR > 0 || (CUDSS_VERSION_MAJOR == 0 && CUDSS_VERSION_MINOR >= 8) cudssReorderingAlg_t reorder_alg = CUDSS_REORDERING_ALG_AMD; @@ -482,27 +487,27 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { raft::common::nvtx::range fun_scope("Barrier: cuDSS Analyze : CUDSS_PHASE_ANALYSIS"); status = cudssExecute(handle, CUDSS_PHASE_REORDERING, solverConfig, solverData, A, cudss_x, cudss_b); - if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { + if (settings_->concurrent_halt != nullptr && *settings_->concurrent_halt == 1) { return CONCURRENT_HALT_RETURN; } if (status != CUDSS_STATUS_SUCCESS) { - settings_.log.printf( + settings_->log.printf( "FAILED: CUDSS call ended unsuccessfully with status = %d, details: cuDSSExecute for " "reordering\n", status); return -1; } f_t reordering_time = toc(start_symbolic); - settings_.log.printf("Reordering time : %.2fs\n", reordering_time); + settings_->log.printf("Reordering time : %.2fs\n", reordering_time); start_symbolic_factor = tic(); status = cudssExecute( handle, CUDSS_PHASE_SYMBOLIC_FACTORIZATION, solverConfig, solverData, A, cudss_x, cudss_b); - if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { + if (settings_->concurrent_halt != nullptr && *settings_->concurrent_halt == 1) { return CONCURRENT_HALT_RETURN; } if (status != CUDSS_STATUS_SUCCESS) { - settings_.log.printf( + settings_->log.printf( "FAILED: CUDSS call ended unsuccessfully with status = %d, details: cuDSSExecute for " "symbolic factorization\n", status); @@ -511,33 +516,43 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { } RAFT_CUDA_TRY(cudaStreamSynchronize(stream)); f_t symbolic_factorization_time = toc(start_symbolic_factor); - settings_.log.printf("Symbolic factorization time : %.2fs\n", symbolic_factorization_time); - settings_.log.printf("Total symbolic time : %.2fs\n", toc(start_symbolic)); + settings_->log.printf("Symbolic factorization time : %.2fs\n", symbolic_factorization_time); + settings_->log.printf("Total symbolic time : %.2fs\n", toc(start_symbolic)); int64_t lu_nz = 0; size_t size_written = 0; CUDSS_CALL_AND_CHECK( cudssDataGet(handle, solverData, CUDSS_DATA_LU_NNZ, &lu_nz, sizeof(int64_t), &size_written), status, "cudssDataGet for LU_NNZ"); - settings_.log.printf("Symbolic nonzeros in factor : %.2e\n", static_cast(lu_nz) / 2.0); + settings_->log.printf("Symbolic nonzeros in factor : %.2e\n", static_cast(lu_nz) / 2.0); // TODO: Is there any way to get nonzeros in the factors? // TODO: Is there any way to get flops for the factorization? RAFT_CUDA_TRY(cudaStreamSynchronize(stream)); handle_ptr_->get_stream().synchronize(); + symbolic_done_ = true; + numeric_factor_valid_ = false; return 0; } i_t factorize(device_csr_matrix_t& Arow) override { raft::common::nvtx::range fun_scope("Factorize: cuDSS"); + if (!symbolic_done_ || !A_created) { + settings_->log.printf( + "Error: cuDSS factorize(device_csr) called before analyze (symbolic_done=%d A_created=%d)\n", + static_cast(symbolic_done_), + static_cast(A_created)); + return -1; + } + // #define PRINT_MATRIX_NORM #ifdef PRINT_MATRIX_NORM cudaStreamSynchronize(stream); csr_matrix_t Arow_host = Arow.to_host(Arow.row_start.stream()); csc_matrix_t A_col(Arow_host.m, Arow_host.n, 1); Arow_host.to_compressed_col(A_col); - settings_.log.printf( + settings_->log.printf( "before factorize || A to factor|| = %.16e hash: %zu\n", A_col.norm1(), A_col.hash()); cudaStreamSynchronize(stream); #endif @@ -546,7 +561,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { auto d_nnz = Arow.row_start.element(Arow.m, Arow.row_start.stream()); if (nnz != d_nnz) { - settings_.log.printf("Error: nnz %d != A_in.col_start[A_in.n] %d\n", nnz, d_nnz); + settings_->log.printf("Error: nnz %d != A_in.col_start[A_in.n] %d\n", nnz, d_nnz); return -1; } @@ -556,11 +571,11 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { f_t start_numeric = tic(); status = cudssExecute( handle, CUDSS_PHASE_FACTORIZATION, solverConfig, solverData, A, cudss_x, cudss_b); - if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { + if (settings_->concurrent_halt != nullptr && *settings_->concurrent_halt == 1) { return CONCURRENT_HALT_RETURN; } if (status != CUDSS_STATUS_SUCCESS) { - settings_.log.printf( + settings_->log.printf( "FAILED: CUDSS call ended unsuccessfully with status = %d, details: cuDSSExecute for " "factorization\n", status); @@ -572,7 +587,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { #endif f_t numeric_time = toc(start_numeric); - if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { + if (settings_->concurrent_halt != nullptr && *settings_->concurrent_halt == 1) { return CONCURRENT_HALT_RETURN; } @@ -586,16 +601,16 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { handle_ptr_->get_stream().synchronize(); RAFT_CUDA_TRY(cudaStreamSynchronize(stream)); if (info != 0) { - settings_.log.printf("Factorization failed info %d\n", info); + settings_->log.printf("Factorization failed info %d\n", info); return -1; } if (first_factor) { - settings_.log.debug("Factorization time : %.2fs\n", numeric_time); + settings_->log.debug("Factorization time : %.2fs\n", numeric_time); first_factor = false; } if (status != CUDSS_STATUS_SUCCESS) { - settings_.log.printf("cuDSS Factorization failed\n"); + settings_->log.printf("cuDSS Factorization failed\n"); return -1; } return 0; @@ -609,15 +624,15 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { FILE* fid = fopen("A.mtx", "w"); A_in.write_matrix_market(fid); fclose(fid); - settings_.log.printf("Wrote A.mtx\n"); + settings_->log.printf("Wrote A.mtx\n"); #endif A_in.to_compressed_row(Arow); #ifdef CHECK_MATRIX - settings_.log.printf("Checking matrices\n"); + settings_->log.printf("Checking matrices\n"); A_in.check_matrix(); Arow.check_matrix(); - settings_.log.printf("Finished checking matrices\n"); + settings_->log.printf("Finished checking matrices\n"); #endif if (A_in.n != n) { printf("Analyze input does not match size %d != %d\n", A_in.n, n); @@ -691,7 +706,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { A_created = true; // Perform symbolic analysis - if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { + if (settings_->concurrent_halt != nullptr && *settings_->concurrent_halt == 1) { return CONCURRENT_HALT_RETURN; } f_t start_analysis = tic(); @@ -701,7 +716,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { "cudssExecute for reordering"); f_t reorder_time = toc(start_analysis); - if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { + if (settings_->concurrent_halt != nullptr && *settings_->concurrent_halt == 1) { return CONCURRENT_HALT_RETURN; } @@ -715,8 +730,8 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { f_t symbolic_time = toc(start_symbolic); f_t analysis_time = toc(start_analysis); - settings_.log.printf("Symbolic factorization time : %.2fs\n", symbolic_time); - if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { + settings_->log.printf("Symbolic factorization time : %.2fs\n", symbolic_time); + if (settings_->concurrent_halt != nullptr && *settings_->concurrent_halt == 1) { RAFT_CUDA_TRY(cudaStreamSynchronize(stream)); handle_ptr_->get_stream().synchronize(); return CONCURRENT_HALT_RETURN; @@ -727,7 +742,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { cudssDataGet(handle, solverData, CUDSS_DATA_LU_NNZ, &lu_nz, sizeof(int64_t), &size_written), status, "cudssDataGet for LU_NNZ"); - settings_.log.printf("Symbolic nonzeros in factor : %.2e\n", static_cast(lu_nz) / 2.0); + settings_->log.printf("Symbolic nonzeros in factor : %.2e\n", static_cast(lu_nz) / 2.0); RAFT_CUDA_TRY(cudaStreamSynchronize(stream)); handle_ptr_->get_stream().synchronize(); // TODO: Is there any way to get nonzeros in the factors? @@ -740,10 +755,10 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { csr_matrix_t Arow(A_in.n, A_in.m, A_in.col_start[A_in.n]); A_in.to_compressed_row(Arow); - if (A_in.n != n) { settings_.log.printf("Error A in n %d != size %d\n", A_in.n, n); } + if (A_in.n != n) { settings_->log.printf("Error A in n %d != size %d\n", A_in.n, n); } if (nnz != A_in.col_start[A_in.n]) { - settings_.log.printf( + settings_->log.printf( "Error: nnz %d != A_in.col_start[A_in.n] %d\n", nnz, A_in.col_start[A_in.n]); return -1; } @@ -767,7 +782,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { "cudssExecute for factorization"); f_t numeric_time = toc(start_numeric); - if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { + if (settings_->concurrent_halt != nullptr && *settings_->concurrent_halt == 1) { return CONCURRENT_HALT_RETURN; } @@ -780,16 +795,16 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { RAFT_CUDA_TRY(cudaStreamSynchronize(stream)); handle_ptr_->get_stream().synchronize(); if (info != 0) { - settings_.log.printf("Factorization failed info %d\n", info); + settings_->log.printf("Factorization failed info %d\n", info); return -1; } if (first_factor) { - settings_.log.debug("Factorization time : %.2fs\n", numeric_time); + settings_->log.debug("Factorization time : %.2fs\n", numeric_time); first_factor = false; } if (status != CUDSS_STATUS_SUCCESS) { - settings_.log.printf("cuDSS Factorization failed\n"); + settings_->log.printf("cuDSS Factorization failed\n"); return -1; } return 0; @@ -818,11 +833,11 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { { handle_ptr_->get_stream().synchronize(); if (static_cast(b.size()) != n) { - settings_.log.printf("Error: b.size() %d != n %d\n", b.size(), n); + settings_->log.printf("Error: b.size() %d != n %d\n", b.size(), n); return -1; } if (static_cast(x.size()) != n) { - settings_.log.printf("Error: x.size() %d != n %d\n", x.size(), n); + settings_->log.printf("Error: x.size() %d != n %d\n", x.size(), n); return -1; } @@ -832,11 +847,11 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { cudssMatrixSetValues(cudss_x, x.data()), status, "cudssMatrixSetValues for x"); status = cudssExecute(handle, CUDSS_PHASE_SOLVE, solverConfig, solverData, A, cudss_x, cudss_b); - if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { + if (settings_->concurrent_halt != nullptr && *settings_->concurrent_halt == 1) { return CONCURRENT_HALT_RETURN; } if (status != CUDSS_STATUS_SUCCESS) { - settings_.log.printf( + settings_->log.printf( "FAILED: CUDSS call ended unsuccessfully with status = %d, details: cuDSSExecute for " "solve\n", status); @@ -852,7 +867,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { raft::copy(b_host.data(), b.data(), n, stream); raft::copy(x_host.data(), x.data(), n, stream); cudaStreamSynchronize(stream); - settings_.log.printf("RHS norm %.16e, hash: %zu, Solution norm %.16e, hash: %zu\n", + settings_->log.printf("RHS norm %.16e, hash: %zu, Solution norm %.16e, hash: %zu\n", vector_norm2(b_host), compute_hash(b_host), vector_norm2(x_host), @@ -867,6 +882,63 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { this->positive_definite = positive_definite; } + void rebind_settings(const simplex_solver_settings_t& settings) override + { + settings_ = &settings; + } + + void invalidate_numeric_factor() override { numeric_factor_valid_ = false; } + + /// Re-point cuDSS CSR wrapper at current device buffers after in-place value refresh. + void rebind_csr_matrix(device_csr_matrix_t& Arow) override + { + if (!symbolic_done_ || !A_created) { return; } + auto d_nnz = Arow.row_start.element(Arow.m, Arow.row_start.stream()); + if (d_nnz != nnz) { return; } + status = cudssMatrixDestroy(A); + if (status != CUDSS_STATUS_SUCCESS) { + settings_->log.printf("cudssMatrixDestroy for A rebind failed: %d\n", status); + return; + } +#if CUDSS_VERSION_MAJOR > 0 || (CUDSS_VERSION_MAJOR == 0 && CUDSS_VERSION_MINOR >= 8) + status = cudssMatrixCreateCsr(&A, + n, + n, + nnz, + Arow.row_start.data(), + nullptr, + Arow.j.data(), + Arow.x.data(), + CUDSS_R_32I, + CUDSS_R_32I, + CUDSS_R_64F, + positive_definite ? CUDSS_MTYPE_SPD : CUDSS_MTYPE_SYMMETRIC, + CUDSS_MVIEW_FULL, + CUDSS_BASE_ZERO); +#else + status = cudssMatrixCreateCsr(&A, + n, + n, + nnz, + Arow.row_start.data(), + nullptr, + Arow.j.data(), + Arow.x.data(), + CUDA_R_32I, + CUDA_R_64F, + positive_definite ? CUDSS_MTYPE_SPD : CUDSS_MTYPE_SYMMETRIC, + CUDSS_MVIEW_FULL, + CUDSS_BASE_ZERO); +#endif + if (status != CUDSS_STATUS_SUCCESS) { + settings_->log.printf("cudssMatrixCreateCsr rebind failed: %d\n", status); + A_created = false; + return; + } + A_created = true; + numeric_factor_valid_ = false; + } + private: raft::handle_t const* handle_ptr_; i_t n; @@ -890,7 +962,9 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { f_t* x_values_d; f_t* b_values_d; - const simplex_solver_settings_t& settings_; + const simplex_solver_settings_t* settings_; + bool symbolic_done_; + bool numeric_factor_valid_; CUgreenCtx barrier_green_ctx; CUstream stream; void* cuGetErrorString_func; diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index d81265358a..ab30fcc2a0 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -351,7 +351,8 @@ template lp_status_t solve_linear_program_with_barrier(const user_problem_t& user_problem, const simplex_solver_settings_t& settings, f_t start_time, - lp_solution_t& solution) + lp_solution_t& solution, + cuopt::cython::lp_solve_session_t* session) { lp_status_t status = lp_status_t::UNSET; lp_problem_t original_lp(user_problem.handle_ptr, 1, 1, 1); @@ -388,7 +389,7 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us lp_solution_t barrier_solution(barrier_lp.num_rows, barrier_lp.num_cols); barrier_solver_t barrier_solver(barrier_lp, presolve_info, barrier_settings); - lp_status_t barrier_status = barrier_solver.solve(start_time, barrier_solution); + lp_status_t barrier_status = barrier_solver.solve(start_time, barrier_solution, session); if (barrier_status == lp_status_t::OPTIMAL) { #ifdef COMPUTE_SCALED_RESIDUALS std::vector scaled_residual = barrier_lp.rhs; @@ -672,10 +673,11 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us template lp_status_t solve_linear_program_with_barrier(const user_problem_t& user_problem, const simplex_solver_settings_t& settings, - lp_solution_t& solution) + lp_solution_t& solution, + cuopt::cython::lp_solve_session_t* session) { f_t start_time = tic(); - return solve_linear_program_with_barrier(user_problem, settings, start_time, solution); + return solve_linear_program_with_barrier(user_problem, settings, start_time, solution, session); } template @@ -817,13 +819,15 @@ template lp_status_t solve_linear_program_with_advanced_basis( template lp_status_t solve_linear_program_with_barrier( const user_problem_t& user_problem, const simplex_solver_settings_t& settings, - lp_solution_t& solution); + lp_solution_t& solution, + cuopt::cython::lp_solve_session_t* session); template lp_status_t solve_linear_program_with_barrier( const user_problem_t& user_problem, const simplex_solver_settings_t& settings, double start_time, - lp_solution_t& solution); + lp_solution_t& solution, + cuopt::cython::lp_solve_session_t* session); template lp_status_t solve_linear_program(const user_problem_t& user_problem, const simplex_solver_settings_t& settings, diff --git a/cpp/src/dual_simplex/solve.hpp b/cpp/src/dual_simplex/solve.hpp index 84eb6faa02..3ad5db550b 100644 --- a/cpp/src/dual_simplex/solve.hpp +++ b/cpp/src/dual_simplex/solve.hpp @@ -17,6 +17,10 @@ namespace cuopt { struct work_limit_context_t; } +namespace cuopt::cython { +class lp_solve_session_t; +} // namespace cuopt::cython + namespace cuopt::linear_programming::dual_simplex { template @@ -88,15 +92,19 @@ lp_status_t solve_linear_program_with_advanced_basis( work_limit_context_t* work_unit_context = nullptr); template -lp_status_t solve_linear_program_with_barrier(const user_problem_t& user_problem, - const simplex_solver_settings_t& settings, - lp_solution_t& solution); +lp_status_t solve_linear_program_with_barrier( + const user_problem_t& user_problem, + const simplex_solver_settings_t& settings, + lp_solution_t& solution, + cuopt::cython::lp_solve_session_t* session = nullptr); template -lp_status_t solve_linear_program_with_barrier(const user_problem_t& user_problem, - const simplex_solver_settings_t& settings, - f_t start_time, - lp_solution_t& solution); +lp_status_t solve_linear_program_with_barrier( + const user_problem_t& user_problem, + const simplex_solver_settings_t& settings, + f_t start_time, + lp_solution_t& solution, + cuopt::cython::lp_solve_session_t* session = nullptr); template lp_status_t solve_linear_program(const user_problem_t& user_problem, diff --git a/cpp/src/pdlp/CMakeLists.txt b/cpp/src/pdlp/CMakeLists.txt index f5f26837b6..7d0d984d5a 100644 --- a/cpp/src/pdlp/CMakeLists.txt +++ b/cpp/src/pdlp/CMakeLists.txt @@ -34,6 +34,7 @@ set(LP_CORE_FILES # C and Python adapter files set(LP_ADAPTER_FILES ${CMAKE_CURRENT_SOURCE_DIR}/utilities/cython_solve.cu + ${CMAKE_CURRENT_SOURCE_DIR}/utilities/lp_solve_session.cu ${CMAKE_CURRENT_SOURCE_DIR}/cuopt_c.cpp ) diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 5d29e97fa0..5b73c2ee10 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -34,6 +34,7 @@ #include #include #include +#include #include #include @@ -61,6 +62,7 @@ #include #include +#include #include #include #include @@ -70,6 +72,52 @@ namespace cuopt::linear_programming { +namespace { + +template +uint64_t fnv1a64_mix(uint64_t hash, uint64_t value) +{ + constexpr uint64_t kFnvPrime = 1099511628211ULL; + constexpr uint64_t kFnvOffset = 14695981039346656037ULL; + if (hash == 0) { hash = kFnvOffset; } + for (int shift = 0; shift < 64; shift += 8) { + hash ^= (value >> shift) & 0xFFULL; + hash *= kFnvPrime; + } + return hash; +} + +template +uint64_t compute_problem_fingerprint(const optimization_problem_t& op) +{ + uint64_t hash = fnv1a64_mix(0, static_cast(op.get_n_variables())); + hash = fnv1a64_mix(hash, static_cast(op.get_n_constraints())); + hash = fnv1a64_mix(hash, static_cast(op.get_nnz())); + + const auto offsets = op.get_constraint_matrix_offsets_host(); + for (i_t off : offsets) { + hash = fnv1a64_mix(hash, static_cast(off)); + } + const auto indices = op.get_constraint_matrix_indices_host(); + for (i_t idx : indices) { + hash = fnv1a64_mix(hash, static_cast(idx)); + } + + if (op.has_quadratic_objective()) { + const auto q_offsets = op.get_quadratic_objective_offsets(); + for (i_t off : q_offsets) { + hash = fnv1a64_mix(hash, static_cast(off)); + } + const auto q_indices = op.get_quadratic_objective_indices(); + for (i_t idx : q_indices) { + hash = fnv1a64_mix(hash, static_cast(idx)); + } + } + return hash; +} + +} // namespace + template extern rmm::device_uvector gpu_cast(const rmm::device_uvector& src, rmm::cuda_stream_view stream); @@ -491,7 +539,8 @@ template std::tuple, dual_simplex::lp_status_t, f_t, f_t, f_t> run_barrier(dual_simplex::user_problem_t& user_problem, pdlp_solver_settings_t const& settings, - const timer_t& timer) + const timer_t& timer, + cuopt::cython::lp_solve_session_t* session = nullptr) { f_t norm_user_objective = dual_simplex::vector_norm2(user_problem.objective); f_t norm_rhs = dual_simplex::vector_norm2(user_problem.rhs); @@ -523,7 +572,7 @@ run_barrier(dual_simplex::user_problem_t& user_problem, dual_simplex::lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); auto status = dual_simplex::solve_linear_program_with_barrier( - user_problem, barrier_settings, timer.get_tic_start(), solution); + user_problem, barrier_settings, timer.get_tic_start(), solution, session); detail::project_barrier_solution_to_model_variables(user_problem, solution); @@ -546,12 +595,13 @@ template optimization_problem_solution_t run_barrier( detail::problem_t& problem, pdlp_solver_settings_t const& settings, - const timer_t& timer) + const timer_t& timer, + cuopt::cython::lp_solve_session_t* session = nullptr) { // Convert data structures to dual simplex format and back dual_simplex::user_problem_t dual_simplex_problem = cuopt_problem_to_user_problem(problem.handle_ptr, problem); - auto sol_dual_simplex = run_barrier(dual_simplex_problem, settings, timer); + auto sol_dual_simplex = run_barrier(dual_simplex_problem, settings, timer, session); return convert_dual_simplex_sol(problem, std::get<0>(sol_dual_simplex), std::get<1>(sol_dual_simplex), @@ -1765,7 +1815,7 @@ optimization_problem_solution_t solve_lp_with_method( if (settings.method == method_t::DualSimplex) { return run_dual_simplex(problem, settings, timer); } else if (settings.method == method_t::Barrier) { - return run_barrier(problem, settings, timer); + return run_barrier(problem, settings, timer, settings.lp_solve_session); } else if (settings.method == method_t::Concurrent) { return run_concurrent(problem, settings, timer, is_batch_mode); } else { @@ -1793,7 +1843,10 @@ optimization_problem_solution_t solve_qcqp( print_version_info(); // Init libraries before to not include it in solve time - init_handler(op_problem.get_handle_ptr()); + { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C02); + init_handler(op_problem.get_handle_ptr()); + } auto qcqp_timer = cuopt::timer_t(settings.time_limit); @@ -1831,10 +1884,14 @@ optimization_problem_solution_t solve_qcqp( CUOPT_LOG_INFO("Writing user problem to file: %s", settings.user_problem_file.c_str()); op_problem.write_to_mps(settings.user_problem_file); } + { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C03); + [[maybe_unused]] const uint64_t fingerprint = compute_problem_fingerprint(op_problem); + } // Convert data structures to dual simplex format and back dual_simplex::user_problem_t dual_simplex_problem = cuopt_optimization_problem_to_user_problem(op_problem.get_handle_ptr(), op_problem); - auto sol_dual_simplex = run_barrier(dual_simplex_problem, settings, qcqp_timer); + auto sol_dual_simplex = run_barrier(dual_simplex_problem, settings, qcqp_timer, settings.lp_solve_session); auto solution = convert_dual_simplex_sol(op_problem, std::get<0>(sol_dual_simplex), std::get<1>(sol_dual_simplex), @@ -1894,7 +1951,10 @@ optimization_problem_solution_t solve_lp( // Init libraries before to not include it in solve time // This needs to be called before pdlp is initialized - init_handler(op_problem.get_handle_ptr()); + { + CUOPT_CACHE_PROFILE_SCOPE(cuopt::linear_programming::cache_profile::cache_id::C02); + init_handler(op_problem.get_handle_ptr()); + } raft::common::nvtx::range fun_scope("Running solver"); diff --git a/cpp/src/pdlp/utilities/cython_solve.cu b/cpp/src/pdlp/utilities/cython_solve.cu index ec3cc2862f..4e49366853 100644 --- a/cpp/src/pdlp/utilities/cython_solve.cu +++ b/cpp/src/pdlp/utilities/cython_solve.cu @@ -18,6 +18,8 @@ #include #include #include +#include +#include #include #include #include @@ -30,11 +32,25 @@ #include #include +#include + #include namespace cuopt { namespace cython { +namespace { + +bool uses_barrier_session_path( + cuopt::linear_programming::solver_settings_t& solver_settings, + cuopt::linear_programming::io::data_model_view_t const& data_model) +{ + if (data_model.has_quadratic_objective() || data_model.has_quadratic_constraints()) { return true; } + return solver_settings.get_pdlp_settings().method == cuopt::linear_programming::method_t::Barrier; +} + +} // namespace + /** * @brief Wrapper for linear_programming to expose the API to cython * @@ -94,24 +110,68 @@ std::unique_ptr call_solve( cuopt::linear_programming::io::data_model_view_t* data_model, cuopt::linear_programming::solver_settings_t* solver_settings, unsigned int flags, - bool is_batch_mode) + bool is_batch_mode, + lp_solve_session_t* session_in) { raft::common::nvtx::range fun_scope("Call Solve"); + namespace cache_profile = cuopt::linear_programming::cache_profile; + if (cache_profile::enabled()) { cache_profile::reset(); } + + cuopt_expects(data_model != nullptr, + error_type_t::ValidationError, + "call_solve: data_model is null."); + cuopt_expects(solver_settings != nullptr, + error_type_t::ValidationError, + "call_solve: solver_settings is null."); + // Determine memory backend based on execution mode auto memory_backend = cuopt::linear_programming::get_memory_backend_type(); solver_ret_t response; + auto& pdlp_settings = solver_settings->get_pdlp_settings(); + const bool session_enabled = pdlp_settings.session_enabled; + const bool barrier_path = uses_barrier_session_path(*solver_settings, *data_model); + const bool want_session = (session_in != nullptr || session_enabled) && barrier_path && + memory_backend == cuopt::linear_programming::memory_backend_t::GPU && + !is_batch_mode; + + std::unique_ptr owned_session; + lp_solve_session_t* active_session = session_in; + pdlp_settings.lp_solve_session = nullptr; + + rmm::cuda_stream ephemeral_stream(static_cast(flags)); + raft::handle_t ephemeral_handle(ephemeral_stream); + raft::handle_t* solve_handle = &ephemeral_handle; + // Create problem instance and CUDA resources based on memory backend if (memory_backend == cuopt::linear_programming::memory_backend_t::GPU) { - // GPU memory backend: Create CUDA resources and GPU problem - rmm::cuda_stream stream(static_cast(flags)); - const raft::handle_t handle_{stream}; + if (want_session) { + if (active_session == nullptr) { + const auto handle_start = std::chrono::steady_clock::now(); + owned_session = lp_solve_session_t::create(flags); + active_session = owned_session.get(); + if (cache_profile::enabled()) { + const double elapsed = + std::chrono::duration(std::chrono::steady_clock::now() - handle_start).count(); + cache_profile::add(cache_profile::cache_id::C01, elapsed); + } + } + solve_handle = active_session->handle_ptr(); + pdlp_settings.lp_solve_session = active_session; + } else { + const auto handle_start = std::chrono::steady_clock::now(); + if (cache_profile::enabled()) { + const double elapsed = + std::chrono::duration(std::chrono::steady_clock::now() - handle_start).count(); + cache_profile::add(cache_profile::cache_id::C01, elapsed); + } + } - auto problem = cuopt::linear_programming::optimization_problem_t(&handle_); + auto problem = cuopt::linear_programming::optimization_problem_t(solve_handle); cuopt::linear_programming::populate_from_data_model_view( - &problem, data_model, solver_settings, &handle_); + &problem, data_model, solver_settings, solve_handle); // Call appropriate solve function and convert to ret struct if (problem.get_problem_category() == linear_programming::problem_category_t::LP) { @@ -140,6 +200,8 @@ std::unique_ptr call_solve( gpu_sols.last_restart_duality_gap_primal_solution_->set_stream(rmm::cuda_stream_per_thread); gpu_sols.last_restart_duality_gap_dual_solution_->set_stream(rmm::cuda_stream_per_thread); + if (owned_session) { response.lp_ret.lp_solve_session = std::move(owned_session); } + } else { // MIP solve auto mip_solution_ptr = @@ -198,6 +260,10 @@ std::unique_ptr call_solve( } } + if (cache_profile::enabled()) { cache_profile::log_summary(); } + + pdlp_settings.lp_solve_session = nullptr; + return std::make_unique(std::move(response)); } @@ -285,7 +351,7 @@ std::pair>, double> call_batch_solve( #pragma omp parallel for num_threads(max_thread) for (std::size_t i = 0; i < size; ++i) - list[i] = call_solve(data_models[i], solver_settings, cudaStreamNonBlocking, is_batch_mode); + list[i] = call_solve(data_models[i], solver_settings, cudaStreamNonBlocking, is_batch_mode, nullptr); auto end = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration_cast(end - start_solver); diff --git a/cpp/src/pdlp/utilities/lp_solve_session.cu b/cpp/src/pdlp/utilities/lp_solve_session.cu new file mode 100644 index 0000000000..ad44d1e926 --- /dev/null +++ b/cpp/src/pdlp/utilities/lp_solve_session.cu @@ -0,0 +1,83 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include + +#include + +#include +#include + +namespace cuopt::cython { + +struct lp_solve_session_t::impl { + std::unique_ptr stream; + std::unique_ptr handle; + std::optional> + symbolic_cache; +}; + +lp_solve_session_t::lp_solve_session_t(std::unique_ptr stream, + std::unique_ptr handle) + : impl_(std::make_unique(impl{std::move(stream), std::move(handle), std::nullopt})) +{ +} + +lp_solve_session_t::~lp_solve_session_t() = default; + +lp_solve_session_t::lp_solve_session_t(lp_solve_session_t&&) noexcept = default; +lp_solve_session_t& lp_solve_session_t::operator=(lp_solve_session_t&&) noexcept = default; + +std::unique_ptr lp_solve_session_t::create(unsigned stream_flags) +{ + auto stream = std::make_unique(static_cast(stream_flags)); + auto handle = std::make_unique(*stream); + return std::unique_ptr( + new lp_solve_session_t(std::move(stream), std::move(handle))); +} + +raft::handle_t* lp_solve_session_t::handle_ptr() +{ + return impl_->handle.get(); +} + +raft::handle_t const* lp_solve_session_t::handle_ptr() const +{ + return impl_->handle.get(); +} + +rmm::cuda_stream_view lp_solve_session_t::stream_view() const +{ + return impl_->stream->view(); +} + +linear_programming::dual_simplex::barrier_symbolic_cache_t* +lp_solve_session_t::symbolic_cache_for_reuse(raft::handle_t const* handle) +{ + if (handle == nullptr || !impl_->symbolic_cache.has_value() || !impl_->symbolic_cache->valid || + impl_->symbolic_cache->handle_ptr != handle) { + return nullptr; + } + return &(*impl_->symbolic_cache); +} + +void lp_solve_session_t::clear_symbolic_cache() +{ + impl_->symbolic_cache.reset(); +} + +void lp_solve_session_t::store_symbolic_cache( + linear_programming::dual_simplex::iteration_data_t& data) +{ + if (!impl_->symbolic_cache.has_value()) { + impl_->symbolic_cache.emplace(impl_->handle->get_stream()); + } + linear_programming::dual_simplex::barrier_store_symbolic_cache_from_iteration_data( + data, *impl_->symbolic_cache); +} + +} // namespace cuopt::cython diff --git a/python/cuopt/cuopt/linear_programming/problem.py b/python/cuopt/cuopt/linear_programming/problem.py index 10600a543b..ccf6256a59 100644 --- a/python/cuopt/cuopt/linear_programming/problem.py +++ b/python/cuopt/cuopt/linear_programming/problem.py @@ -7,7 +7,7 @@ from enum import Enum import numpy as np -from scipy.sparse import coo_matrix +from scipy.sparse import coo_matrix, csr_matrix import cuopt.linear_programming.data_model as data_model from cuopt.linear_programming import ParseMps, Read @@ -1387,9 +1387,10 @@ def getCoefficient(self, var): v_idx = var.index return self.vindex_coeff_dict[v_idx] - def compute_slack(self): + def compute_slack(self, index_to_var=None): # Computes the constraint Slack in the current solution. - index_to_var = {var.index: var for var in self.vars} + if index_to_var is None: + index_to_var = {var.index: var for var in self.vars} lhs = sum( index_to_var[v_idx].Value * coeff for v_idx, coeff in self.vindex_coeff_dict.items() @@ -1459,6 +1460,7 @@ def __init__(self, model_name=""): self.model = None self.solved = False + self._session = None self.rhs = None self.row_sense = None self.constraint_csr_matrix = None @@ -1466,6 +1468,8 @@ def __init__(self, model_name=""): self.lower_bound = None self.upper_bound = None self.var_type = None + self._index_to_var_cache = None + self._constraint_csr_scipy = None class dict_to_object: def __init__(self, mdict): @@ -1557,6 +1561,7 @@ def _to_data_model(self): constr_name = "R" + str(constr.index) self.row_names.append(constr_name) self.constraint_csr_matrix = csr_dict + self._constraint_csr_scipy = None else: for constr in self.constrs: @@ -1591,6 +1596,7 @@ def _to_data_model(self): if self.ObjSense == -1: dm.set_maximize(True) dm.set_constraint_bounds(np.array(self.rhs)) + self.rhs = np.asarray(self.rhs, dtype=np.float64) dm.set_row_types(np.array(self.row_sense, dtype="S1")) dm.set_objective_coefficients(self.objective) dm.set_objective_offset(self.ObjConstant) @@ -1629,6 +1635,92 @@ def _to_data_model(self): self.model = dm + def _index_to_var(self): + if ( + self._index_to_var_cache is None + or len(self._index_to_var_cache) != len(self.vars) + ): + self._index_to_var_cache = {var.index: var for var in self.vars} + return self._index_to_var_cache + + def _invalidate_index_to_var_cache(self): + self._index_to_var_cache = None + self._constraint_csr_scipy = None + + def _constraint_csr_scipy_matrix(self): + csr_dict = self.constraint_csr_matrix + if csr_dict is None or self.rhs is None: + return None + if self._constraint_csr_scipy is not None: + return self._constraint_csr_scipy + n_rows = ( + len(self.rhs) + if isinstance(self.rhs, np.ndarray) + else len(self.rhs) + ) + self._constraint_csr_scipy = csr_matrix( + ( + np.asarray(csr_dict["values"], dtype=np.float64), + np.asarray(csr_dict["column_indices"], dtype=np.int32), + np.asarray(csr_dict["row_pointers"], dtype=np.int32), + ), + shape=(n_rows, len(self.vars)), + ) + return self._constraint_csr_scipy + + def _refresh_data_model_values(self): + """Patch existing DataModel when sparsity structure is unchanged.""" + n = len(self.vars) + if ( + self.model is None + or self.constraint_csr_matrix is None + or self.objective is None + or len(self.objective) != n + or self.rhs is None + ): + self._to_data_model() + return + + for j in range(n): + self.objective[j] = self.vars[j].getObjectiveCoefficient() + + rhs_arr = ( + self.rhs + if isinstance(self.rhs, np.ndarray) + else np.asarray(self.rhs, dtype=np.float64) + ) + linear_row = 0 + for constr in self.constrs: + if constr.is_quadratic: + continue + rhs_arr[linear_row] = constr.RHS + linear_row += 1 + self.rhs = rhs_arr + + self.model.set_objective_coefficients(self.objective) + self.model.set_constraint_bounds(rhs_arr) + + def _populate_slacks_vectorized(self, primal_sol): + """Assign constraint slacks via CSR matvec (structure-stable problems).""" + A = self._constraint_csr_scipy_matrix() + if A is None: + return False + + rhs_arr = ( + self.rhs + if isinstance(self.rhs, np.ndarray) + else np.asarray(self.rhs, dtype=np.float64) + ) + slacks = rhs_arr - A.dot(primal_sol) + + linear_row = 0 + for constr in self.constrs: + if constr.is_quadratic: + continue + constr.Slack = slacks[linear_row] + linear_row += 1 + return True + def update(self): """ Update the problem. This is mandatory if attributes of @@ -1651,7 +1743,9 @@ def reset_solved_values(self): self.constraint_csr_matrix = None self.objective_qmatrix = None self.warmstart_data = None + self._session = None self.solved = False + self._invalidate_index_to_var_cache() def addVariable( self, lb=0.0, ub=float("inf"), obj=0.0, vtype=CONTINUOUS, name="" @@ -1685,6 +1779,10 @@ def addVariable( """ if self.solved: self.reset_solved_values() # Reset all solved values + else: + self.constraint_csr_matrix = None + self.model = None + self._invalidate_index_to_var_cache() n = len(self.vars) var = Variable(lb, ub, obj, vtype, name) var.index = n @@ -1716,6 +1814,10 @@ def addConstraint(self, constr, name=""): """ if self.solved: self.reset_solved_values() # Reset all solved values + else: + self.constraint_csr_matrix = None + self.model = None + self._invalidate_index_to_var_cache() n = len(self.constrs) match constr: case Constraint(): @@ -2186,21 +2288,44 @@ def populate_solution(self, solution): dual_sol = None if not IsMIP: dual_sol = solution.get_dual_solution() - linear_row = 0 - for constr in self.constrs: - if constr.is_quadratic: - continue - if dual_sol is not None and len(dual_sol) > linear_row: - constr.DualValue = dual_sol[linear_row] - constr.Slack = constr.compute_slack() - linear_row += 1 + if not ( + not IsMIP + and len(primal_sol) > 0 + and self._populate_slacks_vectorized(primal_sol) + ): + index_to_var = self._index_to_var() + linear_row = 0 + for constr in self.constrs: + if constr.is_quadratic: + continue + if dual_sol is not None and len(dual_sol) > linear_row: + constr.DualValue = dual_sol[linear_row] + constr.Slack = constr.compute_slack(index_to_var) + linear_row += 1 + else: + linear_row = 0 + for constr in self.constrs: + if constr.is_quadratic: + continue + if dual_sol is not None and len(dual_sol) > linear_row: + constr.DualValue = dual_sol[linear_row] + linear_row += 1 self.solved = True - def solve(self, settings=solver_settings.SolverSettings()): + def solve(self, settings=solver_settings.SolverSettings(), session=None): """ Optimizes the LP or MIP problem with the added variables, constraints and objective. + Parameters + ---------- + settings : SolverSettings + Solver configuration. + session : PyCapsule, optional + Reuse ``cuopt.lp_solve_session`` from a prior solve. When omitted and + ``settings.session_enabled`` is true, the problem keeps the session from + the last solve until the structure changes. + Examples -------- >>> problem = problem.Problem("MIP_model") @@ -2212,10 +2337,13 @@ def solve(self, settings=solver_settings.SolverSettings()): >>> problem.setObjective(x + y, sense=MAXIMIZE) >>> problem.solve() """ - if self.model is None: + if self.model is None or self.constraint_csr_matrix is None: self._to_data_model() - # Call Solver - solution = solver.Solve(self.model, settings) - # Post Solve + else: + self._refresh_data_model_values() + active_session = session if session is not None else self._session + solution = solver.Solve(self.model, settings, session=active_session) + if getattr(settings, "session_enabled", False) and solution.lp_solve_session is not None: + self._session = solution.lp_solve_session self.populate_solution(solution) return solution diff --git a/python/cuopt/cuopt/linear_programming/solution/solution.py b/python/cuopt/cuopt/linear_programming/solution/solution.py index 93d224fdbd..bb3470ca50 100644 --- a/python/cuopt/cuopt/linear_programming/solution/solution.py +++ b/python/cuopt/cuopt/linear_programming/solution/solution.py @@ -123,6 +123,11 @@ class Solution: Note: Applicable to only LP Whether the LP was solved by Dual Simplex, PDLP or Barrier. This is populated by the solver using the values from SolverMethod. + lp_solve_session: optional + GPU barrier/QCQP only: ``PyCapsule`` owning ``cuopt::cython::lp_solve_session_t*`` + (name ``cuopt.lp_solve_session``). Pass to a subsequent solve to reuse the GPU handle + and cuDSS symbolic analysis when sparsity is unchanged. ``None`` on reuse solves or when + session persistence was not enabled. """ def __init__( @@ -168,6 +173,7 @@ def __init__( max_variable_bound_violation=0.0, num_nodes=0, num_simplex_iterations=0, + lp_solve_session=None, ): self.problem_category = problem_category self.primal_solution = primal_solution @@ -210,6 +216,7 @@ def __init__( "nb_iterations": nb_iterations, } self.reduced_cost = reduced_cost + self.lp_solve_session = lp_solve_session self.milp_stats = { "mip_gap": mip_gap, "solution_bound": solution_bound, diff --git a/python/cuopt/cuopt/linear_programming/solver/solver.pxd b/python/cuopt/cuopt/linear_programming/solver/solver.pxd index 5570af1a04..079403c379 100644 --- a/python/cuopt/cuopt/linear_programming/solver/solver.pxd +++ b/python/cuopt/cuopt/linear_programming/solver/solver.pxd @@ -93,6 +93,10 @@ cdef extern from "cuopt/linear_programming/utilities/cython_types.hpp" namespace vector[double] last_restart_duality_gap_primal_solution_ vector[double] last_restart_duality_gap_dual_solution_ +cdef extern from "cuopt/linear_programming/utilities/lp_solve_session.hpp" namespace "cuopt::cython": # noqa + cdef cppclass lp_solve_session_t: + pass + cdef extern from "cuopt/linear_programming/utilities/cython_solve.hpp" namespace "cuopt::cython": # noqa # Unified LP solution struct — solutions_ variant accessed via helpers cdef cppclass linear_programming_ret_t: @@ -117,6 +121,7 @@ cdef extern from "cuopt/linear_programming/utilities/cython_solve.hpp" namespace int nb_iterations_ double solve_time_ method_t solved_by_ + unique_ptr[lp_solve_session_t] lp_solve_session bool is_gpu() # Unified MIP solution struct — solution_ variant accessed via helpers @@ -144,6 +149,9 @@ cdef extern from "cuopt/linear_programming/utilities/cython_solve.hpp" namespace cdef unique_ptr[solver_ret_t] call_solve( data_model_view_t[int, double]* data_model, solver_settings_t[int, double]* solver_settings, + unsigned int flags, + bool is_batch_mode, + lp_solve_session_t* session_in, ) except + nogil cdef pair[vector[unique_ptr[solver_ret_t]], double] call_batch_solve( # noqa diff --git a/python/cuopt/cuopt/linear_programming/solver/solver.py b/python/cuopt/cuopt/linear_programming/solver/solver.py index 4963d0be8a..caec39b538 100644 --- a/python/cuopt/cuopt/linear_programming/solver/solver.py +++ b/python/cuopt/cuopt/linear_programming/solver/solver.py @@ -11,7 +11,7 @@ @catch_cuopt_exception -def Solve(data_model, solver_settings=None): +def Solve(data_model, solver_settings=None, session=None): """ Solve the Linear Program passed as input and returns the solution. @@ -35,6 +35,9 @@ def Solve(data_model, solver_settings=None): solver_settings: SolverSettings Settings to configure solver configurations. By default, it uses default solver settings to solve. + session: PyCapsule, optional + Reuse capsule ``cuopt.lp_solve_session`` from a prior barrier/QCQP solve + (``Solution.lp_solve_session``). Returns ------- @@ -103,6 +106,7 @@ def is_mip(var_types): data_model, solver_settings, mip=is_mip(data_model.get_variable_types()), + session=session, ) if emit_stamps: print(f"CUOPT_SOLVE_RETURN: {time.time()}") diff --git a/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx index 925dd8797b..64a6ffc317 100644 --- a/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx +++ b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx @@ -15,7 +15,7 @@ from dateutil.relativedelta import relativedelta from cuopt.utilities import type_cast -from libc.stdint cimport uintptr_t +from libc.stdint cimport uintptr_t, uint32_t from libc.stdlib cimport free, malloc from libc.string cimport memcpy, strcpy, strlen from libcpp cimport bool @@ -25,6 +25,13 @@ from libcpp.string cimport string from libcpp.utility cimport move from libcpp.vector cimport vector +from cpython.pycapsule cimport ( + PyCapsule_Destructor, + PyCapsule_GetPointer, + PyCapsule_IsValid, + PyCapsule_New, +) + from rmm.pylibrmm.device_buffer cimport DeviceBuffer from cuopt.linear_programming.data_model.data_model cimport data_model_view_t @@ -40,6 +47,7 @@ from cuopt.linear_programming.solver.solver cimport ( linear_programming_ret_t, lp_cpu_solutions_t, lp_gpu_solutions_t, + lp_solve_session_t, mip_ret_t, mip_termination_status_t, pdlp_solver_mode_t, @@ -76,6 +84,25 @@ cdef extern from "cuopt/linear_programming/utilities/internals.hpp" namespace "c cdef cppclass base_solution_callback_t +cdef extern from *: + """ + #include + + static void cuopt_lp_solve_session_capsule_dtor(PyObject *cap) noexcept + { + void *p = PyCapsule_GetPointer(cap, "cuopt.lp_solve_session"); + if (p != nullptr) { + delete reinterpret_cast(p); + } + } + """ + void cuopt_lp_solve_session_capsule_dtor(object cap) noexcept + + +cdef extern from "driver_types.h": + cdef uint32_t cudaStreamNonBlocking + + class MILPTerminationStatus(IntEnum): NoTermination = mip_termination_status_t.NoTermination Optimal = mip_termination_status_t.Optimal @@ -243,6 +270,7 @@ cdef create_solution(unique_ptr[solver_ret_t] sol_ret_ptr, cdef linear_programming_ret_t* lp_ptr cdef lp_gpu_solutions_t* gpu_sols cdef lp_cpu_solutions_t* cpu_sols + cdef object lp_solve_session_capsule = None if sol_ret.problem_type == ProblemCategory.MIP or sol_ret.problem_type == ProblemCategory.IP: # noqa mip_ptr = &sol_ret.mip_ret @@ -276,6 +304,13 @@ cdef create_solution(unique_ptr[solver_ret_t] sol_ret_ptr, else: lp_ptr = &sol_ret.lp_ret + if lp_ptr.lp_solve_session.get() != NULL: + lp_solve_session_capsule = PyCapsule_New( + lp_ptr.lp_solve_session.release(), + b"cuopt.lp_solve_session", + cuopt_lp_solve_session_capsule_dtor, + ) + # Extract solution vectors — branch only for the buffer type if lp_ptr.is_gpu(): gpu_sols = &get_gpu_lp_solutions(lp_ptr[0]) @@ -393,6 +428,7 @@ cdef create_solution(unique_ptr[solver_ret_t] sol_ret_ptr, lp_ptr.gap_, lp_ptr.nb_iterations_, lp_ptr.solved_by_, + lp_solve_session=lp_solve_session_capsule, ) else: return Solution( @@ -415,9 +451,19 @@ cdef create_solution(unique_ptr[solver_ret_t] sol_ret_ptr, ) -def Solve(py_data_model_obj, SolverSettings settings, mip=False): +def Solve(py_data_model_obj, SolverSettings settings, mip=False, session=None): cdef DataModel data_model_obj = py_data_model_obj + cdef lp_solve_session_t* session_in = NULL + + if session is not None: + if not PyCapsule_IsValid(session, b"cuopt.lp_solve_session"): + raise ValueError( + "Expected PyCapsule named 'cuopt.lp_solve_session' for session." + ) + session_in = PyCapsule_GetPointer( + session, b"cuopt.lp_solve_session" + ) data_model_obj.variable_types = type_cast( data_model_obj.variable_types, "S1", "variable_types" @@ -433,6 +479,9 @@ def Solve(py_data_model_obj, SolverSettings settings, mip=False): sol_ret_ptr = move(call_solve( data_model_obj.c_data_model_view.get(), settings.c_solver_settings.get(), + cudaStreamNonBlocking, + False, + session_in, )) return create_solution(move(sol_ret_ptr), data_model_obj) diff --git a/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pxd b/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pxd index 83a87575fe..3bd58fd605 100644 --- a/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pxd +++ b/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pxd @@ -30,6 +30,14 @@ cdef extern from "cuopt/linear_programming/pdlp/solver_settings.hpp" namespace " Barrier "cuopt::linear_programming::method_t::Barrier" # noqa Unset "cuopt::linear_programming::method_t::Unset" # noqa + cdef cppclass pdlp_solver_settings_t[i_t, f_t]: + bool session_enabled + lp_solve_session_t* lp_solve_session + +cdef extern from "cuopt/linear_programming/utilities/lp_solve_session.hpp" namespace "cuopt::cython": # noqa + cdef cppclass lp_solve_session_t: + pass + cdef extern from "cuopt/linear_programming/solver_settings.hpp" namespace "cuopt::linear_programming": # noqa cdef cppclass solver_settings_t[i_t, f_t]: @@ -90,9 +98,12 @@ cdef extern from "cuopt/linear_programming/solver_settings.hpp" namespace "cuopt void load_parameters_from_file(const string& path) except + + pdlp_solver_settings_t[i_t, f_t]& get_pdlp_settings() + cdef class SolverSettings: cdef unique_ptr[solver_settings_t[int, double]] c_solver_settings cdef public dict settings_dict cdef public object pdlp_warm_start_data cdef public list mip_callbacks + cdef public bint session_enabled diff --git a/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx b/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx index a5dcc78d18..be00999d5d 100644 --- a/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx +++ b/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx @@ -117,6 +117,7 @@ cdef class SolverSettings: self.settings_dict = {} self.pdlp_warm_start_data = None self.mip_callbacks = [] + self.session_enabled = False def to_base_type(self, value): """Convert a string to a base type. @@ -459,6 +460,16 @@ cdef class SolverSettings: warm_start_data.iterations_since_last_restart # noqa ) + c_solver_settings.get_pdlp_settings().session_enabled = self.session_enabled + + def set_session_enabled(self, enabled): + """Enable GPU barrier session persistence for repeated solves with the same sparsity.""" + self.session_enabled = True if enabled else False + + def get_session_enabled(self): + """Return whether GPU barrier session persistence is enabled.""" + return self.session_enabled + def dump_parameters_to_file(self, path, hyperparameters_only=True): """Apply ``settings_dict`` / warm start to C++, then dump parameters to *path*. diff --git a/python/cuopt/cuopt/tests/linear_programming/session_cache_helpers.py b/python/cuopt/cuopt/tests/linear_programming/session_cache_helpers.py new file mode 100644 index 0000000000..4302b4151d --- /dev/null +++ b/python/cuopt/cuopt/tests/linear_programming/session_cache_helpers.py @@ -0,0 +1,350 @@ +# SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 + +"""Shared helpers for solver-session / symbolic-cache pytest and benchmarks.""" + +from __future__ import annotations + +import io +import os +import re +import sys +from contextlib import contextmanager + +import numpy as np + +from cuopt.linear_programming import solver_settings +from cuopt.linear_programming.problem import ( + LinearExpression, + MINIMIZE, + Problem, + QuadraticExpression, +) +from cuopt.linear_programming.solver.solver_parameters import ( + CUOPT_AUGMENTED, + CUOPT_METHOD, + CUOPT_PRESOLVE, +) +from cuopt.linear_programming.solver_settings import SolverMethod + +_REUSE_SYMBOLIC_LINE = re.compile( + r"Barrier: reusing cuDSS symbolic analysis \(sparsity hash match\)" +) +_REBUILT_SYMBOLIC_LINE = re.compile( + r"Barrier: rebuilt cuDSS symbolic analysis" +) +_STORE_AUGMENTED_LINE = re.compile( + r"Barrier: stored augmented symbolic cache hash=0x[0-9a-f]+" +) +_STORE_ADAT_LINE = re.compile( + r"Barrier: stored ADAT symbolic cache hash=0x[0-9a-f]+" +) +_STORE_HASH_LINE = re.compile( + r"Barrier: stored (?:ADAT|augmented) symbolic cache hash=(0x[0-9a-f]+)" +) +_CLEAR_CACHE_LINE = re.compile( + r"Barrier: hash match but numeric refresh failed \(CUDA\); clearing symbolic cache" +) +_LINEAR_SYSTEM_ADAT = re.compile(r"Linear system\s+:\s+ADAT") +_LINEAR_SYSTEM_AUGMENTED = re.compile(r"Linear system\s+:\s+augmented") +_CACHE_PROFILE_LINE = re.compile( + r"^Cache profile: (C\d+) .+? ([0-9]+(?:\.[0-9]+)?)\s*$" +) + + +def count_log_matches(text: str, pattern: re.Pattern[str]) -> int: + return len(pattern.findall(text)) + + +def parse_cache_profile(text: str) -> dict[str, float]: + """Parse the last ``=== Solver cache profile ===`` block from cuOpt logs.""" + profiles: list[dict[str, float]] = [] + current: dict[str, float] | None = None + for line in text.splitlines(): + if "=== Solver cache profile" in line: + current = {} + continue + if line.strip() == "=== End solver cache profile ===": + if current is not None: + profiles.append(current) + current = None + continue + if current is None: + continue + m = _CACHE_PROFILE_LINE.match(line.strip()) + if m: + current[m.group(1)] = float(m.group(2)) + return profiles[-1] if profiles else {} + + +@contextmanager +def capture_solver_output(): + """Capture solver stdout and stderr (C++ logs use both).""" + read_out, write_out = os.pipe() + read_err, write_err = os.pipe() + saved_out = os.dup(1) + saved_err = os.dup(2) + capture = io.StringIO() + try: + os.dup2(write_out, 1) + os.dup2(write_err, 2) + os.close(write_out) + os.close(write_err) + yield capture + finally: + os.dup2(saved_out, 1) + os.dup2(saved_err, 2) + os.close(saved_out) + os.close(saved_err) + with os.fdopen(read_out, "r", encoding="utf-8", errors="replace") as reader_out: + out_text = reader_out.read() + with os.fdopen(read_err, "r", encoding="utf-8", errors="replace") as reader_err: + err_text = reader_err.read() + text = out_text + err_text + capture.write(text) + capture.seek(0) + if text: + sys.stdout.write(out_text) + sys.stdout.flush() + sys.stderr.write(err_text) + sys.stderr.flush() + + +def session_barrier_settings( + *, + session_enabled: bool = True, + augmented: int = -1, +) -> solver_settings.SolverSettings: + """Barrier settings with optional session and augmented/ADAT override.""" + ss = solver_settings.SolverSettings() + ss.set_parameter(CUOPT_METHOD, SolverMethod.Barrier) + ss.set_parameter(CUOPT_PRESOLVE, 0) + if augmented != -1: + ss.set_parameter(CUOPT_AUGMENTED, augmented) + if session_enabled: + ss.set_session_enabled(True) + return ss + + +def lp_problem_dims() -> tuple[int, int, int]: + """Return (n_vars, n_rows, nnz_per_row) sized for meaningful symbolic work.""" + small = os.environ.get("CUOPT_SESSION_TEST_SMALL", "").lower() in ("1", "true", "yes") + if small: + return 300, 150, 6 + return 1200, 600, 8 + + +def build_sparse_lp( + *, + seed: int = 42, + n: int | None = None, + m: int | None = None, + nnz_per_row: int | None = None, +) -> tuple[Problem, list, np.ndarray]: + """ + Random sparse LP (no quadratic term) suitable for the ADAT barrier path. + + Returns ``(problem, variables, objective_coefficients)``. + """ + default_n, default_m, default_nnz = lp_problem_dims() + n = n if n is not None else default_n + m = m if m is not None else default_m + nnz_per_row = nnz_per_row if nnz_per_row is not None else default_nnz + + rng = np.random.default_rng(seed) + prob = Problem("sparse_lp_session") + xs = [prob.addVariable(lb=0.0, name=f"x{i}") for i in range(n)] + k = min(nnz_per_row, n) + for j in range(m): + cols = rng.choice(n, size=k, replace=False) + coeffs = rng.uniform(0.5, 2.0, size=k) + rhs = float(rng.uniform(50.0, 200.0)) + prob.addConstraint( + LinearExpression([xs[i] for i in cols], coeffs.tolist(), 0.0) <= rhs + ) + c = rng.uniform(1.0, 10.0, size=n) + prob.setObjective(LinearExpression(xs, c.tolist(), 0.0), sense=MINIMIZE) + return prob, xs, c + + +def perturb_lp_values(prob: Problem, xs: list, c: np.ndarray, seed: int) -> None: + """Value-only update: objective coefficients and one constraint RHS.""" + rng = np.random.default_rng(seed) + c[:] = c * (1.0 + 0.002 * rng.standard_normal(c.size)) + lin = LinearExpression(xs, c.tolist(), 0.0) + prob.setObjective(lin, sense=MINIMIZE) + if prob.constrs: + prob.constrs[0].RHS = float(prob.constrs[0].RHS) * ( + 1.0 + 0.001 * rng.standard_normal() + ) + + +def rewire_lp_row_sparsity( + prob: Problem, + xs: list, + *, + row_idx: int = 0, + seed: int, + nnz_per_row: int | None = None, +) -> None: + """ + Change which variables appear in one constraint (same m, n; new A pattern). + + Does not add or remove variables or constraints — only the row's column indices. + """ + rng = np.random.default_rng(seed) + n = len(xs) + constr = prob.constrs[row_idx] + old_cols = set(constr.vindex_coeff_dict.keys()) + k = nnz_per_row if nnz_per_row is not None else max(len(old_cols), 1) + k = min(k, n) + for _ in range(20): + new_cols = rng.choice(n, size=k, replace=False) + if set(new_cols) != old_cols: + break + else: + raise RuntimeError("could not pick a different sparsity pattern for constraint row") + + constr.vindex_coeff_dict.clear() + for j in new_cols: + constr.vindex_coeff_dict[int(j)] = float(rng.uniform(0.5, 2.0)) + constr.RHS = float(rng.uniform(50.0, 200.0)) + prob.solved = False + prob.warmstart_data = None + prob.model = None + prob.constraint_csr_matrix = None + prob._invalidate_index_to_var_cache() + + +def build_augmented_qp( + *, + seed: int = 11, + n: int = 60, + m: int = 30, + nnz_per_row: int = 4, +) -> tuple[Problem, list, np.ndarray]: + """Small QP with off-diagonal ``Q`` (augmented KKT path).""" + rng = np.random.default_rng(seed) + prob = Problem("augmented_qp_session") + xs = [prob.addVariable(lb=0.0, ub=5.0, name=f"x{i}") for i in range(n)] + k = min(nnz_per_row, n) + for i in range(m): + cols = rng.choice(n, size=k, replace=False) + coeffs = rng.uniform(0.5, 2.0, size=k) + rhs = float(rng.uniform(10.0, 50.0)) + prob.addConstraint( + LinearExpression([xs[j] for j in cols], coeffs.tolist(), 0.0) <= rhs + ) + c = rng.uniform(-5.0, -1.0, size=n) + qv1, qv2, qc = [], [], [] + for i in range(n): + qv1.append(xs[i]) + qv2.append(xs[i]) + qc.append(float(rng.uniform(0.5, 2.0))) + for i in range(min(4, n - 1)): + qv1.append(xs[i]) + qv2.append(xs[i + 1]) + qc.append(float(rng.uniform(0.1, 0.4))) + prob.setObjective( + QuadraticExpression( + qvars1=qv1, + qvars2=qv2, + qcoefficients=qc, + vars=[], + coefficients=[], + constant=0.0, + ) + + LinearExpression(xs, c.tolist(), 0.0), + sense=MINIMIZE, + ) + return prob, xs, c + + +def perturb_qp_values(prob: Problem, xs: list, c: np.ndarray, seed: int) -> None: + rng = np.random.default_rng(seed) + c[:] = c * (1.0 + 0.002 * rng.standard_normal(c.size)) + for i, var in enumerate(xs): + var.setObjectiveCoefficient(float(c[i])) + if prob.constrs: + prob.constrs[0].RHS = float(prob.constrs[0].RHS) * ( + 1.0 + 0.001 * rng.standard_normal() + ) + prob.solved = False + prob.warmstart_data = None + + +def solve_with_log(prob: Problem, settings, session=None): + """Run ``prob.solve`` and return ``(solution, log_text, cache_profile)``.""" + with capture_solver_output() as capture: + solution = prob.solve(settings, session=session) + log_text = capture.getvalue() + return solution, log_text, parse_cache_profile(log_text) + + +def assert_optimal(solution) -> None: + assert solution.get_termination_reason() == "Optimal" + + +def assert_warm_symbolic_reuse( + cold_log: str, + warm_log: str, + cold_profile: dict[str, float], + warm_profile: dict[str, float], + *, + expect_adat: bool = False, + expect_augmented: bool = False, +) -> None: + """Cold run stores cache; warm run reuses symbolic factorization.""" + if expect_adat: + assert count_log_matches(cold_log, _LINEAR_SYSTEM_ADAT) >= 1 + assert count_log_matches(cold_log, _STORE_ADAT_LINE) >= 1 + if expect_augmented: + assert count_log_matches(cold_log, _LINEAR_SYSTEM_AUGMENTED) >= 1 + assert count_log_matches(cold_log, _STORE_AUGMENTED_LINE) >= 1 + + assert count_log_matches(cold_log, _REUSE_SYMBOLIC_LINE) == 0 + + warm_reuse = count_log_matches(warm_log, _REUSE_SYMBOLIC_LINE) + c07_c = cold_profile.get("C07", 0.0) + c07_w = warm_profile.get("C07", 0.0) + + if warm_reuse >= 1: + assert count_log_matches(warm_log, _REBUILT_SYMBOLIC_LINE) == 0 + else: + # Logs may be buffered when stderr is piped; fall back to C07 timing. + assert c07_c > 0.0, "expected non-zero C07 on cold symbolic factorization" + assert c07_w <= max(1.0, 0.05 * c07_c), ( + f"expected warm C07 near zero with reuse (cold={c07_c:.2f} ms, warm={c07_w:.2f} ms)" + ) + + +def stored_sparsity_hashes(text: str) -> list[str]: + return _STORE_HASH_LINE.findall(text) + + +def assert_full_symbolic_reanalyze( + log_text: str, + profile: dict[str, float], + *, + cold_c07: float | None = None, + allow_cache_clear: bool = False, +) -> None: + """ + Warm solve after sparsity change must not reuse cached symbolic analysis. + + When ``allow_cache_clear`` is true, a false hash match that clears the cache + (numeric refresh CUDA failure) is accepted; caller should retry the solve. + """ + if count_log_matches(log_text, _REUSE_SYMBOLIC_LINE) > 0: + raise AssertionError("unexpected symbolic reuse log after sparsity change") + if allow_cache_clear and count_log_matches(log_text, _CLEAR_CACHE_LINE) >= 1: + return + c07 = profile.get("C07", 0.0) + rebuilt = count_log_matches(log_text, _REBUILT_SYMBOLIC_LINE) >= 1 + if cold_c07 is not None and cold_c07 > 5.0: + assert rebuilt or c07 >= 0.1 * cold_c07, ( + f"expected full re-analyze (rebuilt log or C07; got C07={c07:.2f} ms, " + f"cold ref={cold_c07:.2f} ms, rebuilt={rebuilt})" + ) + else: + assert c07 > 0.0 or rebuilt diff --git a/python/cuopt/cuopt/tests/linear_programming/test_lp_solve_session.py b/python/cuopt/cuopt/tests/linear_programming/test_lp_solve_session.py new file mode 100644 index 0000000000..160c808450 --- /dev/null +++ b/python/cuopt/cuopt/tests/linear_programming/test_lp_solve_session.py @@ -0,0 +1,221 @@ +# SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 + +""" +Solver-session cache tests: ADAT symbolic reuse and sparsity-hash mismatch paths. + +Requires a CUDA GPU and ``CUOPT_CACHE_PROFILE=1`` for C07 timing assertions +(log-based checks are used when reuse lines are captured). +""" + +from __future__ import annotations + +import importlib.util +import os +from pathlib import Path + +import numpy as np + +from cuopt.linear_programming.problem import LinearExpression, MINIMIZE + +_helpers_path = Path(__file__).resolve().with_name("session_cache_helpers.py") +_spec = importlib.util.spec_from_file_location("session_cache_helpers", _helpers_path) +assert _spec and _spec.loader +_helpers = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(_helpers) + +assert_full_symbolic_reanalyze = _helpers.assert_full_symbolic_reanalyze +assert_optimal = _helpers.assert_optimal +assert_warm_symbolic_reuse = _helpers.assert_warm_symbolic_reuse +build_augmented_qp = _helpers.build_augmented_qp +build_sparse_lp = _helpers.build_sparse_lp +count_log_matches = _helpers.count_log_matches +perturb_lp_values = _helpers.perturb_lp_values +perturb_qp_values = _helpers.perturb_qp_values +rewire_lp_row_sparsity = _helpers.rewire_lp_row_sparsity +session_barrier_settings = _helpers.session_barrier_settings +solve_with_log = _helpers.solve_with_log +stored_sparsity_hashes = _helpers.stored_sparsity_hashes +_CLEAR_CACHE_LINE = _helpers._CLEAR_CACHE_LINE +_REUSE_SYMBOLIC_LINE = _helpers._REUSE_SYMBOLIC_LINE + +os.environ.setdefault("CUOPT_CACHE_PROFILE", "1") + + +def _resolve_after_cache_clear(prob, settings, session): + """If adopt hits a stale hash and clears cache, one retry must succeed.""" + sol, log_text, profile = solve_with_log(prob, settings, session=session) + if count_log_matches(log_text, _CLEAR_CACHE_LINE) >= 1: + sol, log_text, profile = solve_with_log(prob, settings, session=session) + return sol, log_text, profile + + +def test_adat_session_warm_reuse(): + """Sparse LP on ADAT path: store on cold, reuse on value-only warm solve.""" + prob, xs, c = build_sparse_lp(seed=7) + settings = session_barrier_settings(session_enabled=True, augmented=0) + + sol_cold, cold_log, cold_profile = solve_with_log(prob, settings) + assert_optimal(sol_cold) + session = sol_cold.lp_solve_session + assert session is not None + + perturb_lp_values(prob, xs, c, seed=101) + sol_warm, warm_log, warm_profile = solve_with_log(prob, settings, session=session) + assert_optimal(sol_warm) + + assert_warm_symbolic_reuse( + cold_log, warm_log, cold_profile, warm_profile, expect_adat=True + ) + + +def test_augmented_session_warm_reuse(): + """QP with off-diagonal Q (augmented KKT): store on cold, reuse on warm solve.""" + prob, xs, c = build_augmented_qp(seed=11) + settings = session_barrier_settings(session_enabled=True, augmented=-1) + + sol_cold, cold_log, cold_profile = solve_with_log(prob, settings) + assert_optimal(sol_cold) + session = sol_cold.lp_solve_session + + perturb_qp_values(prob, xs, c, seed=202) + sol_warm, warm_log, warm_profile = solve_with_log(prob, settings, session=session) + assert_optimal(sol_warm) + + assert_warm_symbolic_reuse( + cold_log, warm_log, cold_profile, warm_profile, expect_augmented=True + ) + + +def test_sparsity_hash_mismatch_add_constraint(): + """Adding a constraint changes sparsity; solver must not reuse symbolic cache.""" + prob, xs, c = build_sparse_lp(seed=13) + settings = session_barrier_settings(session_enabled=True, augmented=0) + + sol_cold, cold_log, cold_profile = solve_with_log(prob, settings) + assert_optimal(sol_cold) + session = sol_cold.lp_solve_session + cold_c07 = cold_profile.get("C07", 0.0) + cold_hashes = stored_sparsity_hashes(cold_log) + + perturb_lp_values(prob, xs, c, seed=303) + prob.addConstraint( + LinearExpression([xs[0], xs[1]], [1.0, 1.0], 0.0) <= 5.0, name="extra_cap" + ) + sol_warm, warm_log, warm_profile = _resolve_after_cache_clear( + prob, settings, session + ) + assert_optimal(sol_warm) + assert count_log_matches(warm_log, _REUSE_SYMBOLIC_LINE) == 0 + warm_hash = stored_sparsity_hashes(warm_log) + if warm_hash and cold_hashes: + assert warm_hash[-1] != cold_hashes[-1] + assert_full_symbolic_reanalyze(warm_log, warm_profile, cold_c07=cold_c07) + + +def test_sparsity_hash_mismatch_rewire_row_pattern(): + """Same m/n but a different A row pattern must not reuse symbolic cache.""" + prob, xs, c = build_sparse_lp(seed=37) + settings = session_barrier_settings(session_enabled=True, augmented=0) + + sol_cold, cold_log, cold_profile = solve_with_log(prob, settings) + assert_optimal(sol_cold) + n_before = len(prob.vars) + m_before = len(prob.constrs) + session = sol_cold.lp_solve_session + cold_c07 = cold_profile.get("C07", 0.0) + cold_hashes = stored_sparsity_hashes(cold_log) + + rewire_lp_row_sparsity(prob, xs, row_idx=0, seed=707) + assert len(prob.vars) == n_before + assert len(prob.constrs) == m_before + + sol_warm, warm_log, warm_profile = _resolve_after_cache_clear( + prob, settings, session + ) + assert_optimal(sol_warm) + assert count_log_matches(warm_log, _REUSE_SYMBOLIC_LINE) == 0 + warm_hashes = stored_sparsity_hashes(warm_log) + if cold_hashes and warm_hashes: + assert warm_hashes[-1] != cold_hashes[-1] + assert_full_symbolic_reanalyze(warm_log, warm_profile, cold_c07=cold_c07) + + +def test_sparsity_hash_mismatch_add_variable(): + """Adding a variable and constraint changes the KKT pattern.""" + prob, xs, c = build_sparse_lp(seed=17) + settings = session_barrier_settings(session_enabled=True, augmented=0) + + sol_cold, cold_log, cold_profile = solve_with_log(prob, settings) + assert_optimal(sol_cold) + session = sol_cold.lp_solve_session + cold_c07 = cold_profile.get("C07", 0.0) + cold_hashes = stored_sparsity_hashes(cold_log) + + z = prob.addVariable(lb=0.0, name="z_extra") + prob.addConstraint(LinearExpression([xs[0], z], [1.0, 1.0], 0.0) <= 10.0) + c = np.append(c, 3.0) + prob.setObjective( + LinearExpression(xs + [z], c.tolist(), 0.0), sense=MINIMIZE + ) + + sol_warm, warm_log, warm_profile = _resolve_after_cache_clear( + prob, settings, session + ) + assert_optimal(sol_warm) + assert count_log_matches(warm_log, _REUSE_SYMBOLIC_LINE) == 0 + warm_hashes = stored_sparsity_hashes(warm_log) + if cold_hashes and warm_hashes: + assert warm_hashes[-1] != cold_hashes[-1] + assert_full_symbolic_reanalyze(warm_log, warm_profile, cold_c07=cold_c07) + + +def test_cross_system_mismatch_augmented_then_adat(): + """Augmented cache in session must not be reused for a different ADAT problem.""" + prob_qp, _, _ = build_augmented_qp(seed=19) + settings_aug = session_barrier_settings(session_enabled=True, augmented=-1) + sol_qp, _, _ = solve_with_log(prob_qp, settings_aug) + assert_optimal(sol_qp) + session = sol_qp.lp_solve_session + + prob_lp, xs, c = build_sparse_lp(seed=23) + settings_adat = session_barrier_settings(session_enabled=True, augmented=0) + sol_lp, lp_log, lp_profile = _resolve_after_cache_clear( + prob_lp, settings_adat, session + ) + assert_optimal(sol_lp) + assert count_log_matches(lp_log, _REUSE_SYMBOLIC_LINE) == 0 + assert_full_symbolic_reanalyze(lp_log, lp_profile) + + perturb_lp_values(prob_lp, xs, c, seed=404) + _, warm_log, warm_profile = solve_with_log(prob_lp, settings_adat, session=session) + assert_warm_symbolic_reuse( + lp_log, warm_log, lp_profile, warm_profile, expect_adat=True + ) + + +def test_different_lp_structures_same_session(): + """Two LP instances with different sparsity patterns must not share symbolic cache.""" + prob_a, xs_a, c_a = build_sparse_lp(seed=29, n=400, m=200, nnz_per_row=5) + prob_b, xs_b, c_b = build_sparse_lp(seed=31, n=900, m=450, nnz_per_row=9) + settings = session_barrier_settings(session_enabled=True, augmented=0) + + sol_a, cold_log, profile_a = solve_with_log(prob_a, settings) + assert_optimal(sol_a) + session = sol_a.lp_solve_session + cold_c07 = profile_a.get("C07", 0.0) + cold_hashes = stored_sparsity_hashes(cold_log) + + sol_b, log_b, profile_b = _resolve_after_cache_clear(prob_b, settings, session) + assert_optimal(sol_b) + assert count_log_matches(log_b, _REUSE_SYMBOLIC_LINE) == 0 + warm_hashes = stored_sparsity_hashes(log_b) + if cold_hashes and warm_hashes: + assert warm_hashes[-1] != cold_hashes[-1] + assert_full_symbolic_reanalyze(log_b, profile_b, cold_c07=cold_c07) + + perturb_lp_values(prob_b, xs_b, c_b, seed=505) + _, warm_log, warm_profile = solve_with_log(prob_b, settings, session=session) + assert_warm_symbolic_reuse( + log_b, warm_log, profile_b, warm_profile, expect_adat=True + ) diff --git a/script_perf_eval.py b/script_perf_eval.py new file mode 100644 index 0000000000..f717ca6221 --- /dev/null +++ b/script_perf_eval.py @@ -0,0 +1,927 @@ +#!/usr/bin/env python3 +# SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 +""" +Portfolio QP benchmark for solver-session cache evaluation. + +Same problem structure as ``benchmark_cvxpylayers`` (portfolio QP with +``gamma`` quadratic in ``x`` and ``y``, transaction-cost ``z`` linearization). + +Measures cold (first) vs warm (value-only) ``prob.solve()`` times on a fixed structure. +Parameter updates are applied outside the timer. +Per-cache breakdown (C01–C09) is measured in cuOpt when ``CUOPT_CACHE_PROFILE=1`` +(see ``docs/solver_session_cache.md``). + +Environment: + CUOPT_CACHE_PROFILE=1 -> emit per-cache timings from cuOpt (default on) + CUOPT_PORTFOLIO_BENCHMARK_SMALL=1 -> k=5, scale=2 (smoke test) + CUOPT_CACHE_BENCHMARK_WARM=5 -> number of warm iterations (default 5) + CUOPT_BENCHMARK_MODE=baseline|session|all -> single mode or two-process run (default all) + +Run with conda env ``py313`` activated. + +Examples: + python script_perf_eval.py --mode baseline # session disabled, one process + python script_perf_eval.py --mode session # session enabled, one process + python script_perf_eval.py --mode all # baseline then session (new process each) +""" + +from __future__ import annotations + +import argparse +import gc +import io +import json +import os +import re +import subprocess +import sys +import time +from contextlib import contextmanager +from pathlib import Path + +import numpy as np + +from cuopt.linear_programming import solver_settings +from cuopt.linear_programming.problem import ( + LinearExpression, + MINIMIZE, + Problem, + QuadraticExpression, +) +from cuopt.linear_programming.solver.solver_parameters import CUOPT_METHOD +from cuopt.linear_programming.solver_settings import SolverMethod + + +CACHE_ITEM_NAMES = { + "C01": "raft::handle_t + stream", + "C02": "cuBLAS / cuSparse warmup", + "C03": "Problem fingerprint", + "C04": "Augmented vs ADAT choice", + "C05": "ADAT / augmented sparsity pattern", + "C06": "cuDSS handle + config", + "C07": "cuDSS symbolic factorization", + "C08": "Dense-column / SOC layout metadata", + "C09": "Device buffer allocation sizes", +} + +_CACHE_PROFILE_LINE = re.compile( + r"^Cache profile: (C\d+) .+? ([0-9]+(?:\.[0-9]+)?)\s*$" +) +_REUSE_SYMBOLIC_LINE = re.compile( + r"Barrier: reusing cuDSS symbolic analysis \(sparsity hash match\)" +) +_REBUILT_SYMBOLIC_LINE = re.compile( + r"Barrier: rebuilt cuDSS symbolic analysis" +) +_STORE_SYMBOLIC_HASH_LINE = re.compile( + r"Barrier: stored augmented symbolic cache hash=0x[0-9a-f]+" +) + + +def _parse_cache_profile(text: str) -> dict[str, float]: + """Parse the last ``=== Solver cache profile ===`` block from cuOpt logs.""" + profiles: list[dict[str, float]] = [] + current: dict[str, float] | None = None + for line in text.splitlines(): + if "=== Solver cache profile" in line: + current = {} + continue + if line.strip() == "=== End solver cache profile ===": + if current is not None: + profiles.append(current) + current = None + continue + if current is None: + continue + m = _CACHE_PROFILE_LINE.match(line.strip()) + if m: + current[m.group(1)] = float(m.group(2)) + return profiles[-1] if profiles else {} + + +def _avg_cache_profiles(profiles: list[dict[str, float]]) -> dict[str, float]: + if not profiles: + return {} + keys = sorted({k for p in profiles for k in p}) + return {k: float(np.mean([p.get(k, 0.0) for p in profiles])) for k in keys} + + +@contextmanager +def _capture_solver_output(): + """Capture C++ solver logs written directly to stderr (fd 2).""" + read_fd, write_fd = os.pipe() + saved_stderr = os.dup(2) + capture = io.StringIO() + try: + os.dup2(write_fd, 2) + os.close(write_fd) + yield capture + finally: + os.dup2(saved_stderr, 2) + os.close(saved_stderr) + with os.fdopen(read_fd, "r", encoding="utf-8", errors="replace") as reader: + text = reader.read() + capture.write(text) + capture.seek(0) + if text: + sys.stderr.write(text) + sys.stderr.flush() + + +def _objective_snapshot( + prob: Problem, + info: dict, + solution, + mu: np.ndarray, + d: float, + x0: np.ndarray, + tc: float, + label: str, +) -> dict: + x_np = portfolio_x(prob) + eval_info = {**info, "mu": mu, "x0": x0, "d": d, "tc_rate": tc} + formula_obj = float(portfolio_objective(eval_info, x_np)) + cuopt_obj = float(solution.get_primal_objective()) + return { + "label": label, + "cuopt_primal_objective": cuopt_obj, + "formula_objective": formula_obj, + "cuopt_vs_formula_delta": abs(cuopt_obj - formula_obj), + "termination": str(solution.termination_status), + "x_norm": float(np.linalg.norm(x_np)), + } + + +def timed_solve( + prob: Problem, + settings, + info: dict, + mu: np.ndarray, + d: float, + x0: np.ndarray, + tc: float, + label: str, +) -> tuple[float, dict[str, float], str, dict]: + """Run ``prob.solve``; return elapsed ms, cache profile, log text, objective snapshot.""" + t0 = time.perf_counter() + with _capture_solver_output() as capture: + solution = prob.solve(settings) + elapsed_ms = (time.perf_counter() - t0) * 1000.0 + log_text = capture.getvalue() + profile = _parse_cache_profile(log_text) + snapshot = _objective_snapshot(prob, info, solution, mu, d, x0, tc, label) + return elapsed_ms, profile, log_text, snapshot + + +def _count_log_matches(text: str, pattern: re.Pattern[str]) -> int: + return len(pattern.findall(text)) + + +def generate_portfolio_data(k=50, scale=100, seed=42): + """Identical to ``benchmark_cvxpylayers.generate_portfolio_data``.""" + np.random.seed(seed) + n = k * scale + d_diag = np.random.rand(n) * np.sqrt(k) + f_mat = np.random.randn(n, k) * (np.random.rand(n, k) < 0.5) + omega_temp = np.random.randn(k, k) + omega = (omega_temp @ omega_temp.T) / k + mu = (3.0 + 9.0 * np.random.rand(n)) / 100.0 + x0 = np.zeros(n) + gamma = 1.0 + d = 1.0 + tc_rate = 0.002 + return { + "n": n, + "k": k, + "D_diag": d_diag, + "F": f_mat, + "Omega": omega, + "mu": mu, + "x0": x0, + "gamma": gamma, + "d": d, + "tc_rate": tc_rate, + } + + +def portfolio_objective(info, x_np): + """Same evaluation as ``benchmark_cvxpylayers.portfolio_objective``.""" + y = info["F"].T @ x_np + z = np.abs(x_np - info["x0"]) + d_matrix = np.diag(info["D_diag"]) + return ( + -info["mu"] @ x_np + + info["gamma"] + * (x_np @ d_matrix @ x_np + y @ info["Omega"] @ y) + + info["tc_rate"] * np.sum(z) + ) + + +def build_cuopt_portfolio_problem(info: dict) -> Problem: + """Build portfolio QP and stash handles for value-only updates.""" + n = info["n"] + k = info["k"] + gamma = info["gamma"] + d_diag = info["D_diag"] + F = info["F"] + omega = info["Omega"] + mu = info["mu"] + x0 = info["x0"] + d_rhs = float(info["d"]) + tc_rate = float(info["tc_rate"]) + + prob = Problem("portfolio_qp_cuopt") + xs = [prob.addVariable(lb=0.0, ub=0.1, name=f"x{i}") for i in range(n)] + ys = [prob.addVariable(lb=0.0, ub=0.1, name=f"y{j}") for j in range(k)] + zs = [prob.addVariable(lb=0.0, ub=1.0e6, name=f"z{i}") for i in range(n)] + + budget_constr = prob.addConstraint(LinearExpression(xs, [1.0] * n, 0.0) == d_rhs) + + for j in range(k): + c = [-float(F[i, j]) for i in range(n)] + [1.0] + v = xs + [ys[j]] + prob.addConstraint(LinearExpression(v, c, 0.0) == 0.0) + + x0_le_constrs = [] + x0_ge_constrs = [] + for i in range(n): + x0i = float(x0[i]) + x0_le_constrs.append( + prob.addConstraint(LinearExpression([xs[i], zs[i]], [1.0, -1.0], 0.0) <= x0i) + ) + x0_ge_constrs.append( + prob.addConstraint(LinearExpression([xs[i], zs[i]], [1.0, 1.0], 0.0) >= x0i) + ) + + qv1, qv2, qc = [], [], [] + for i in range(n): + qv1.append(xs[i]) + qv2.append(xs[i]) + qc.append(gamma * float(d_diag[i])) + for p in range(k): + for q in range(k): + o = float(omega[p, q]) + if o != 0.0: + qv1.append(ys[p]) + qv2.append(ys[q]) + qc.append(gamma * o) + + quad = QuadraticExpression( + qvars1=qv1, + qvars2=qv2, + qcoefficients=qc, + vars=[], + coefficients=[], + constant=0.0, + ) + lin_vars = xs + ys + zs + lin_c = np.concatenate([-mu, np.zeros(k, dtype=np.float64), np.full(n, tc_rate)]) + linear = LinearExpression(lin_vars, lin_c.tolist(), 0.0) + prob.setObjective(quad + linear, sense=MINIMIZE) + + prob._portfolio_n = n + prob._portfolio_k = k + prob._portfolio_xs = xs + prob._portfolio_zs = zs + prob._portfolio_budget_constr = budget_constr + prob._portfolio_x0_le = x0_le_constrs + prob._portfolio_x0_ge = x0_ge_constrs + # Keep Q matrix: Problem.reset_solved_values() clears objective_qmatrix. + prob._portfolio_qmatrix = prob.objective_qmatrix + return prob + + +def _invalidate_for_resolve(prob: Problem) -> None: + """Mark problem dirty for re-solve without dropping structure or DataModel.""" + prob.solved = False + prob.warmstart_data = None + prob.objective_qmatrix = prob._portfolio_qmatrix + + +def update_portfolio_values( + prob: Problem, + mu: np.ndarray, + d: float, + x0: np.ndarray, + tc: float, +) -> None: + """Value-only update: linear objective coeffs and selected constraint RHS.""" + n = prob._portfolio_n + for i in range(n): + prob._portfolio_xs[i].setObjectiveCoefficient(-float(mu[i])) + prob._portfolio_zs[i].setObjectiveCoefficient(float(tc)) + prob._portfolio_x0_le[i].RHS = float(x0[i]) + prob._portfolio_x0_ge[i].RHS = float(x0[i]) + prob._portfolio_budget_constr.RHS = float(d) + _invalidate_for_resolve(prob) + + +def portfolio_x(prob: Problem) -> np.ndarray: + """First ``n`` primal values (portfolio weights ``x``).""" + n = prob._portfolio_n + return np.array([float(prob.vars[i].Value) for i in range(n)]) + + +def solve_portfolio( + prob: Problem, + mu: np.ndarray, + d: float, + x0: np.ndarray, + tc: float, + settings, +) -> np.ndarray: + """Update values and solve; returns primal ``x`` (first n variables).""" + update_portfolio_values(prob, mu, d, x0, tc) + prob.solve(settings) + return portfolio_x(prob) + + +def _perturb_portfolio_params(info: dict, seed: int) -> tuple[np.ndarray, float, np.ndarray, float]: + """Small value-only perturbation; structure of the Problem is unchanged.""" + rng = np.random.default_rng(seed) + mu = info["mu"] * (1.0 + 0.001 * rng.standard_normal(info["n"])) + d = float(info["d"]) * (1.0 + 0.001 * rng.standard_normal()) + x0 = info["x0"] + 1e-4 * rng.standard_normal(info["n"]) + tc = float(info["tc_rate"]) * (1.0 + 0.001 * rng.standard_normal()) + return mu, d, x0, tc + + +def _base_portfolio_params(info: dict) -> tuple[np.ndarray, float, np.ndarray, float]: + """Unperturbed parameters from ``info``.""" + return info["mu"], float(info["d"]), info["x0"], float(info["tc_rate"]) + + +def _fresh_settings(*, use_session: bool) -> solver_settings.SolverSettings: + ss = solver_settings.SolverSettings() + ss.set_parameter(CUOPT_METHOD, SolverMethod.Barrier) + if use_session: + ss.set_session_enabled(True) + return ss + + +def run_cache_benchmark( + prob: Problem, + info: dict, + settings, + n_warm: int = 5, + *, + same_values: bool = False, + use_session: bool = False, +): + """ + Cold (first) vs warm (subsequent) solve times. + + If ``same_values`` is True, every solve uses identical ``mu, d, x0, tc`` + (no parameter perturbation). Otherwise warm solves use perturbed values. + + Timings cover ``prob.solve()`` only; parameter updates run outside the timer. + """ + if same_values: + mu, d, x0, tc = _base_portfolio_params(info) + else: + mu, d, x0, tc = _perturb_portfolio_params(info, seed=0) + + if use_session: + prob._session = None + + update_portfolio_values(prob, mu, d, x0, tc) + cold_ms, cold_profile, cold_log, cold_obj = timed_solve( + prob, settings, info, mu, d, x0, tc, "cold" + ) + cold_reuse = _count_log_matches(cold_log, _REUSE_SYMBOLIC_LINE) + cold_rebuild = _count_log_matches(cold_log, _REBUILT_SYMBOLIC_LINE) + cold_store_hash = _count_log_matches(cold_log, _STORE_SYMBOLIC_HASH_LINE) + objective_records = [cold_obj] + + warm_ms = [] + warm_profiles = [] + warm_reuse_hits = [] + warm_rebuild_hits = [] + warm_store_hash_hits = [] + for i in range(n_warm): + if not same_values: + mu, d, x0, tc = _perturb_portfolio_params(info, seed=100 + i) + update_portfolio_values(prob, mu, d, x0, tc) + ms, profile, log_text, obj_snap = timed_solve( + prob, settings, info, mu, d, x0, tc, f"warm_{i}" + ) + warm_ms.append(ms) + objective_records.append(obj_snap) + warm_reuse_hits.append(_count_log_matches(log_text, _REUSE_SYMBOLIC_LINE)) + warm_rebuild_hits.append(_count_log_matches(log_text, _REBUILT_SYMBOLIC_LINE)) + warm_store_hash_hits.append(_count_log_matches(log_text, _STORE_SYMBOLIC_HASH_LINE)) + if profile: + warm_profiles.append(profile) + + warm_avg = float(np.mean(warm_ms)) + warm_min = float(np.min(warm_ms)) + saved_ms = cold_ms - warm_avg + saved_pct = (saved_ms / cold_ms * 100.0) if cold_ms > 0 else 0.0 + warm_profile_avg = _avg_cache_profiles(warm_profiles) + cache_save_ms = {} + for cid in CACHE_ITEM_NAMES: + cold_val = cold_profile.get(cid, 0.0) + warm_val = warm_profile_avg.get(cid, 0.0) + cache_save_ms[cid] = max(cold_val - warm_val, 0.0) + return { + "cold_ms": cold_ms, + "warm_ms": warm_ms, + "warm_avg_ms": warm_avg, + "warm_min_ms": warm_min, + "saved_ms": saved_ms, + "saved_pct": saved_pct, + "same_values": same_values, + "use_session": use_session, + "cold_cache_ms": cold_profile, + "warm_cache_ms_avg": warm_profile_avg, + "cache_save_ms": cache_save_ms, + "cold_reuse_log_count": cold_reuse, + "cold_rebuild_log_count": cold_rebuild, + "cold_store_hash_count": cold_store_hash, + "warm_reuse_log_counts": warm_reuse_hits, + "warm_rebuild_log_counts": warm_rebuild_hits, + "warm_store_hash_counts": warm_store_hash_hits, + "session_after_cold": prob._session is not None, + "objective_records": objective_records, + } + + +def _param_schedule(info: dict, n_warm: int, *, same_values: bool) -> list[tuple[str, np.ndarray, float, np.ndarray, float]]: + schedule: list[tuple[str, np.ndarray, float, np.ndarray, float]] = [] + if same_values: + mu, d, x0, tc = _base_portfolio_params(info) + for label in ["cold"] + [f"warm_{i}" for i in range(n_warm)]: + schedule.append((label, mu, d, x0, tc)) + else: + mu, d, x0, tc = _perturb_portfolio_params(info, seed=0) + schedule.append(("cold", mu, d, x0, tc)) + for i in range(n_warm): + mu, d, x0, tc = _perturb_portfolio_params(info, seed=100 + i) + schedule.append((f"warm_{i}", mu, d, x0, tc)) + return schedule + + +def _solve_one_objective( + info: dict, + settings, + label: str, + mu: np.ndarray, + d: float, + x0: np.ndarray, + tc: float, + *, + use_session: bool, + session=None, +) -> dict: + prob = build_cuopt_portfolio_problem(info) + if use_session and session is not None: + prob._session = session + update_portfolio_values(prob, mu, d, x0, tc) + with _capture_solver_output(): + solution = prob.solve(settings, session=session if use_session else None) + snap = _objective_snapshot(prob, info, solution, mu, d, x0, tc, label) + snap["use_session"] = use_session + snap["session_capsule"] = solution.lp_solve_session if use_session else None + return snap + + +def verify_objectives_across_modes( + info: dict, + n_warm: int, + *, + same_values: bool, + obj_atol: float = 1e-6, +) -> dict: + """ + For each parameter setting in the benchmark schedule, solve once without + session and once with session (cold then warm reuse) and compare objectives. + """ + schedule = _param_schedule(info, n_warm, same_values=same_values) + ss_base = _fresh_settings(use_session=False) + ss_sess = _fresh_settings(use_session=True) + + paired = [] + session = None + for label, mu, d, x0, tc in schedule: + base = _solve_one_objective( + info, ss_base, label, mu, d, x0, tc, use_session=False + ) + if label == "cold": + sess = _solve_one_objective( + info, ss_sess, label, mu, d, x0, tc, use_session=True + ) + session = sess.get("session_capsule") + else: + sess = _solve_one_objective( + info, + ss_sess, + label, + mu, + d, + x0, + tc, + use_session=True, + session=session, + ) + delta = abs(base["cuopt_primal_objective"] - sess["cuopt_primal_objective"]) + paired.append( + { + "label": label, + "baseline_cuopt": base["cuopt_primal_objective"], + "session_cuopt": sess["cuopt_primal_objective"], + "baseline_formula": base["formula_objective"], + "session_formula": sess["formula_objective"], + "cuopt_delta": delta, + "cuopt_match": delta <= obj_atol, + "baseline_termination": base["termination"], + "session_termination": sess["termination"], + } + ) + + all_cuopt_match = all(r["cuopt_match"] for r in paired) + return { + "paired": paired, + "all_cuopt_match": all_cuopt_match, + "obj_atol": obj_atol, + } + + +def _print_objective_records(mode: str, records: list[dict]) -> None: + print(f"\n Objectives ({mode}):") + print( + " {:>8} {:>18} {:>18} {:>12} {}".format( + "run", "cuOpt primal", "formula", "cuOpt-form", "status" + ) + ) + for rec in records: + print( + " {:>8} {:>18.10f} {:>18.10f} {:>12.2e} {}".format( + rec["label"], + rec["cuopt_primal_objective"], + rec["formula_objective"], + rec["cuopt_vs_formula_delta"], + rec["termination"], + ) + ) + + +def _print_objective_cross_check(cross: dict) -> None: + print("\n" + "=" * 50) + print("Objective cross-check (baseline vs session, same parameters)") + print("=" * 50) + print( + f" Tolerance (cuOpt primal): {cross['obj_atol']:.1e} " + f"Overall: {'PASS' if cross['all_cuopt_match'] else 'FAIL'}" + ) + print( + " {:>8} {:>18} {:>18} {:>12} {}".format( + "run", "baseline cuOpt", "session cuOpt", "delta", "match" + ) + ) + for row in cross["paired"]: + mark = "OK" if row["cuopt_match"] else "FAIL" + print( + " {:>8} {:>18.10f} {:>18.10f} {:>12.2e} {}".format( + row["label"], + row["baseline_cuopt"], + row["session_cuopt"], + row["cuopt_delta"], + mark, + ) + ) + + +def _write_results( + results_path: str, + *, + k: int, + scale: int, + n_assets: int, + n_warm: int, + build_time_s: float, + bench: dict, +) -> None: + os.makedirs(os.path.dirname(results_path), exist_ok=True) + cold_cache = bench.get("cold_cache_ms") or {} + warm_cache = bench.get("warm_cache_ms_avg") or {} + cache_save = bench.get("cache_save_ms") or {} + has_measured_cache = bool(cold_cache or warm_cache) + with open(results_path, "w", encoding="utf-8") as f: + f.write("# Cache performance results\n\n") + f.write("Generated by `script_perf_eval.py` (conda env `py313`).\n\n") + if has_measured_cache: + f.write( + "Per-cache rows are measured by cuOpt (`CUOPT_CACHE_PROFILE=1`) as " + "cold minus warm average component time.\n\n" + ) + else: + f.write( + "**Note:** Rebuild cuOpt with cache profiling enabled and set " + "`CUOPT_CACHE_PROFILE=1` for per-cache rows.\n\n" + ) + f.write(f"- Problem: portfolio QP, k={k}, scale={scale}, n={n_assets}\n") + f.write(f"- Solver: Barrier (QP)\n") + f.write(f"- Parameters: {'identical every solve' if bench.get('same_values') else 'perturbed on warm solves'}\n") + f.write(f"- Warm iterations: {n_warm}\n\n") + f.write("## Measured (aggregate)\n\n") + f.write("| Metric | ms |\n|--------|----|\n") + f.write(f"| Structure build | {build_time_s * 1000:.2f} |\n") + f.write(f"| Cold solve (1st) | {bench['cold_ms']:.2f} |\n") + f.write(f"| Warm solve (avg) | {bench['warm_avg_ms']:.2f} |\n") + f.write(f"| Warm solve (min) | {bench['warm_min_ms']:.2f} |\n") + f.write( + f"| **Aggregate save** | **{bench['saved_ms']:.2f}** " + f"({bench['saved_pct']:.1f}%) |\n\n" + ) + f.write("## Per-cache component times (ms)\n\n") + f.write( + "| ID | Cache item | Cold | Warm (avg) | Est. save (cold−warm) | " + "% of cold solve |\n" + ) + f.write("|----|------------|------|------------|----------------------|" + "-----------------|\n") + for cid, name in CACHE_ITEM_NAMES.items(): + cold_ms = cold_cache.get(cid, 0.0) + warm_ms = warm_cache.get(cid, 0.0) + save_ms = cache_save.get(cid, 0.0) + pct = (save_ms / bench["cold_ms"] * 100.0) if bench["cold_ms"] > 0 else 0.0 + f.write( + f"| {cid} | {name} | {cold_ms:.2f} | {warm_ms:.2f} | " + f"{save_ms:.2f} | {pct:.1f}% |\n" + ) + f.write("\nSee `docs/solver_session_cache.md` for cache definitions.\n") + + +def _print_bench_summary(label: str, bench: dict, n_warm: int) -> None: + print(f"\n{label}") + print("-" * 50) + print(f" session_enabled: {bench.get('use_session', False)}") + print(f" session after cold: {bench.get('session_after_cold', False)}") + print(f" Cold solve (1st): {bench['cold_ms']:.2f} ms") + print(f" Warm solve (avg): {bench['warm_avg_ms']:.2f} ms (n={n_warm})") + print(f" Warm solve (best): {bench['warm_min_ms']:.2f} ms") + print(f" Aggregate save: {bench['saved_ms']:.2f} ms ({bench['saved_pct']:.1f}%)") + c07_c = bench.get("cold_cache_ms", {}).get("C07", 0.0) + c07_w = bench.get("warm_cache_ms_avg", {}).get("C07", 0.0) + print(f" C07 symbolic (cold/warm avg): {c07_c:.2f} / {c07_w:.2f} ms") + print( + f" Sparsity hash reuse logs: cold={bench.get('cold_reuse_log_count', 0)}, " + f"warm={bench.get('warm_reuse_log_counts', [])}" + ) + print( + f" Symbolic rebuild logs: cold={bench.get('cold_rebuild_log_count', 0)}, " + f"warm={bench.get('warm_rebuild_log_counts', [])}" + ) + print( + f" Store hash logs: cold={bench.get('cold_store_hash_count', 0)}, " + f"warm={bench.get('warm_store_hash_counts', [])}" + ) + if bench.get("use_session"): + warm_reuse = bench.get("warm_reuse_log_counts", []) + if warm_reuse and all(c >= 1 for c in warm_reuse): + print(" Fingerprint/sparsity gate: PASS (warm runs reused symbolic analysis)") + elif warm_reuse and any(c >= 1 for c in warm_reuse): + print(" Fingerprint/sparsity gate: PARTIAL (some warm runs reused symbolic)") + else: + print(" Fingerprint/sparsity gate: FAIL (no warm reuse log; check session wiring)") + if bench.get("objective_records"): + _print_objective_records( + "session" if bench.get("use_session") else "baseline", + bench["objective_records"], + ) + + +def _assert_objectives(records: list[dict], *, formula_atol: float) -> None: + for rec in records: + if "Optimal" not in rec["termination"] and rec["termination"] != "1": + raise AssertionError( + f"{rec['label']}: expected Optimal termination, got {rec['termination']}" + ) + if rec["cuopt_vs_formula_delta"] > formula_atol: + raise AssertionError( + f"{rec['label']}: cuOpt vs formula delta {rec['cuopt_vs_formula_delta']:.2e} " + f"> {formula_atol:.1e}" + ) + + +def _assert_log_expectations(mode: str, bench: dict) -> None: + warm_reuse = bench.get("warm_reuse_log_counts", []) + c07_c = bench.get("cold_cache_ms", {}).get("C07", 0.0) + c07_w = bench.get("warm_cache_ms_avg", {}).get("C07", 0.0) + + if mode == "baseline": + if bench.get("cold_reuse_log_count", 0) != 0 or any(c != 0 for c in warm_reuse): + raise AssertionError( + f"baseline: unexpected sparsity reuse log " + f"(cold={bench.get('cold_reuse_log_count')}, warm={warm_reuse})" + ) + if c07_w < 50.0: + raise AssertionError( + f"baseline warm: expected C07 symbolic > 50 ms without session, got {c07_w:.2f}" + ) + print( + f" Log/fingerprint check: PASS (no reuse logs; warm C07={c07_w:.1f} ms each solve)" + ) + elif mode == "session": + if bench.get("cold_reuse_log_count", 0) != 0: + raise AssertionError("session cold: unexpected sparsity reuse log on first solve") + if c07_c < 50.0: + raise AssertionError( + f"session cold: expected C07 symbolic on first solve, got {c07_c:.2f} ms" + ) + if c07_w > 1.0: + raise AssertionError( + f"session warm: expected C07 symbolic ~0 ms with reuse, got {c07_w:.2f} ms" + ) + if any(c > 0 for c in warm_reuse): + if not all(c >= 1 for c in warm_reuse): + raise AssertionError( + f"session warm: partial reuse logs in capture {warm_reuse}" + ) + print( + f" Log/fingerprint check: PASS (reuse logs on all warm; " + f"C07 cold={c07_c:.1f} ms warm avg={c07_w:.2f} ms)" + ) + else: + print( + f" Log/fingerprint check: PASS via C07 " + f"(cold={c07_c:.1f} ms, warm avg={c07_w:.2f} ms; " + "barrier reuse lines may be stderr-buffered when piped)" + ) + else: + raise ValueError(f"unknown benchmark mode {mode!r}") + + +def _write_bench_artifact(path: Path, mode: str, bench: dict) -> None: + payload = { + "mode": mode, + "objective_records": bench["objective_records"], + "cold_ms": bench["cold_ms"], + "warm_avg_ms": bench["warm_avg_ms"], + } + path.write_text(json.dumps(payload, indent=2), encoding="utf-8") + + +def _compare_mode_objectives( + baseline_path: Path, + session_path: Path, + *, + obj_atol: float, +) -> None: + baseline = json.loads(baseline_path.read_text(encoding="utf-8")) + session = json.loads(session_path.read_text(encoding="utf-8")) + print("\n" + "=" * 50) + print("Cross-process objective check (baseline vs session)") + print("=" * 50) + all_ok = True + for b_rec, s_rec in zip(baseline["objective_records"], session["objective_records"]): + if b_rec["label"] != s_rec["label"]: + raise AssertionError("objective record label mismatch between processes") + delta = abs(b_rec["cuopt_primal_objective"] - s_rec["cuopt_primal_objective"]) + ok = delta <= obj_atol + all_ok = all_ok and ok + print( + f" {b_rec['label']:>8}: delta={delta:.2e} " + f"baseline={b_rec['cuopt_primal_objective']:.10f} " + f"session={s_rec['cuopt_primal_objective']:.10f} " + f"{'OK' if ok else 'FAIL'}" + ) + if not all_ok: + raise AssertionError("cross-process objective verification failed") + print(" Cross-process objective check: PASS") + + +def run_single_benchmark(mode: str) -> dict: + """Run cold vs warm benchmark for one mode in the current process.""" + use_session = mode == "session" + label = ( + "Session disabled (new handle each solve)" + if mode == "baseline" + else "Session enabled (reuse handle + symbolic cache)" + ) + + if os.environ.get("CUOPT_CACHE_PROFILE", "1").lower() not in ("0", "false", "no"): + os.environ["CUOPT_CACHE_PROFILE"] = "1" + + k, scale = 50, 100 + if os.environ.get("CUOPT_PORTFOLIO_BENCHMARK_SMALL", "").lower() in ("1", "true", "yes"): + k, scale = 5, 2 + + print("=" * 70) + print(f"Portfolio QP benchmark — mode={mode}") + print("=" * 70) + print(f"\nProblem settings: k={k}, scale={scale} -> n={k * scale} stocks, {k} factors") + + info = generate_portfolio_data(k=k, scale=scale, seed=42) + same_values = os.environ.get("CUOPT_CACHE_SAME_VALUES", "").lower() in ( + "1", + "true", + "yes", + ) + n_warm = int(os.environ.get("CUOPT_CACHE_BENCHMARK_WARM", "5")) + obj_atol = float(os.environ.get("CUOPT_OBJ_VERIFY_ATOL", "1e-6")) + formula_atol = float(os.environ.get("CUOPT_FORMULA_VERIFY_ATOL", "1e-4")) + + print("\nBuilding cuOpt Problem...") + t_build = time.time() + prob = build_cuopt_portfolio_problem(info) + print(f"Structure build time: {time.time() - t_build:.2f} s") + print( + "\nMode: " + + ( + "identical parameters on every solve" + if same_values + else "value-only parameter updates on warm solves" + ) + ) + + settings = _fresh_settings(use_session=use_session) + bench = run_cache_benchmark( + prob, + info, + settings, + n_warm=n_warm, + same_values=same_values, + use_session=use_session, + ) + _print_bench_summary(label, bench, n_warm) + _assert_objectives(bench["objective_records"], formula_atol=formula_atol) + print(" Objective check: PASS (Optimal; cuOpt vs formula within tolerance)") + _assert_log_expectations(mode, bench) + + artifact = Path(os.environ.get("CUOPT_BENCHMARK_ARTIFACT", f"/tmp/cuopt_bench_{mode}.json")) + _write_bench_artifact(artifact, mode, bench) + print(f" Wrote artifact: {artifact}") + return bench + + +def main(): + parser = argparse.ArgumentParser(description="Portfolio QP session cache benchmark") + parser.add_argument( + "--mode", + choices=("baseline", "session", "all"), + default=os.environ.get("CUOPT_BENCHMARK_MODE", "all"), + help="baseline=session off; session=session on (fresh process when invoked via --mode all)", + ) + args = parser.parse_args() + + try: + import torch + except ImportError: + torch = None + + if torch is not None and not torch.cuda.is_available(): + print("CUDA required.", file=sys.stderr) + sys.exit(1) + + script = Path(__file__).resolve() + env = os.environ.copy() + + if args.mode == "all": + print("=" * 70) + print("Running two-process benchmark: baseline then session") + print("=" * 70) + baseline_artifact = Path("/tmp/cuopt_bench_baseline.json") + session_artifact = Path("/tmp/cuopt_bench_session.json") + env["CUOPT_BENCHMARK_ARTIFACT"] = str(baseline_artifact) + print("\n>>> Process 1/2: baseline (session disabled)") + rc = subprocess.run( + ["stdbuf", "-eL", "-oL", sys.executable, str(script), "--mode", "baseline"], + env=env, + check=False, + ) + if rc.returncode != 0: + sys.exit(rc.returncode) + + env["CUOPT_BENCHMARK_ARTIFACT"] = str(session_artifact) + print("\n>>> Process 2/2: session (session enabled, new process)") + rc = subprocess.run( + ["stdbuf", "-eL", "-oL", sys.executable, str(script), "--mode", "session"], + env=env, + check=False, + ) + if rc.returncode != 0: + sys.exit(rc.returncode) + + obj_atol = float(os.environ.get("CUOPT_OBJ_VERIFY_ATOL", "1e-6")) + _compare_mode_objectives(baseline_artifact, session_artifact, obj_atol=obj_atol) + + baseline_bench = json.loads(baseline_artifact.read_text(encoding="utf-8")) + session_bench = json.loads(session_artifact.read_text(encoding="utf-8")) + delta = baseline_bench["warm_avg_ms"] - session_bench["warm_avg_ms"] + print("\n" + "=" * 50) + print("Timing summary (separate processes)") + print("=" * 50) + print(f" Baseline warm avg: {baseline_bench['warm_avg_ms']:.2f} ms") + print(f" Session warm avg: {session_bench['warm_avg_ms']:.2f} ms") + print(f" Session savings: {delta:.2f} ms") + print("\nAll benchmarks and assertions PASSED.") + return + + run_single_benchmark(args.mode) + print("\nAll assertions PASSED.") + + +if __name__ == "__main__": + main() diff --git a/script_session_cache_tests.py b/script_session_cache_tests.py new file mode 100644 index 0000000000..f5a31e457f --- /dev/null +++ b/script_session_cache_tests.py @@ -0,0 +1,150 @@ +#!/usr/bin/env python3 +# SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 +""" +Runnable session-cache checks: ADAT warm reuse and sparsity-hash mismatch. + + CUOPT_CACHE_PROFILE=1 python script_session_cache_tests.py + CUOPT_SESSION_TEST_SMALL=1 python script_session_cache_tests.py +""" + +from __future__ import annotations + +import os +import sys +import time + +sys.path.insert( + 0, + os.path.join( + os.path.dirname(__file__), + "python", + "cuopt", + "cuopt", + "tests", + "linear_programming", + ), +) + +from cuopt.linear_programming.problem import LinearExpression +from session_cache_helpers import ( # noqa: E402 + _CLEAR_CACHE_LINE, + _REUSE_SYMBOLIC_LINE, + _STORE_ADAT_LINE, + _STORE_AUGMENTED_LINE, + assert_full_symbolic_reanalyze, + assert_optimal, + assert_warm_symbolic_reuse, + build_augmented_qp, + build_sparse_lp, + count_log_matches, + perturb_lp_values, + perturb_qp_values, + session_barrier_settings, + solve_with_log, +) + +os.environ.setdefault("CUOPT_CACHE_PROFILE", "1") + + +def _resolve_after_cache_clear(prob, settings, session): + sol, log_text, profile = solve_with_log(prob, settings, session=session) + if count_log_matches(log_text, _CLEAR_CACHE_LINE) >= 1: + sol, log_text, profile = solve_with_log(prob, settings, session=session) + return sol, log_text, profile + + +def _run_case(name: str, fn) -> None: + print("\n" + "=" * 60) + print(name) + print("=" * 60) + t0 = time.perf_counter() + fn() + print(f"PASS ({(time.perf_counter() - t0):.1f} s)") + + +def _case_adat_warm_reuse() -> None: + prob, xs, c = build_sparse_lp(seed=7) + settings = session_barrier_settings(session_enabled=True, augmented=0) + sol_c, cold_log, cold_prof = solve_with_log(prob, settings) + assert_optimal(sol_c) + session = sol_c.lp_solve_session + perturb_lp_values(prob, xs, c, seed=101) + sol_w, warm_log, warm_prof = solve_with_log(prob, settings, session=session) + assert_optimal(sol_w) + assert_warm_symbolic_reuse(cold_log, warm_log, cold_prof, warm_prof, expect_adat=True) + print( + f" cold C07={cold_prof.get('C07', 0):.2f} ms " + f"warm C07={warm_prof.get('C07', 0):.2f} ms " + f"reuse_logs={count_log_matches(warm_log, _REUSE_SYMBOLIC_LINE)} " + f"ADAT_store={count_log_matches(cold_log, _STORE_ADAT_LINE)}" + ) + + +def _case_augmented_warm_reuse() -> None: + prob, xs, c = build_augmented_qp(seed=11) + settings = session_barrier_settings(session_enabled=True, augmented=-1) + sol_c, cold_log, cold_prof = solve_with_log(prob, settings) + assert_optimal(sol_c) + session = sol_c.lp_solve_session + perturb_qp_values(prob, xs, c, seed=202) + sol_w, warm_log, warm_prof = solve_with_log(prob, settings, session=session) + assert_optimal(sol_w) + assert_warm_symbolic_reuse( + cold_log, warm_log, cold_prof, warm_prof, expect_augmented=True + ) + print( + f" cold C07={cold_prof.get('C07', 0):.2f} ms " + f"warm C07={warm_prof.get('C07', 0):.2f} ms " + f"aug_store={count_log_matches(cold_log, _STORE_AUGMENTED_LINE)}" + ) + + +def _case_hash_mismatch_add_constraint() -> None: + prob, xs, c = build_sparse_lp(seed=13) + settings = session_barrier_settings(session_enabled=True, augmented=0) + sol_c, _, cold_prof = solve_with_log(prob, settings) + assert_optimal(sol_c) + session = sol_c.lp_solve_session + cold_c07 = cold_prof.get("C07", 0.0) + perturb_lp_values(prob, xs, c, seed=303) + prob.addConstraint(LinearExpression([xs[0], xs[1]], [1.0, 1.0], 0.0) <= 5.0) + _, warm_log, warm_prof = _resolve_after_cache_clear(prob, settings, session) + assert_full_symbolic_reanalyze(warm_log, warm_prof, cold_c07=cold_c07) + print( + f" warm C07={warm_prof.get('C07', 0):.2f} ms (cold ref={cold_c07:.2f}) " + f"reuse_logs={count_log_matches(warm_log, _REUSE_SYMBOLIC_LINE)}" + ) + + +def _case_cross_augmented_to_adat() -> None: + prob_qp, _, _ = build_augmented_qp(seed=19) + settings_aug = session_barrier_settings(session_enabled=True, augmented=-1) + sol_qp, _, _ = solve_with_log(prob_qp, settings_aug) + assert_optimal(sol_qp) + session = sol_qp.lp_solve_session + prob_lp, xs, c = build_sparse_lp(seed=23) + settings_adat = session_barrier_settings(session_enabled=True, augmented=0) + _, lp_log, lp_prof = _resolve_after_cache_clear(prob_lp, settings_adat, session) + assert_full_symbolic_reanalyze(lp_log, lp_prof) + perturb_lp_values(prob_lp, xs, c, seed=404) + _, warm_log, warm_prof = solve_with_log(prob_lp, settings_adat, session=session) + assert_warm_symbolic_reuse(lp_log, warm_log, lp_prof, warm_prof, expect_adat=True) + print(" augmented session -> ADAT LP re-analyzed, then ADAT warm reuse OK") + + +def main() -> None: + cases = [ + ("ADAT session warm reuse", _case_adat_warm_reuse), + ("Augmented QP warm reuse", _case_augmented_warm_reuse), + ("Sparsity mismatch (add constraint)", _case_hash_mismatch_add_constraint), + ("Cross-system mismatch (augmented -> ADAT)", _case_cross_augmented_to_adat), + ] + print("Session cache verification (GPU required)") + for name, fn in cases: + _run_case(name, fn) + print("\nAll session cache checks passed.") + + +if __name__ == "__main__": + main()