diff --git a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp index bd07aaac32..addd678e23 100644 --- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp @@ -22,6 +22,11 @@ #include + +namespace cuopt::cython { +class barrier_cache_t; +} + namespace cuopt { namespace CUOPT_EXPORT mathematical_optimization { @@ -365,6 +370,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 a ``barrier_cache_t`` capsule. */ + bool sequence_solve{false}; + /** Non-owning cache pointer set by ``call_solve`` for barrier symbolic reuse. */ + cuopt::cython::barrier_cache_t* barrier_cache{nullptr}; private: /** Initial primal solution */ diff --git a/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp b/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp new file mode 100644 index 0000000000..4b3643f80d --- /dev/null +++ b/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp @@ -0,0 +1,87 @@ +/* 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::mathematical_optimization::barrier { +template +class iteration_data_t; + +void destroy_iteration_data(iteration_data_t* data); + +void apply_barrier_linear_objective(iteration_data_t& data, + double const* barrier_c, + int n); +} // namespace cuopt::mathematical_optimization::barrier + +namespace cuopt { +namespace CUOPT_EXPORT cython { + +struct barrier_transform_t; + +/** + * @brief GPU solve cache owned by DataModel when sequence_solve is on. + * + * After an Optimal full solve, holds iteration_data_t and the user↔barrier transform. + * update_q crushes the new linear objective and sets c_dirty so the next Solve + * reuses that workspace (skip convert/presolve/scaling). + */ +class barrier_cache_t { + public: + static std::unique_ptr create(unsigned stream_flags); + + barrier_cache_t(barrier_cache_t&&) noexcept; + barrier_cache_t& operator=(barrier_cache_t&&) noexcept; + ~barrier_cache_t(); + + [[nodiscard]] raft::handle_t* handle_ptr(); + [[nodiscard]] raft::handle_t const* handle_ptr() const; + [[nodiscard]] rmm::cuda_stream_view stream_view() const; + + /** Drop cached iteration workspace and transform (handle/stream stay). */ + void clear(); + + /** + * @brief Take ownership of barrier iteration workspace. @p data may be null (clears). + */ + void store_iteration_data( + mathematical_optimization::barrier::iteration_data_t* data); + + /** + * @brief Release ownership of cached iteration workspace; caller must delete or wrap it. + */ + mathematical_optimization::barrier::iteration_data_t* release_iteration_data(); + + void store_transform(std::unique_ptr transform); + [[nodiscard]] barrier_transform_t* transform(); + [[nodiscard]] barrier_transform_t const* transform() const; + void set_c_dirty(bool dirty); + [[nodiscard]] bool c_dirty() const; + + /** + * Crush user-space linear objective into cached iteration_data_t.c / d_c_ and set c_dirty. + * Requires a stored transform and iteration_data from an Optimal solve. + */ + void update_linear_objective(double const* c, int n); + + private: + barrier_cache_t(std::unique_ptr stream, + std::unique_ptr handle); + + struct impl; + std::unique_ptr impl_; +}; + +} // namespace CUOPT_EXPORT cython +} // namespace cuopt diff --git a/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp b/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp index f84119a8dc..ec7736b5a5 100644 --- a/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp +++ b/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp @@ -12,7 +12,6 @@ #include #include #include - #include #include #include @@ -56,7 +55,8 @@ std::unique_ptr call_solve( cuopt::mathematical_optimization::io::data_model_view_t*, mathematical_optimization::solver_settings_t*, unsigned int flags = cudaStreamNonBlocking, - bool is_batch_mode = false); + bool is_batch_mode = false, + barrier_cache_t* cache_in = nullptr); std::pair>, double> solve_batch_remote( std::vector*>, diff --git a/cpp/include/cuopt/mathematical_optimization/utilities/cython_types.hpp b/cpp/include/cuopt/mathematical_optimization/utilities/cython_types.hpp index 69d6f91604..f047ca9637 100644 --- a/cpp/include/cuopt/mathematical_optimization/utilities/cython_types.hpp +++ b/cpp/include/cuopt/mathematical_optimization/utilities/cython_types.hpp @@ -11,6 +11,7 @@ #include #include #include +#include #include @@ -86,6 +87,9 @@ struct linear_programming_ret_t { double solve_time_{}; mathematical_optimization::method_t solved_by_{}; + /** GPU barrier cache (stream + handle + iteration workspace); moved to Python capsule when set. */ + std::unique_ptr barrier_cache; + bool is_gpu() const { return std::holds_alternative(solutions_); } }; diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index 6e1e045d19..0c34e2ca53 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -29,6 +29,12 @@ #include +#include + +#include +#include + +#include #include #include @@ -37,6 +43,7 @@ #include #include +#include #include #include #include @@ -299,8 +306,8 @@ class iteration_data_t { A(lp.A), Q(Qin), cusparse_Q_view_(lp.handle_ptr, Q), - cusparse_view_(lp.handle_ptr, lp.A), - cusparse_info(lp.handle_ptr), + cusparse_view_(lp.handle_ptr, A), + 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()), @@ -384,6 +391,9 @@ class iteration_data_t { transform_reduce_helper_(lp.handle_ptr->get_stream()), transform_reduce_pair_helper_(lp.handle_ptr->get_stream()), sum_reduce_helper_(lp.handle_ptr->get_stream()), + d_scalar_batch_(kNumScalarBatchSlots, lp.handle_ptr->get_stream()), + h_scalar_batch_(kNumScalarBatchSlots), + d_reduce_tmp_(0, lp.handle_ptr->get_stream()), indefinite_Q(false), Q_diagonal(false), symbolic_status(0), @@ -553,7 +563,7 @@ class iteration_data_t { if (n_dense_rows > 0) { settings.log.printf("Dense rows : %d\n", n_dense_rows); } - settings.log.printf("Density estimator time : %.2fs\n", column_density_time); + settings.log.printf("Density estimator time : %.3fs\n", column_density_time); if ((settings.augmented != 0) && (n_dense_columns > 50 || n_dense_rows > 10 || lp.A.m == 0 /* handle case with no constraints */ || @@ -685,32 +695,102 @@ class iteration_data_t { raft::common::nvtx::range scope("Barrier: LP Data: Cholesky init"); i_t factorization_size = use_augmented ? augmented_system_size(lp.num_cols, lp.num_rows) : lp.num_rows; - chol = std::make_unique>( + + if (use_augmented) { + raft::common::nvtx::range form_scope("Barrier: LP Data: form augmented"); + form_augmented(true); + } else { + raft::common::nvtx::range form_scope("Barrier: LP Data: form ADAT"); + form_adat(true); + } + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } + + chol = std::make_shared>( 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 + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } 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); + { + 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_augmented); + } + } + } + + device_csr_matrix_t& augmented_system() { return aug_mat(); } + const device_csr_matrix_t& augmented_system() const { return aug_mat(); } + + // Attach this solve's settings and rewind iterate-dependent state so IPM can + // Mehrotra-start with the new c. A and Q are unchanged; the previous solve + // left D and the KKT values at its last iterate. Reuse is QP-only (no cones), + // so form_*(false) updates values in the existing CSR; no symbolic rebuild. + bool prepare_for_reuse(const simplex_solver_settings_t& settings) + { + if (chol == nullptr || symbolic_status != 0) { return false; } + settings_ = settings; + if (chol != nullptr) { + static_cast*>(chol.get())->rebind_settings(settings_); + } + + { + raft::common::nvtx::range fun_scope("Barrier: reset diagonal scaling"); + const bool has_Q = Q.n > 0; + + if (has_cones()) { + primal_perturb = 1e-8; + dual_perturb = 1e-8; } else { - { - raft::common::nvtx::range form_scope("Barrier: LP Data: form ADAT"); - form_adat(true); + primal_perturb = 1e-6; + dual_perturb = 0; + } + + diag.set_scalar(1.0); + for (i_t k = 0; k < n_upper_bounds; k++) { + diag[upper_bounds[k]] = 2.0; + } + if (has_Q && !use_augmented) { + for (i_t j = 0; j < Q.n; j++) { + diag[j] += Qdiag[j]; } - if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } - symbolic_status = chol->analyze(device_ADAT); } + + inv_diag.set_scalar(1.0); + if (n_upper_bounds > 0 || (has_Q && !use_augmented)) { diag.inverse(inv_diag); } + raft::copy(d_diag_.data(), diag.data(), diag.size(), stream_view_); + 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) { + form_augmented(false); + } else { + form_adat(false); + handle_ptr->sync_stream(); + if (chol != nullptr) { chol->rebind_csr_matrix(adat_mat()); } } + if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { return false; } + + 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(); } @@ -852,7 +932,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_), @@ -866,7 +946,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; @@ -914,9 +994,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()); } { @@ -948,10 +1028,10 @@ 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_); @@ -960,8 +1040,11 @@ class iteration_data_t { if (first_call) { 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; @@ -971,21 +1054,21 @@ 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) { - settings_.log.printf("ADAT time : %.2fs\n", adat_time); + settings_.log.printf("ADAT time : %.3fs\n", adat_time); settings_.log.printf("ADAT nonzeros : %.2e\n", static_cast(adat_nnz)); 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))); } } @@ -1410,7 +1493,7 @@ class iteration_data_t { solution.z = z_tilde; dense_vector_t dual_res = z_tilde; - dual_res.axpy(-1.0, lp.objective, 1.0); + dual_res.axpy(-1.0, c, 1.0); cusparse_view.transpose_spmv(1.0, solution.y, 1.0, dual_res); if (Q.n > 0) { matrix_vector_multiply(Q, -1.0, x, 1.0, dual_res); } f_t dual_residual_norm = vector_norm_inf(dual_res, stream_view_); @@ -1995,12 +2078,27 @@ class iteration_data_t { dense_matrix_t AD_dense; dense_matrix_t H; dense_matrix_t Hchol; - const csc_matrix_t& A; + csc_matrix_t A; - const csc_matrix_t& Q; + csc_matrix_t Q; std::vector Qdiag; bool Q_diagonal; rmm::device_uvector d_augmented_diagonal_indices_; + + device_csr_matrix_t& aug_mat() { return device_augmented; } + const device_csr_matrix_t& aug_mat() const { return device_augmented; } + device_csr_matrix_t& adat_mat() { return device_ADAT; } + const device_csr_matrix_t& adat_mat() const { return device_ADAT; } + device_csr_matrix_t& a_mat() { return device_A; } + device_csc_matrix_t& ad_mat() { return device_AD; } + rmm::device_uvector& original_a_values() { return d_original_A_values; } + rmm::device_uvector& a_x_values() { return device_A_x_values; } + cusparse_info_t& spgemm_info() + { + cuopt_assert(cusparse_info_ != nullptr, "spgemm_info: cusparse workspace unset"); + return *cusparse_info_; + } + cone_kkt_data_t cone_kkt_data_; bool indefinite_Q; cusparse_view_t cusparse_Q_view_; @@ -2018,13 +2116,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_; pdlp::cusparse_dn_vec_descr_wrapper_t cusparse_tmp4_; pdlp::cusparse_dn_vec_descr_wrapper_t cusparse_h_; @@ -2111,12 +2209,20 @@ class iteration_data_t { transform_reduce_pair_helper_t transform_reduce_pair_helper_; sum_reduce_helper_t sum_reduce_helper_; + // Staging area for compute_residual_norms_mu_and_objective: several independent GPU + // reductions/dot-products write into slots of d_scalar_batch_, then a single copy into + // h_scalar_batch_ + one stream sync reads them all back at once instead of one sync each. + static constexpr i_t kNumScalarBatchSlots = 12; + rmm::device_uvector d_scalar_batch_; + pinned_dense_vector_t h_scalar_batch_; + rmm::device_buffer d_reduce_tmp_; + bool cone_combined_step_; f_t cone_sigma_mu_; rmm::cuda_stream_view stream_view_; - const simplex_solver_settings_t& settings_; + simplex_solver_settings_t settings_; }; // Move the Cholesky debug logic to a reusable function. @@ -2276,13 +2382,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) { @@ -2293,7 +2399,7 @@ int barrier_solver_t::initial_point(iteration_data_t& data) data.has_solve_info = false; // rhs_x <- b - dense_vector_t rhs_x(lp.rhs); + dense_vector_t rhs_x(data.b); dense_vector_t Fu(lp.num_cols); data.gather_upper_bounds(lp.upper, Fu); @@ -2386,7 +2492,7 @@ int barrier_solver_t::initial_point(iteration_data_t& data) // Verify A*x = b dense_vector_t init_primal_residual(lp.num_rows); - init_primal_residual = lp.rhs; + init_primal_residual = data.b; data.cusparse_view_.spmv(1.0, data.x, -1.0, init_primal_residual); data.handle_ptr->get_stream().synchronize(); #ifdef PRINT_INFO @@ -2418,7 +2524,7 @@ int barrier_solver_t::initial_point(iteration_data_t& data) // y = 0 data.y.set_scalar(0.0); - f_t epsilon = 1.0 + vector_norm1(lp.objective); + f_t epsilon = 1.0 + vector_norm1(data.c); // A^T y + z - E^T v - Q x = c // when y = 0, z - E^T v = c + Q x @@ -2483,7 +2589,7 @@ int barrier_solver_t::initial_point(iteration_data_t& data) // First compute rhs = A*Dinv*c dense_vector_t rhs(lp.num_rows); dense_vector_t Dinvc(lp.num_cols); - data.inv_diag.pairwise_product(lp.objective, Dinvc); + data.inv_diag.pairwise_product(data.c, Dinvc); // rhs = 1.0 * A * Dinv * c data.cusparse_view_.spmv(1.0, Dinvc, 0.0, rhs); @@ -2849,7 +2955,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 @@ -2869,7 +2975,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; @@ -4057,6 +4163,155 @@ void barrier_solver_t::compute_primal_dual_objective(iteration_data_t< #endif } +// Hot-loop fusion of compute_residual_norms + compute_mu + compute_primal_dual_objective: the +// three functions above each read every reduction/dot-product result back individually via a +// blocking rmm::device_scalar::value(stream) (cudaMemcpyAsync + full stream synchronize), even +// though none of these values are used until well after all of them have been computed. This +// version issues every reduction/dot-product kernel into a slot of data.d_scalar_batch_ and +// defers the host readback to a single copy + a single sync at the end. +template +void barrier_solver_t::compute_residual_norms_mu_and_objective( + iteration_data_t& data, + f_t& primal_residual_norm, + f_t& dual_residual_norm, + f_t& complementarity_residual_norm, + f_t& mu, + f_t& primal_objective, + f_t& dual_objective) +{ + raft::common::nvtx::range fun_scope("Barrier: compute_residual_norms_mu_and_objective"); + + gpu_compute_residuals(data.d_w_, data.d_x_, data.d_y_, data.d_v_, data.d_z_, data); + + constexpr i_t kSlotPrimalResidual = 0; + constexpr i_t kSlotBoundResidual = 1; + constexpr i_t kSlotDualResidual = 2; + constexpr i_t kSlotComplXzLinear = 3; + constexpr i_t kSlotComplWv = 4; + constexpr i_t kSlotComplCone = 5; + constexpr i_t kSlotMuXzSum = 6; + constexpr i_t kSlotMuWvSum = 7; + constexpr i_t kSlotCx = 8; + constexpr i_t kSlotBy = 9; + constexpr i_t kSlotUv = 10; + constexpr i_t kSlotXQx = 11; + + f_t* d_batch = data.d_scalar_batch_.data(); + + const bool has_soc = data.has_cones(); + const i_t linear_xz_size = data.linear_xz_size(data.d_complementarity_xz_residual_.size()); + auto linear_xz_span = + raft::device_span(data.d_complementarity_xz_residual_.data(), linear_xz_size); + + // All enqueue calls below must stay on stream_view_: correctness relies on strict + // single-stream FIFO ordering, so that the single sync at the bottom is enough for every + // result to be ready on the host. + enqueue_norm_inf_into(data.d_primal_residual_.data(), + data.d_primal_residual_.size(), + d_batch + kSlotPrimalResidual, + data.d_reduce_tmp_, + stream_view_); + enqueue_norm_inf_into(data.d_bound_residual_.data(), + data.d_bound_residual_.size(), + d_batch + kSlotBoundResidual, + data.d_reduce_tmp_, + stream_view_); + enqueue_norm_inf_into(data.d_dual_residual_.data(), + data.d_dual_residual_.size(), + d_batch + kSlotDualResidual, + data.d_reduce_tmp_, + stream_view_); + enqueue_norm_inf_into(linear_xz_span.data(), + linear_xz_span.size(), + d_batch + kSlotComplXzLinear, + data.d_reduce_tmp_, + stream_view_); + enqueue_norm_inf_into(data.d_complementarity_wv_residual_.data(), + data.d_complementarity_wv_residual_.size(), + d_batch + kSlotComplWv, + data.d_reduce_tmp_, + stream_view_); + + if (has_soc) { + raft::device_span cone_dot = data.cones().scratch.template get_slot<0>(); + data.cones().segmented_sum( + data.d_complementarity_xz_residual_.data() + data.cone_start(), cone_dot, stream_view_); + enqueue_max_into( + cone_dot.data(), cone_dot.size(), d_batch + kSlotComplCone, data.d_reduce_tmp_, stream_view_); + } + + enqueue_sum_into(data.d_complementarity_xz_residual_.data(), + data.d_complementarity_xz_residual_.size(), + d_batch + kSlotMuXzSum, + data.d_reduce_tmp_, + stream_view_); + enqueue_sum_into(data.d_complementarity_wv_residual_.data(), + data.d_complementarity_wv_residual_.size(), + d_batch + kSlotMuWvSum, + data.d_reduce_tmp_, + stream_view_); + + RAFT_CUBLAS_TRY(raft::linalg::detail::cublasdot(lp.handle_ptr->get_cublas_handle(), + data.d_c_.size(), + data.d_c_.data(), + 1, + data.d_x_.data(), + 1, + d_batch + kSlotCx, + stream_view_)); + RAFT_CUBLAS_TRY(raft::linalg::detail::cublasdot(lp.handle_ptr->get_cublas_handle(), + data.d_b_.size(), + data.d_b_.data(), + 1, + data.d_y_.data(), + 1, + d_batch + kSlotBy, + stream_view_)); + RAFT_CUBLAS_TRY(raft::linalg::detail::cublasdot(lp.handle_ptr->get_cublas_handle(), + data.d_restrict_u_.size(), + data.d_restrict_u_.data(), + 1, + data.d_v_.data(), + 1, + d_batch + kSlotUv, + stream_view_)); + if (data.Q.n > 0) { + auto cusparse_d_x = data.cusparse_view_.create_vector(data.d_x_); + auto cusparse_Qx = data.cusparse_view_.create_vector(data.d_Qx_); + data.cusparse_Q_view_.spmv(1.0, cusparse_d_x, 0.0, cusparse_Qx); + RAFT_CUBLAS_TRY(raft::linalg::detail::cublasdot(lp.handle_ptr->get_cublas_handle(), + data.d_Qx_.size(), + data.d_Qx_.data(), + 1, + data.d_x_.data(), + 1, + d_batch + kSlotXQx, + stream_view_)); + } + + raft::copy(data.h_scalar_batch_.data(), + data.d_scalar_batch_.data(), + data.kNumScalarBatchSlots, + stream_view_); + stream_view_.synchronize(); + + const f_t* h = data.h_scalar_batch_.data(); + + primal_residual_norm = std::max(h[kSlotPrimalResidual], h[kSlotBoundResidual]); + dual_residual_norm = h[kSlotDualResidual]; + complementarity_residual_norm = std::max(h[kSlotComplXzLinear], h[kSlotComplWv]); + if (has_soc) { + complementarity_residual_norm = std::max(complementarity_residual_norm, h[kSlotComplCone]); + } + + const f_t mu_denom = data.complementarity_degree(data.x.size(), data.n_upper_bounds); + mu = (h[kSlotMuXzSum] + h[kSlotMuWvSum]) / mu_denom; + + const f_t quad_objective = (data.Q.n > 0) ? 0.5 * h[kSlotXQx] : f_t(0); + primal_objective = h[kSlotCx] + quad_objective; + dual_objective = h[kSlotBy] - h[kSlotUv] - quad_objective; +} + template lp_status_t barrier_solver_t::check_for_suboptimal_solution( iteration_data_t& data, @@ -4090,7 +4345,7 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( solution); settings.log.printf("\n"); settings.log.printf( - "Suboptimal solution found in %d iterations and %.2f seconds\n", iter, toc(start_time)); + "Suboptimal solution found in %d iterations and %.3f seconds\n", iter, toc(start_time)); settings.log.printf("Objective %+.8e\n", compute_user_objective(lp, primal_objective)); settings.log.printf("Primal infeasibility (abs/rel): %8.2e/%8.2e\n", primal_residual_norm, @@ -4129,7 +4384,7 @@ lp_status_t barrier_solver_t::check_for_suboptimal_solution( solution); settings.log.printf("\n"); settings.log.printf( - "Suboptimal solution found in %d iterations and %.2f seconds\n", iter, toc(start_time)); + "Suboptimal solution found in %d iterations and %.3f seconds\n", iter, toc(start_time)); settings.log.printf("Objective %+.8e\n", compute_user_objective(lp, primal_objective_save)); settings.log.printf("Primal infeasibility (abs/rel): %8.2e/%8.2e\n", data.primal_residual_norm_save, @@ -4154,56 +4409,92 @@ 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::barrier_advanced_solve( + f_t start_time, lp_solution_t& solution, cuopt::cython::barrier_cache_t* cache) { - settings.log.printf("Barrier solver started at %.2f seconds\n", toc(start_time)); + settings.log.printf("Barrier solver started at %.3f seconds\n", toc(start_time)); try { - raft::common::nvtx::range fun_scope("Barrier: solve"); + raft::common::nvtx::range fun_scope("Barrier: barrier_advanced_solve"); i_t n = lp.num_cols; i_t m = lp.num_rows; - solution.resize(m, n); settings.log.printf( "Barrier solver: %d constraints, %d variables, %ld nonzeros\n", m, n, lp.A.col_start[n]); - - settings.log.printf("\n"); - if (lp.Q.n > 0) { settings.log.printf("Quadratic objective matrix : %d nonzeros\n", lp.Q.row_start[lp.Q.n]); } - if (lp.second_order_cone_dims.size() > 0) { - settings.log.printf("Second-order cones : %d\n", - static_cast(lp.second_order_cone_dims.size())); - } - // Compute the number of free variables - i_t num_free_variables = presolve_info.free_variable_pairs.size() / 2; - if (num_free_variables > 0) { - settings.log.printf("Free variables : %d\n", num_free_variables); + std::unique_ptr> owned_data; + if (cache != nullptr) { + if (auto* cached = cache->release_iteration_data()) { owned_data.reset(cached); } } - - // Compute the number of upper bounds - i_t num_upper_bounds = 0; - for (i_t j = 0; j < n; j++) { - if (lp.upper[j] < inf) { num_upper_bounds++; } + if (!owned_data) { + if (cache != nullptr) { cache->clear(); } + settings.log.printf( + "Barrier: cache reuse failed; cached iteration_data is missing or invalid\n"); + return lp_status_t::NUMERICAL_ISSUES; + } + try { + if (!owned_data->prepare_for_reuse(settings)) { + owned_data.reset(); + if (cache != nullptr) { cache->clear(); } + settings.log.printf( + "Barrier: cache reuse failed; cached iteration_data is missing or invalid\n"); + return lp_status_t::NUMERICAL_ISSUES; + } + } catch (const raft::cuda_error&) { + owned_data.reset(); + if (cache != nullptr) { cache->clear(); } + settings.log.printf( + "Barrier: cache reuse failed; cached iteration_data is missing or invalid\n"); + return lp_status_t::NUMERICAL_ISSUES; } - csc_matrix_t Q(lp.num_cols, 0, 0); - if (lp.Q.n > 0) { create_Q(lp, Q); } - - iteration_data_t data( - lp, num_upper_bounds, presolve_info.direct_free_variables, Q, settings); + iteration_data_t& data = *owned_data; + auto fail_reuse = [&](lp_status_t status) { + if (cache != nullptr) { cache->clear(); } + return status; + }; if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Barrier solver halted\n"); - return lp_status_t::CONCURRENT_LIMIT; + return fail_reuse(lp_status_t::CONCURRENT_LIMIT); } - if (data.indefinite_Q) { return lp_status_t::NUMERICAL_ISSUES; } + if (data.indefinite_Q) { return fail_reuse(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 fail_reuse(lp_status_t::NUMERICAL_ISSUES); } + settings.log.printf("Barrier setup complete at %.3f seconds\n", toc(start_time)); + return run_ipm(start_time, solution, cache, owned_data); + } catch (const raft::cuda_error& e) { + settings.log.printf("Error in barrier_solver_t: %s\n", e.what()); + return lp_status_t::NUMERICAL_ISSUES; + } catch (const std::bad_alloc& e) { + settings.log.printf("Out of memory in barrier_solver_t: %s\n", e.what()); + return lp_status_t::NUMERICAL_ISSUES; + } +} +template +lp_status_t barrier_solver_t::run_ipm( + f_t start_time, + lp_solution_t& solution, + cuopt::cython::barrier_cache_t* cache, + std::unique_ptr>& owned_data) +{ + auto finish_cache = [&](lp_status_t status) -> lp_status_t { + if (cache != nullptr && owned_data) { + if (status == lp_status_t::OPTIMAL) { + cache->store_iteration_data(owned_data.release()); + } else { + cache->clear(); + } + } + return status; + }; + iteration_data_t& data = *owned_data; + { data.cusparse_dual_residual_ = data.cusparse_view_.create_vector(data.d_dual_residual_); data.cusparse_r1_ = data.cusparse_view_.create_vector(data.d_r1_); data.cusparse_tmp4_ = data.cusparse_view_.create_vector(data.d_tmp4_); @@ -4211,118 +4502,310 @@ 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 finish_cache(lp_status_t::TIME_LIMIT); + } + + // Handle automatic adaptive regularization (-1: auto, 0: off, 1: on). + // Policy is already applied to data.dual_perturb during construction + // (before form_augmented / initial_point). + const bool adaptive_regularization = + should_use_adaptive_regularization(settings, data.has_cones()); + if (settings.barrier_adaptive_regularization == -1 && adaptive_regularization) { + settings.log.printf("Adaptive regularization enabled\n"); + } + + i_t initial_status = initial_point(data); + if (toc(start_time) > settings.time_limit) { + settings.log.printf("Barrier time limit exceeded\n"); + return finish_cache(lp_status_t::TIME_LIMIT); + } + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + settings.log.printf("Barrier solver halted\n"); + return finish_cache(lp_status_t::CONCURRENT_LIMIT); + } + if (initial_status != 0) { + settings.log.printf("Unable to compute initial point\n"); + return finish_cache(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_); + data.d_complementarity_wv_rhs_.resize(data.n_upper_bounds, stream_view_); + data.d_x_.resize(data.x.size(), stream_view_); + raft::copy(data.d_x_.data(), data.x.data(), data.x.size(), stream_view_); + data.d_y_.resize(data.y.size(), stream_view_); + raft::copy(data.d_y_.data(), data.y.data(), data.y.size(), stream_view_); + data.d_z_.resize(data.z.size(), stream_view_); + raft::copy(data.d_z_.data(), data.z.data(), data.z.size(), stream_view_); + data.d_w_.resize(data.w.size(), stream_view_); + raft::copy(data.d_w_.data(), data.w.data(), data.w.size(), stream_view_); + data.d_v_.resize(data.v.size(), stream_view_); + raft::copy(data.d_v_.data(), data.v.data(), data.v.size(), stream_view_); + data.d_upper_bounds_.resize(data.upper_bounds.size(), stream_view_); + raft::copy( + data.d_upper_bounds_.data(), data.upper_bounds.data(), data.upper_bounds.size(), stream_view_); + data.d_upper_.resize(lp.upper.size(), stream_view_); + raft::copy(data.d_upper_.data(), lp.upper.data(), lp.upper.size(), stream_view_); + data.d_bound_residual_.resize(data.n_upper_bounds, stream_view_); + + f_t primal_residual_norm, dual_residual_norm, complementarity_residual_norm; + gpu_compute_residual_norms(data.d_w_, + data.d_x_, + data.d_y_, + data.d_v_, + data.d_z_, + data, + primal_residual_norm, + dual_residual_norm, + complementarity_residual_norm); + f_t mu; + compute_mu(data, mu); + + f_t norm_b = vector_norm_inf(data.b, stream_view_); + f_t norm_c = vector_norm_inf(data.c, stream_view_); + + f_t quad_objective = 0.0; + if (data.Q.n > 0) { + dense_vector_t Qx(data.Q.n); + matrix_vector_multiply(data.Q, 1.0, data.x, 0.0, Qx); + quad_objective = 0.5 * data.x.inner_product(Qx); + } + f_t primal_objective = data.c.inner_product(data.x) + quad_objective; + + f_t relative_primal_residual = primal_residual_norm / (1.0 + norm_b); + f_t relative_dual_residual = dual_residual_norm / (1.0 + norm_c); + f_t relative_complementarity_residual = + complementarity_residual_norm / + (1.0 + + std::min(std::abs(compute_user_objective(lp, primal_objective)), std::abs(primal_objective))); + + dense_vector_t upper(lp.upper); + data.gather_upper_bounds(upper, data.restrict_u_); + data.d_restrict_u_.resize(data.restrict_u_.size(), stream_view_); + raft::copy( + data.d_restrict_u_.data(), data.restrict_u_.data(), data.restrict_u_.size(), stream_view_); + f_t dual_objective = + data.b.inner_product(data.y) - data.restrict_u_.inner_product(data.v) - quad_objective; + + f_t objective_gap_abs = std::abs(primal_objective - dual_objective); + f_t objective_gap_rel = + objective_gap_abs / + std::max(f_t(1), std::min(std::abs(primal_objective), std::abs(dual_objective))); + + data.w_save = data.w; + data.x_save = data.x; + data.y_save = data.y; + data.v_save = data.v; + data.z_save = data.z; + + i_t iter = 0; + settings.log.printf("\n"); + settings.log.printf( + " Objective Infeasibility Time\n"); + settings.log.printf( + "Iter Primal Dual Primal Dual Compl. Elapsed\n"); + float64_t elapsed_time = toc(start_time); + settings.log.printf("%3d %+.12e %+.12e %.2e %.2e %.2e %.3f\n", + iter, + compute_user_objective(lp, primal_objective), + compute_user_objective(lp, dual_objective), + relative_primal_residual, + relative_dual_residual, + relative_complementarity_residual, + elapsed_time); + + bool converged = primal_residual_norm < settings.barrier_relative_feasibility_tol && + dual_residual_norm < settings.barrier_relative_optimality_tol && + complementarity_residual_norm < settings.barrier_relative_complementarity_tol; + + const i_t iteration_limit = settings.iteration_limit; + + // Adaptive regularization for the augmented system. + f_t dual_perturb = (settings.barrier_dual_regularization >= 0) + ? settings.barrier_dual_regularization + : (adaptive_regularization ? 1e-8 : 0); + f_t primal_perturb = (settings.barrier_primal_regularization >= 0) + ? settings.barrier_primal_regularization + : (data.has_cones() ? 1e-8 : 1e-6); + + while (iter < iteration_limit) { + raft::common::nvtx::range fun_scope("Barrier: iteration"); if (toc(start_time) > settings.time_limit) { settings.log.printf("Barrier time limit exceeded\n"); - return lp_status_t::TIME_LIMIT; + return finish_cache(lp_status_t::TIME_LIMIT); } + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + settings.log.printf("Barrier solver halted\n"); + return finish_cache(lp_status_t::CONCURRENT_LIMIT); + } + + // Compute the affine step. This is the call that (re)factorizes the + // augmented system, so the IR residual here drives the adaptation of + // dual_perturb / primal_perturb for the next iteration's matrix. + compute_affine_rhs(data); + f_t max_affine_residual = 0.0; - // Handle automatic adaptive regularization (-1: auto, 0: off, 1: on). - // Policy is already applied to data.dual_perturb during construction - // (before form_augmented / initial_point). - const bool adaptive_regularization = - should_use_adaptive_regularization(settings, data.has_cones()); - if (settings.barrier_adaptive_regularization == -1 && adaptive_regularization) { - settings.log.printf("Adaptive regularization enabled\n"); + i_t status; + { + raft::common::nvtx::range fun_scope("Barrier: search_direction (affine)"); + status = + gpu_compute_search_direction(data, dual_perturb, primal_perturb, max_affine_residual); + } + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + settings.log.printf("Barrier solver halted\n"); + return finish_cache(lp_status_t::CONCURRENT_LIMIT); } - i_t initial_status = initial_point(data); + if (status < 0) { + return finish_cache(check_for_suboptimal_solution(data, + start_time, + iter, + primal_objective, + primal_residual_norm, + dual_residual_norm, + complementarity_residual_norm, + relative_primal_residual, + relative_dual_residual, + relative_complementarity_residual, + solution)); + } if (toc(start_time) > settings.time_limit) { settings.log.printf("Barrier time limit exceeded\n"); - return lp_status_t::TIME_LIMIT; + return finish_cache(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; - } - if (initial_status != 0) { - settings.log.printf("Unable to compute initial point\n"); - return lp_status_t::NUMERICAL_ISSUES; + return finish_cache(lp_status_t::CONCURRENT_LIMIT); } - // Upload initial point to device and compute initial residuals/norms on GPU - data.d_complementarity_wv_residual_.resize(data.n_upper_bounds, stream_view_); - data.d_complementarity_wv_rhs_.resize(data.n_upper_bounds, stream_view_); - data.d_x_.resize(data.x.size(), stream_view_); - raft::copy(data.d_x_.data(), data.x.data(), data.x.size(), stream_view_); - data.d_y_.resize(data.y.size(), stream_view_); - raft::copy(data.d_y_.data(), data.y.data(), data.y.size(), stream_view_); - data.d_z_.resize(data.z.size(), stream_view_); - raft::copy(data.d_z_.data(), data.z.data(), data.z.size(), stream_view_); - data.d_w_.resize(data.w.size(), stream_view_); - raft::copy(data.d_w_.data(), data.w.data(), data.w.size(), stream_view_); - data.d_v_.resize(data.v.size(), stream_view_); - raft::copy(data.d_v_.data(), data.v.data(), data.v.size(), stream_view_); - data.d_upper_bounds_.resize(data.upper_bounds.size(), stream_view_); - raft::copy(data.d_upper_bounds_.data(), - data.upper_bounds.data(), - data.upper_bounds.size(), - stream_view_); - data.d_upper_.resize(lp.upper.size(), stream_view_); - raft::copy(data.d_upper_.data(), lp.upper.data(), lp.upper.size(), stream_view_); - data.d_bound_residual_.resize(data.n_upper_bounds, stream_view_); + f_t mu_aff, sigma, new_mu; + compute_target_mu(data, mu, mu_aff, sigma, new_mu); + + compute_cc_rhs(data, new_mu); - f_t primal_residual_norm, dual_residual_norm, complementarity_residual_norm; - gpu_compute_residual_norms(data.d_w_, - data.d_x_, - data.d_y_, - data.d_v_, - data.d_z_, - data, - primal_residual_norm, - dual_residual_norm, - complementarity_residual_norm); - f_t mu; - compute_mu(data, mu); - - f_t norm_b = vector_norm_inf(data.b, stream_view_); - f_t norm_c = vector_norm_inf(data.c, stream_view_); - - f_t quad_objective = 0.0; - if (data.Q.n > 0) { - dense_vector_t Qx(data.Q.n); - matrix_vector_multiply(data.Q, 1.0, data.x, 0.0, Qx); - quad_objective = 0.5 * data.x.inner_product(Qx); + // Corrector / centering step: reuses the factorization built by the + // affine call above, so the perturbation is fixed for this solve + f_t max_corrector_residual = 0.0; + + { + raft::common::nvtx::range fun_scope("Barrier: search_direction (corrector)"); + status = + gpu_compute_search_direction(data, dual_perturb, primal_perturb, max_corrector_residual); + } + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + settings.log.printf("Barrier solver halted\n"); + return finish_cache(lp_status_t::CONCURRENT_LIMIT); + } + if (status < 0) { + return finish_cache(check_for_suboptimal_solution(data, + start_time, + iter, + primal_objective, + primal_residual_norm, + dual_residual_norm, + complementarity_residual_norm, + relative_primal_residual, + relative_dual_residual, + relative_complementarity_residual, + solution)); + } + data.has_factorization = false; + data.has_solve_info = false; + if (toc(start_time) > settings.time_limit) { + settings.log.printf("Barrier time limit exceeded\n"); + return finish_cache(lp_status_t::TIME_LIMIT); + } + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + settings.log.printf("Barrier solver halted\n"); + return finish_cache(lp_status_t::CONCURRENT_LIMIT); } - f_t primal_objective = data.c.inner_product(data.x) + quad_objective; - f_t relative_primal_residual = primal_residual_norm / (1.0 + norm_b); - f_t relative_dual_residual = dual_residual_norm / (1.0 + norm_c); - f_t relative_complementarity_residual = + compute_final_direction(data); + f_t step_primal, step_dual; + compute_primal_dual_step_length(data, settings.barrier_step_scale, step_primal, step_dual); + + compute_next_iterate(data, settings.barrier_step_scale, step_primal, step_dual); + + compute_residual_norms_mu_and_objective(data, + primal_residual_norm, + dual_residual_norm, + complementarity_residual_norm, + mu, + primal_objective, + dual_objective); + + relative_primal_residual = primal_residual_norm / (1.0 + norm_b); + relative_dual_residual = dual_residual_norm / (1.0 + norm_c); + relative_complementarity_residual = complementarity_residual_norm / (1.0 + std::min(std::abs(compute_user_objective(lp, primal_objective)), std::abs(primal_objective))); - dense_vector_t upper(lp.upper); - data.gather_upper_bounds(upper, data.restrict_u_); - data.d_restrict_u_.resize(data.restrict_u_.size(), stream_view_); - raft::copy( - data.d_restrict_u_.data(), data.restrict_u_.data(), data.restrict_u_.size(), stream_view_); - f_t dual_objective = - data.b.inner_product(data.y) - data.restrict_u_.inner_product(data.v) - quad_objective; - - f_t objective_gap_abs = std::abs(primal_objective - dual_objective); - f_t objective_gap_rel = + objective_gap_abs = std::abs(primal_objective - dual_objective); + objective_gap_rel = objective_gap_abs / std::max(f_t(1), std::min(std::abs(primal_objective), std::abs(dual_objective))); - data.w_save = data.w; - data.x_save = data.x; - data.y_save = data.y; - data.v_save = data.v; - data.z_save = data.z; + if (relative_primal_residual < settings.barrier_relaxed_feasibility_tol && + relative_dual_residual < settings.barrier_relaxed_optimality_tol && + relative_complementarity_residual < settings.barrier_relaxed_complementarity_tol) { + if (relative_primal_residual < data.relative_primal_residual_save && + relative_dual_residual < data.relative_dual_residual_save && + relative_complementarity_residual < data.relative_complementarity_residual_save && + primal_objective == primal_objective && dual_objective == dual_objective) { + settings.log.debug( + "Saving solution at iter %d: feasibility %.2e, optimality %.2e, complementarity " + "%.2e\n", + iter, + relative_primal_residual, + relative_dual_residual, + relative_complementarity_residual); + raft::copy(data.w.data(), data.d_w_.data(), data.d_w_.size(), stream_view_); + raft::copy(data.x.data(), data.d_x_.data(), data.d_x_.size(), stream_view_); + raft::copy(data.y.data(), data.d_y_.data(), data.d_y_.size(), stream_view_); + raft::copy(data.v.data(), data.d_v_.data(), data.d_v_.size(), stream_view_); + raft::copy(data.z.data(), data.d_z_.data(), data.d_z_.size(), stream_view_); + RAFT_CUDA_TRY(cudaStreamSynchronize(stream_view_)); + data.w_save = data.w; + data.x_save = data.x; + data.y_save = data.y; + data.v_save = data.v; + data.z_save = data.z; + data.relative_primal_residual_save = relative_primal_residual; + data.relative_dual_residual_save = relative_dual_residual; + data.relative_complementarity_residual_save = relative_complementarity_residual; + data.primal_residual_norm_save = primal_residual_norm; + data.dual_residual_norm_save = dual_residual_norm; + data.complementarity_residual_norm_save = complementarity_residual_norm; + } + } + + iter++; + elapsed_time = toc(start_time); + + if (primal_objective != primal_objective || dual_objective != dual_objective) { + settings.log.printf("Numerical error in objective\n"); + return finish_cache(check_for_suboptimal_solution(data, + start_time, + iter, + primal_objective, + primal_residual_norm, + dual_residual_norm, + complementarity_residual_norm, + relative_primal_residual, + relative_dual_residual, + relative_complementarity_residual, + solution)); + } - i_t iter = 0; - settings.log.printf("\n"); - settings.log.printf( - " Objective Infeasibility Time\n"); - settings.log.printf( - "Iter Primal Dual Primal Dual Compl. Elapsed\n"); - float64_t elapsed_time = toc(start_time); - settings.log.printf("%3d %+.12e %+.12e %.2e %.2e %.2e %.1f\n", + settings.log.printf("%3d %+.12e %+.12e %.2e %.2e %.2e %.3f\n", iter, compute_user_objective(lp, primal_objective), compute_user_objective(lp, dual_objective), @@ -4331,283 +4814,151 @@ lp_status_t barrier_solver_t::solve(f_t start_time, lp_solution_t= 0) - ? settings.barrier_dual_regularization - : (adaptive_regularization ? 1e-8 : 0); - f_t primal_perturb = (settings.barrier_primal_regularization >= 0) - ? settings.barrier_primal_regularization - : (data.has_cones() ? 1e-8 : 1e-6); - - while (iter < iteration_limit) { - raft::common::nvtx::range fun_scope("Barrier: iteration"); - - if (toc(start_time) > settings.time_limit) { - settings.log.printf("Barrier time limit exceeded\n"); - return 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; - } - - // Compute the affine step. This is the call that (re)factorizes the - // augmented system, so the IR residual here drives the adaptation of - // dual_perturb / primal_perturb for the next iteration's matrix. - compute_affine_rhs(data); - f_t max_affine_residual = 0.0; - - i_t status; - { - raft::common::nvtx::range fun_scope("Barrier: search_direction (affine)"); - status = - gpu_compute_search_direction(data, dual_perturb, primal_perturb, max_affine_residual); - } - if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { - settings.log.printf("Barrier solver halted\n"); - return lp_status_t::CONCURRENT_LIMIT; - } - - if (status < 0) { - return check_for_suboptimal_solution(data, - start_time, - iter, - primal_objective, - primal_residual_norm, - dual_residual_norm, - complementarity_residual_norm, - relative_primal_residual, - relative_dual_residual, - relative_complementarity_residual, - solution); - } - if (toc(start_time) > settings.time_limit) { - settings.log.printf("Barrier time limit exceeded\n"); - return 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; - } - - f_t mu_aff, sigma, new_mu; - compute_target_mu(data, mu, mu_aff, sigma, new_mu); + bool primal_feasible = relative_primal_residual < settings.barrier_relative_feasibility_tol; + bool dual_feasible = relative_dual_residual < settings.barrier_relative_optimality_tol; + bool small_gap = + relative_complementarity_residual < settings.barrier_relative_complementarity_tol; + bool small_objective_gap = + !data.has_cones() || objective_gap_rel < settings.barrier_relaxed_complementarity_tol; - compute_cc_rhs(data, new_mu); + converged = primal_feasible && dual_feasible && small_gap && small_objective_gap; - // Corrector / centering step: reuses the factorization built by the - // affine call above, so the perturbation is fixed for this solve - f_t max_corrector_residual = 0.0; + if (converged) { + settings.log.printf("\n"); + settings.log.printf( + "Optimal solution found in %d iterations and %.3fs\n", iter, toc(start_time)); + settings.log.printf("Objective %+.8e\n", compute_user_objective(lp, primal_objective)); + settings.log.printf("Primal infeasibility (abs/rel): %8.2e/%8.2e\n", + primal_residual_norm, + relative_primal_residual); + settings.log.printf("Dual infeasibility (abs/rel): %8.2e/%8.2e\n", + dual_residual_norm, + relative_dual_residual); + settings.log.printf("Complementarity gap (abs/rel): %8.2e/%8.2e\n", + complementarity_residual_norm, + relative_complementarity_residual); + settings.log.printf("\n"); + raft::copy(data.x.data(), data.d_x_.data(), data.d_x_.size(), stream_view_); + raft::copy(data.y.data(), data.d_y_.data(), data.d_y_.size(), stream_view_); + raft::copy(data.z.data(), data.d_z_.data(), data.d_z_.size(), stream_view_); + raft::copy(data.v.data(), data.d_v_.data(), data.d_v_.size(), stream_view_); + RAFT_CUDA_TRY(cudaStreamSynchronize(stream_view_)); + data.to_solution(lp, + iter, + primal_objective, + compute_user_objective(lp, primal_objective), + primal_residual_norm, + data.cusparse_view_, + solution); + return finish_cache(lp_status_t::OPTIMAL); + } - { - raft::common::nvtx::range fun_scope("Barrier: search_direction (corrector)"); - status = - gpu_compute_search_direction(data, dual_perturb, primal_perturb, max_corrector_residual); - } - if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { - settings.log.printf("Barrier solver halted\n"); - return lp_status_t::CONCURRENT_LIMIT; + // Check if the solution is getting worse + if (data.Q.n > 0 && + ((!primal_feasible && + relative_primal_residual > 100 * data.relative_primal_residual_save) || + (!dual_feasible && relative_dual_residual > 100 * data.relative_dual_residual_save) || + (!small_gap && relative_complementarity_residual > + 10000 * data.relative_complementarity_residual_save))) { + if (data.relative_primal_residual_save < settings.barrier_relaxed_feasibility_tol && + data.relative_dual_residual_save < settings.barrier_relaxed_optimality_tol && + data.relative_complementarity_residual_save < + settings.barrier_relaxed_complementarity_tol) { + return finish_cache(check_for_suboptimal_solution(data, + start_time, + iter, + primal_objective, + primal_residual_norm, + dual_residual_norm, + complementarity_residual_norm, + relative_primal_residual, + relative_dual_residual, + relative_complementarity_residual, + solution)); } - if (status < 0) { - return check_for_suboptimal_solution(data, - start_time, - iter, - primal_objective, - primal_residual_norm, - dual_residual_norm, - complementarity_residual_norm, - relative_primal_residual, - relative_dual_residual, - relative_complementarity_residual, - solution); - } - data.has_factorization = false; - data.has_solve_info = false; - if (toc(start_time) > settings.time_limit) { - settings.log.printf("Barrier time limit exceeded\n"); - return 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; - } - - compute_final_direction(data); - f_t step_primal, step_dual; - compute_primal_dual_step_length(data, settings.barrier_step_scale, step_primal, step_dual); + } + } + raft::copy(data.x.data(), data.d_x_.data(), data.d_x_.size(), stream_view_); + raft::copy(data.y.data(), data.d_y_.data(), data.d_y_.size(), stream_view_); + raft::copy(data.z.data(), data.d_z_.data(), data.d_z_.size(), stream_view_); + raft::copy(data.v.data(), data.d_v_.data(), data.d_v_.size(), stream_view_); + RAFT_CUDA_TRY(cudaStreamSynchronize(stream_view_)); + data.to_solution(lp, + iter, + primal_objective, + compute_user_objective(lp, primal_objective), + primal_residual_norm, + data.cusparse_view_, + solution); + return finish_cache(lp_status_t::ITERATION_LIMIT); +} - compute_next_iterate(data, settings.barrier_step_scale, step_primal, step_dual); +template +lp_status_t barrier_solver_t::solve(f_t start_time, + lp_solution_t& solution, + cuopt::cython::barrier_cache_t* cache) +{ + settings.log.printf("Barrier solver started at %.3f seconds\n", toc(start_time)); + try { + raft::common::nvtx::range fun_scope("Barrier: solve"); - compute_residual_norms( - data, primal_residual_norm, dual_residual_norm, complementarity_residual_norm); + i_t n = lp.num_cols; + i_t m = lp.num_rows; - compute_mu(data, mu); + solution.resize(m, n); + settings.log.printf( + "Barrier solver: %d constraints, %d variables, %ld nonzeros\n", m, n, lp.A.col_start[n]); - compute_primal_dual_objective(data, primal_objective, dual_objective); + settings.log.printf("\n"); - relative_primal_residual = primal_residual_norm / (1.0 + norm_b); - relative_dual_residual = dual_residual_norm / (1.0 + norm_c); - relative_complementarity_residual = - complementarity_residual_norm / - (1.0 + std::min(std::abs(compute_user_objective(lp, primal_objective)), - std::abs(primal_objective))); + if (lp.Q.n > 0) { + settings.log.printf("Quadratic objective matrix : %d nonzeros\n", lp.Q.row_start[lp.Q.n]); + } + if (lp.second_order_cone_dims.size() > 0) { + settings.log.printf("Second-order cones : %d\n", + static_cast(lp.second_order_cone_dims.size())); + } - objective_gap_abs = std::abs(primal_objective - dual_objective); - objective_gap_rel = - objective_gap_abs / - std::max(f_t(1), std::min(std::abs(primal_objective), std::abs(dual_objective))); + i_t num_free_variables = presolve_info.free_variable_pairs.size() / 2; + if (num_free_variables > 0) { + settings.log.printf("Free variables : %d\n", num_free_variables); + } - if (relative_primal_residual < settings.barrier_relaxed_feasibility_tol && - relative_dual_residual < settings.barrier_relaxed_optimality_tol && - relative_complementarity_residual < settings.barrier_relaxed_complementarity_tol) { - if (relative_primal_residual < data.relative_primal_residual_save && - relative_dual_residual < data.relative_dual_residual_save && - relative_complementarity_residual < data.relative_complementarity_residual_save && - primal_objective == primal_objective && dual_objective == dual_objective) { - settings.log.debug( - "Saving solution at iter %d: feasibility %.2e, optimality %.2e, complementarity " - "%.2e\n", - iter, - relative_primal_residual, - relative_dual_residual, - relative_complementarity_residual); - raft::copy(data.w.data(), data.d_w_.data(), data.d_w_.size(), stream_view_); - raft::copy(data.x.data(), data.d_x_.data(), data.d_x_.size(), stream_view_); - raft::copy(data.y.data(), data.d_y_.data(), data.d_y_.size(), stream_view_); - raft::copy(data.v.data(), data.d_v_.data(), data.d_v_.size(), stream_view_); - raft::copy(data.z.data(), data.d_z_.data(), data.d_z_.size(), stream_view_); - RAFT_CUDA_TRY(cudaStreamSynchronize(stream_view_)); - data.w_save = data.w; - data.x_save = data.x; - data.y_save = data.y; - data.v_save = data.v; - data.z_save = data.z; - data.relative_primal_residual_save = relative_primal_residual; - data.relative_dual_residual_save = relative_dual_residual; - data.relative_complementarity_residual_save = relative_complementarity_residual; - data.primal_residual_norm_save = primal_residual_norm; - data.dual_residual_norm_save = dual_residual_norm; - data.complementarity_residual_norm_save = complementarity_residual_norm; - } - } + i_t num_upper_bounds = 0; + for (i_t j = 0; j < n; j++) { + if (lp.upper[j] < inf) { num_upper_bounds++; } + } - iter++; - elapsed_time = toc(start_time); - - if (primal_objective != primal_objective || dual_objective != dual_objective) { - settings.log.printf("Numerical error in objective\n"); - return check_for_suboptimal_solution(data, - start_time, - iter, - primal_objective, - primal_residual_norm, - dual_residual_norm, - complementarity_residual_norm, - relative_primal_residual, - relative_dual_residual, - relative_complementarity_residual, - solution); - } + csc_matrix_t Q(lp.num_cols, 0, 0); + std::unique_ptr> owned_data; - settings.log.printf("%3d %+.12e %+.12e %.2e %.2e %.2e %.1f\n", - iter, - compute_user_objective(lp, primal_objective), - compute_user_objective(lp, dual_objective), - relative_primal_residual, - relative_dual_residual, - relative_complementarity_residual, - elapsed_time); - - bool primal_feasible = relative_primal_residual < settings.barrier_relative_feasibility_tol; - bool dual_feasible = relative_dual_residual < settings.barrier_relative_optimality_tol; - bool small_gap = - relative_complementarity_residual < settings.barrier_relative_complementarity_tol; - bool small_objective_gap = - !data.has_cones() || objective_gap_rel < settings.barrier_relaxed_complementarity_tol; - - converged = primal_feasible && dual_feasible && small_gap && small_objective_gap; - - if (converged) { - settings.log.printf("\n"); - settings.log.printf( - "Optimal solution found in %d iterations and %.3fs\n", iter, toc(start_time)); - settings.log.printf("Objective %+.8e\n", compute_user_objective(lp, primal_objective)); - settings.log.printf("Primal infeasibility (abs/rel): %8.2e/%8.2e\n", - primal_residual_norm, - relative_primal_residual); - settings.log.printf("Dual infeasibility (abs/rel): %8.2e/%8.2e\n", - dual_residual_norm, - relative_dual_residual); - settings.log.printf("Complementarity gap (abs/rel): %8.2e/%8.2e\n", - complementarity_residual_norm, - relative_complementarity_residual); - settings.log.printf("\n"); - raft::copy(data.x.data(), data.d_x_.data(), data.d_x_.size(), stream_view_); - raft::copy(data.y.data(), data.d_y_.data(), data.d_y_.size(), stream_view_); - raft::copy(data.z.data(), data.d_z_.data(), data.d_z_.size(), stream_view_); - raft::copy(data.v.data(), data.d_v_.data(), data.d_v_.size(), stream_view_); - RAFT_CUDA_TRY(cudaStreamSynchronize(stream_view_)); - data.to_solution(lp, - iter, - primal_objective, - compute_user_objective(lp, primal_objective), - primal_residual_norm, - data.cusparse_view_, - solution); - return lp_status_t::OPTIMAL; - } + if (cache != nullptr) { cache->store_iteration_data(nullptr); } + if (lp.Q.n > 0) { create_Q(lp, Q); } + owned_data = std::make_unique>( + lp, num_upper_bounds, presolve_info.direct_free_variables, Q, settings); + iteration_data_t& data = *owned_data; - // Check if the solution is getting worse - if (data.Q.n > 0 && - ((!primal_feasible && - relative_primal_residual > 100 * data.relative_primal_residual_save) || - (!dual_feasible && relative_dual_residual > 100 * data.relative_dual_residual_save) || - (!small_gap && relative_complementarity_residual > - 10000 * data.relative_complementarity_residual_save))) { - if (data.relative_primal_residual_save < settings.barrier_relaxed_feasibility_tol && - data.relative_dual_residual_save < settings.barrier_relaxed_optimality_tol && - data.relative_complementarity_residual_save < - settings.barrier_relaxed_complementarity_tol) { - return check_for_suboptimal_solution(data, - start_time, - iter, - primal_objective, - primal_residual_norm, - dual_residual_norm, - complementarity_residual_norm, - relative_primal_residual, - relative_dual_residual, - relative_complementarity_residual, - solution); - } - } + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + settings.log.printf("Barrier solver halted\n"); + if (cache != nullptr) { cache->clear(); } + return lp_status_t::CONCURRENT_LIMIT; } - raft::copy(data.x.data(), data.d_x_.data(), data.d_x_.size(), stream_view_); - raft::copy(data.y.data(), data.d_y_.data(), data.d_y_.size(), stream_view_); - raft::copy(data.z.data(), data.d_z_.data(), data.d_z_.size(), stream_view_); - raft::copy(data.v.data(), data.d_v_.data(), data.d_v_.size(), stream_view_); - RAFT_CUDA_TRY(cudaStreamSynchronize(stream_view_)); - data.to_solution(lp, - iter, - primal_objective, - compute_user_objective(lp, primal_objective), - primal_residual_norm, - data.cusparse_view_, - solution); - return lp_status_t::ITERATION_LIMIT; + if (data.indefinite_Q) { + if (cache != nullptr) { cache->clear(); } + return lp_status_t::NUMERICAL_ISSUES; + } + if (data.symbolic_status != 0) { + settings.log.printf("Error in symbolic analysis\n"); + if (cache != nullptr) { cache->clear(); } + return lp_status_t::NUMERICAL_ISSUES; + } + + settings.log.printf("Barrier setup complete at %.3f seconds\n", toc(start_time)); + return run_ipm(start_time, solution, cache, owned_data); } catch (const raft::cuda_error& e) { settings.log.printf("Error in barrier_solver_t: %s\n", e.what()); return lp_status_t::NUMERICAL_ISSUES; } catch (const std::bad_alloc& e) { - // Covers rmm::out_of_memory and any other allocation failure. The barrier sizes its normal - // equations from the problem, so a shape it cannot hold is a property of the input rather - // than a defect, and the solvers running concurrently with it are unaffected. settings.log.printf("Out of memory in barrier_solver_t: %s\n", e.what()); return lp_status_t::NUMERICAL_ISSUES; } @@ -4622,4 +4973,20 @@ template class sparse_cholesky_cudss_t; template class iteration_data_t; #endif +void destroy_iteration_data(iteration_data_t* data) { delete data; } + +void apply_barrier_linear_objective(iteration_data_t& data, + double const* barrier_c, + int n) +{ + if (barrier_c == nullptr || static_cast(data.c.size()) != n || + static_cast(data.d_c_.size()) != n) { + throw std::invalid_argument( + "update_q: barrier linear objective size does not match cached iteration_data_t."); + } + std::copy(barrier_c, barrier_c + n, data.c.data()); + raft::copy( + data.d_c_.data(), data.c.data(), static_cast(n), data.handle_ptr->get_stream()); +} + } // namespace cuopt::mathematical_optimization::barrier diff --git a/cpp/src/barrier/barrier.hpp b/cpp/src/barrier/barrier.hpp index f72b2ff728..efd7355bd0 100644 --- a/cpp/src/barrier/barrier.hpp +++ b/cpp/src/barrier/barrier.hpp @@ -17,9 +17,16 @@ #include #include +#include + #include #include + +namespace cuopt::cython { +class barrier_cache_t; +} // namespace cuopt::cython + namespace cuopt::mathematical_optimization::barrier { /** Validates SOC layout on an simplex::lp_problem_t before barrier presolve/solve. */ @@ -36,9 +43,21 @@ class barrier_solver_t { barrier_solver_t(const simplex::lp_problem_t& lp, const simplex::presolve_info_t& presolve, const simplex::simplex_solver_settings_t& settings); - simplex::lp_status_t solve(f_t start_time, simplex::lp_solution_t& solution); + simplex::lp_status_t solve(f_t start_time, + simplex::lp_solution_t& solution, + cuopt::cython::barrier_cache_t* cache = nullptr); + // Cache reuse: cached iteration_data_t already has the updated linear objective. + // Prepare the workspace, compute a new initial point, run IPM. Same status/solution contract as + // solve(). + simplex::lp_status_t barrier_advanced_solve(f_t start_time, + simplex::lp_solution_t& solution, + cuopt::cython::barrier_cache_t* cache); private: + simplex::lp_status_t run_ipm(f_t start_time, + simplex::lp_solution_t& solution, + cuopt::cython::barrier_cache_t* cache, + std::unique_ptr>& owned_data); void my_pop_range(bool debug) const; void create_Q(const simplex::lp_problem_t& lp, csc_matrix_t& Q); int initial_point(iteration_data_t& data); @@ -65,6 +84,13 @@ class barrier_solver_t { void compute_primal_dual_objective(iteration_data_t& data, f_t& primal_objective, f_t& dual_objective); + void compute_residual_norms_mu_and_objective(iteration_data_t& data, + f_t& primal_residual_norm, + f_t& dual_residual_norm, + f_t& complementarity_residual_norm, + f_t& mu, + f_t& primal_objective, + f_t& dual_objective); // To be able to directly pass lambdas to transform functions public: diff --git a/cpp/src/barrier/cusparse_view.cu b/cpp/src/barrier/cusparse_view.cu index 477200c5e9..e67eb7bfe0 100644 --- a/cpp/src/barrier/cusparse_view.cu +++ b/cpp/src/barrier/cusparse_view.cu @@ -246,6 +246,16 @@ cusparse_view_t::~cusparse_view_t() if (A_T_ != nullptr) { 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 pdlp::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 ea6bf363b9..1815b8a9b7 100644 --- a/cpp/src/barrier/cusparse_view.hpp +++ b/cpp/src/barrier/cusparse_view.hpp @@ -57,6 +57,8 @@ class cusparse_view_t { f_t beta, pdlp::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 974e2b0f4a..7272b9efa0 100644 --- a/cpp/src/barrier/device_sparse_matrix.cuh +++ b/cpp/src/barrier/device_sparse_matrix.cuh @@ -179,6 +179,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), @@ -318,6 +322,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 01045847d1..216e8a7ebe 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::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,9 +160,10 @@ 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::get_driver_entry_point("cuGetErrorString"); + // 1. Set up the GPU resources CUdevResource initial_device_GPU_resources = {}; auto cuDeviceGetDevResource_func = cuopt::get_driver_entry_point("cuDeviceGetDevResource"); @@ -277,18 +283,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( @@ -307,7 +313,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 +379,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { 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::get_driver_entry_point("cuStreamDestroy"); CU_CHECK(reinterpret_cast(cuStreamDestroy_func)(stream), reinterpret_cast(cuGetErrorString_func)); @@ -396,9 +402,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 +413,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 +488,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 : %.3fs\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,32 +517,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("Symbolic factorization time : %.3fs\n", symbolic_factorization_time); 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 @@ -545,7 +562,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; } @@ -555,11 +572,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); @@ -571,7 +588,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; } @@ -585,16 +602,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 : %.3fs\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; @@ -608,15 +625,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); @@ -690,7 +707,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(); @@ -700,7 +717,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; } @@ -714,8 +731,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 : %.3fs\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; @@ -726,7 +743,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? @@ -739,10 +756,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; } @@ -766,7 +783,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; } @@ -779,16 +796,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 : %.3fs\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; @@ -817,11 +834,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; } @@ -831,11 +848,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); @@ -851,11 +868,11 @@ 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", - vector_norm2(b_host), - compute_hash(b_host), - vector_norm2(x_host), - compute_hash(x_host)); + 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), + compute_hash(x_host)); #endif return 0; @@ -866,6 +883,63 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { this->positive_definite = positive_definite; } + void rebind_settings(const simplex::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; @@ -889,7 +963,10 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t { f_t* x_values_d; f_t* b_values_d; - const simplex::simplex_solver_settings_t& settings_; + bool symbolic_done_; + bool numeric_factor_valid_; + const simplex::simplex_solver_settings_t* settings_; + 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 388bb43b35..719d8cc45a 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -26,17 +26,75 @@ #include #include +#include +#include + #include #include #include +#include #include +#include #include namespace cuopt::mathematical_optimization::simplex { namespace { +template +void unscale_uncrush_barrier_to_user(const user_problem_t& user_problem, + const raft::handle_t* handle_ptr, + i_t original_num_rows, + i_t original_num_cols, + const lp_problem_t& barrier_lp, + const presolve_info_t& presolve_info, + const std::vector& column_scales, + const std::vector& row_scales, + const simplex_solver_settings_t& barrier_settings, + const lp_solution_t& barrier_solution, + lp_solution_t& solution) +{ + std::vector unscaled_x(barrier_lp.num_cols); + std::vector unscaled_y(barrier_lp.num_rows); + std::vector unscaled_z(barrier_lp.num_cols); + unscale_solution(column_scales, + row_scales, + barrier_solution.x, + barrier_solution.y, + barrier_solution.z, + unscaled_x, + unscaled_y, + unscaled_z); + + // Dummy converted LP: sizes only. Bound-free=0 so uncrush_solution never reads A. + lp_problem_t converted(handle_ptr, original_num_rows, original_num_cols, 0); + lp_solution_t lp_solution(original_num_rows, original_num_cols); + uncrush_solution(presolve_info, + barrier_settings, + converted, + unscaled_x, + unscaled_y, + unscaled_z, + lp_solution.x, + lp_solution.y, + lp_solution.z); + + uncrush_primal_solution(user_problem, converted, lp_solution.x, solution.x); + uncrush_dual_solution( + user_problem, converted, lp_solution.y, lp_solution.z, solution.y, solution.z); + solution.objective = + barrier_solution.user_objective / user_problem.obj_scale - user_problem.obj_constant; + solution.user_objective = barrier_solution.user_objective; + solution.l2_primal_residual = barrier_solution.l2_primal_residual; + solution.l2_dual_residual = barrier_solution.l2_dual_residual; + solution.iterations = barrier_solution.iterations; +} + +} // namespace + +namespace { + template void write_matlab(const std::string& filename, const simplex::lp_problem_t& lp) { @@ -358,16 +416,57 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us const simplex_solver_settings_t& settings, f_t start_time, lp_solution_t& solution, + cuopt::cython::barrier_cache_t* cache, const raft::handle_t* handle_ptr) { - lp_status_t status = lp_status_t::UNSET; + lp_status_t status = lp_status_t::UNSET; + simplex_solver_settings_t barrier_settings = settings; + + auto const* xf = + (cache != nullptr && cache->c_dirty()) ? cache->transform() : nullptr; + const bool reuse_c_only = + xf != nullptr && xf->barrier_lp != nullptr && !user_problem.Q_values.empty() && + user_problem.second_order_cone_dims.empty() && xf->second_order_cone_dims.empty() && + xf->barrier_lp->second_order_cone_dims.empty() && + settings.barrier_presolve_bound_free_variables == 0 && + user_problem.num_cols == xf->user_num_cols && + user_problem.num_rows == xf->user_num_rows; + + if (reuse_c_only) { + settings.log.printf("Barrier: reusing cache (skip convert/presolve/scaling)\n"); + lp_solution_t barrier_solution(xf->barrier_lp->num_rows, + xf->barrier_lp->num_cols); + barrier::barrier_solver_t barrier_solver( + *xf->barrier_lp, xf->presolve_info, barrier_settings); + lp_status_t barrier_status = + barrier_solver.barrier_advanced_solve(start_time, barrier_solution, cache); + if (barrier_status == lp_status_t::OPTIMAL) { + unscale_uncrush_barrier_to_user(user_problem, + cache->handle_ptr(), + xf->original_num_rows, + xf->original_num_cols, + *xf->barrier_lp, + xf->presolve_info, + xf->column_scales, + xf->row_scales, + barrier_settings, + barrier_solution, + solution); + cache->set_c_dirty(false); + } else { + cache->clear(); + } + return barrier_status; + } + lp_problem_t original_lp(handle_ptr, 1, 1, 1); // Convert the user problem to a linear program with only equality constraints std::vector new_slacks; - simplex_solver_settings_t barrier_settings = settings; dualize_info_t dualize_info; - convert_user_problem(user_problem, barrier_settings, original_lp, new_slacks, dualize_info); + { + convert_user_problem(user_problem, barrier_settings, original_lp, new_slacks, dualize_info); + } if (!barrier::validate_barrier_cone_layout(original_lp, barrier_settings)) { return lp_status_t::NUMERICAL_ISSUES; } @@ -377,7 +476,10 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us // Presolve the linear program presolve_info_t presolve_info; lp_problem_t presolved_lp(handle_ptr, 1, 1, 1); - const i_t ok = presolve(original_lp, barrier_settings, presolved_lp, presolve_info); + i_t ok; + { + ok = presolve(original_lp, barrier_settings, presolved_lp, presolve_info); + } if (ok == CONCURRENT_HALT_RETURN) { return lp_status_t::CONCURRENT_LIMIT; } if (ok == TIME_LIMIT_RETURN) { return lp_status_t::TIME_LIMIT; } if (ok == -1) { return lp_status_t::INFEASIBLE; } @@ -389,13 +491,56 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us presolved_lp.A.col_start[presolved_lp.num_cols]); std::vector column_scales; std::vector row_scales; - scaling(presolved_lp, barrier_settings, barrier_lp, column_scales, row_scales); + { + scaling(presolved_lp, barrier_settings, barrier_lp, column_scales, row_scales); + } // Solve using barrier lp_solution_t barrier_solution(barrier_lp.num_rows, barrier_lp.num_cols); barrier::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, cache); + + if (cache != nullptr) { + if (barrier_status == lp_status_t::OPTIMAL) { + auto xf = std::make_unique(); + xf->user_num_cols = user_problem.num_cols; + xf->user_num_rows = user_problem.num_rows; + xf->original_num_cols = original_lp.num_cols; + xf->original_num_rows = original_lp.num_rows; + xf->obj_scale = user_problem.obj_scale; + xf->obj_constant = user_problem.obj_constant; + xf->row_sense = user_problem.row_sense; + xf->cone_var_start = user_problem.cone_var_start; + xf->second_order_cone_dims = user_problem.second_order_cone_dims; + xf->expanded_original_num_cols = user_problem.original_num_cols; + xf->original_col_to_expanded_col = user_problem.original_col_to_expanded_col; + xf->presolve_info = presolve_info; + xf->column_scales = column_scales; + xf->row_scales = row_scales; + xf->barrier_lp = std::make_unique>(barrier_lp); + { + try { + auto crushed = cuopt::cython::crush_user_linear_objective( + *xf, user_problem.objective.data(), user_problem.num_cols); + xf->linear_obj_shift.resize(static_cast(barrier_lp.num_cols), 0.0); + if (static_cast(crushed.size()) == barrier_lp.num_cols) { + for (int j = 0; j < barrier_lp.num_cols; ++j) { + xf->linear_obj_shift[static_cast(j)] = + barrier_lp.objective[static_cast(j)] - + crushed[static_cast(j)]; + } + } + } catch (std::exception const&) { + xf->linear_obj_shift.assign(static_cast(barrier_lp.num_cols), 0.0); + } + } + cache->store_transform(std::move(xf)); + } else { + cache->clear(); + } + } + if (barrier_status == lp_status_t::OPTIMAL) { #ifdef COMPUTE_SCALED_RESIDUALS std::vector scaled_residual = barrier_lp.rhs; @@ -688,20 +833,23 @@ 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, - f_t start_time, - lp_solution_t& solution) + lp_solution_t& solution, + cuopt::cython::barrier_cache_t* cache) { + f_t start_time = tic(); return solve_linear_program_with_barrier( - user_problem, settings, start_time, solution, user_problem.handle_ptr); + user_problem, settings, start_time, solution, cache, user_problem.handle_ptr); } 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) + f_t start_time, + lp_solution_t& solution, + cuopt::cython::barrier_cache_t* cache) { - 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, cache, user_problem.handle_ptr); } template @@ -845,19 +993,22 @@ 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::barrier_cache_t* cache); 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::barrier_cache_t* cache); 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, + cuopt::cython::barrier_cache_t* cache, const raft::handle_t* handle_ptr); template lp_status_t solve_linear_program(const user_problem_t& user_problem, diff --git a/cpp/src/dual_simplex/solve.hpp b/cpp/src/dual_simplex/solve.hpp index 308c462de5..69046b6998 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 barrier_cache_t; +} // namespace cuopt::cython + namespace cuopt::mathematical_optimization::simplex { template @@ -91,21 +95,26 @@ 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::barrier_cache_t* cache = 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::barrier_cache_t* cache = 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, + cuopt::cython::barrier_cache_t* cache, const raft::handle_t* handle_ptr); template diff --git a/cpp/src/linear_algebra/vector_math.cuh b/cpp/src/linear_algebra/vector_math.cuh index ac9d24001b..fe3e11bd99 100644 --- a/cpp/src/linear_algebra/vector_math.cuh +++ b/cpp/src/linear_algebra/vector_math.cuh @@ -14,6 +14,7 @@ #include #include +#include #include #include #include @@ -68,6 +69,71 @@ f_t device_custom_vector_norm_inf(InputIteratorT in, i_t size, rmm::cuda_stream_ return d_out.value(stream_view); } +// Same reduction as device_custom_vector_norm_inf, but writes into a caller-supplied device +// pointer (and reuses a caller-supplied temp-storage buffer) instead of allocating a private +// rmm::device_scalar and blocking on .value(). Lets callers batch several reductions and defer +// the host readback to a single copy + sync. +template +void enqueue_norm_inf_into( + InputIteratorT in, i_t size, f_t* out, rmm::device_buffer& tmp, rmm::cuda_stream_view stream_view) +{ + if (size == 0) { + RAFT_CUDA_TRY(cudaMemsetAsync(out, 0, sizeof(f_t), stream_view.value())); + return; + } + size_t temp_storage_bytes = 0; + f_t init = 0; + auto custom_op = norm_inf_max{}; + cub::DeviceReduce::Reduce( + nullptr, temp_storage_bytes, in, out, size, custom_op, init, stream_view); + + tmp.resize(temp_storage_bytes, stream_view); + + cub::DeviceReduce::Reduce( + tmp.data(), temp_storage_bytes, in, out, size, custom_op, init, stream_view); +} + +// Sum reduction into a caller-supplied device pointer/temp-storage buffer, deferring the host +// readback (see enqueue_norm_inf_into). +template +void enqueue_sum_into( + InputIteratorT in, i_t size, f_t* out, rmm::device_buffer& tmp, rmm::cuda_stream_view stream_view) +{ + if (size == 0) { + RAFT_CUDA_TRY(cudaMemsetAsync(out, 0, sizeof(f_t), stream_view.value())); + return; + } + size_t temp_storage_bytes = 0; + cub::DeviceReduce::Sum(nullptr, temp_storage_bytes, in, out, size, stream_view); + + tmp.resize(temp_storage_bytes, stream_view); + + cub::DeviceReduce::Sum(tmp.data(), temp_storage_bytes, in, out, size, stream_view); +} + +// Max reduction (with a floor of 0, matching this codebase's existing +// thrust::reduce(..., f_t(0), thrust::maximum()) usage) into a caller-supplied device +// pointer/temp-storage buffer, deferring the host readback (see enqueue_norm_inf_into). +template +void enqueue_max_into( + InputIteratorT in, i_t size, f_t* out, rmm::device_buffer& tmp, rmm::cuda_stream_view stream_view) +{ + if (size == 0) { + RAFT_CUDA_TRY(cudaMemsetAsync(out, 0, sizeof(f_t), stream_view.value())); + return; + } + size_t temp_storage_bytes = 0; + f_t init = 0; + auto custom_op = thrust::maximum{}; + cub::DeviceReduce::Reduce( + nullptr, temp_storage_bytes, in, out, size, custom_op, init, stream_view); + + tmp.resize(temp_storage_bytes, stream_view); + + cub::DeviceReduce::Reduce( + tmp.data(), temp_storage_bytes, in, out, size, custom_op, init, stream_view); +} + template f_t device_vector_norm_inf(const rmm::device_uvector& in, rmm::cuda_stream_view stream_view) { diff --git a/cpp/src/pdlp/CMakeLists.txt b/cpp/src/pdlp/CMakeLists.txt index 2f90f94872..9673d283dc 100644 --- a/cpp/src/pdlp/CMakeLists.txt +++ b/cpp/src/pdlp/CMakeLists.txt @@ -39,6 +39,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/barrier_cache.cu ${CMAKE_CURRENT_SOURCE_DIR}/cuopt_c.cpp ) diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index d08a36d178..f5f8192f81 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -35,6 +35,8 @@ #include #include #include +#include +#include #include #include @@ -62,6 +64,7 @@ #include #include +#include #include #include #include @@ -72,6 +75,33 @@ namespace cuopt::mathematical_optimization { +namespace { + +template +simplex::user_problem_t user_problem_from_transform( + raft::handle_t const* handle_ptr, + optimization_problem_t& model, + cuopt::cython::barrier_transform_t const& xf) +{ + simplex::user_problem_t user_problem(handle_ptr); + user_problem.num_rows = xf.user_num_rows; + user_problem.num_cols = xf.user_num_cols; + user_problem.objective = model.get_objective_coefficients_host(); + user_problem.row_sense = xf.row_sense; + user_problem.rhs.assign(static_cast(xf.user_num_rows), f_t(0)); + user_problem.obj_scale = static_cast(xf.obj_scale); + user_problem.obj_constant = static_cast(xf.obj_constant); + // Nonempty Q so the cache-reuse path accepts this as a QP (it rejects empty Q). + user_problem.Q_values.assign(1, f_t(1)); + user_problem.cone_var_start = xf.cone_var_start; + user_problem.second_order_cone_dims = xf.second_order_cone_dims; + user_problem.original_num_cols = xf.expanded_original_num_cols; + user_problem.original_col_to_expanded_col = xf.original_col_to_expanded_col; + return user_problem; +} + +} // namespace + template extern rmm::device_uvector gpu_cast(const rmm::device_uvector& src, rmm::cuda_stream_view stream); @@ -492,7 +522,8 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t const simplex::user_problem_t& user_problem, pdlp_solver_settings_t const& settings, const timer_t& timer, - const raft::handle_t* handle_ptr) + const raft::handle_t* handle_ptr, + cuopt::cython::barrier_cache_t* cache = nullptr) { f_t norm_user_objective = vector_norm2(user_problem.objective); f_t norm_rhs = vector_norm2(user_problem.rhs); @@ -533,7 +564,7 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t simplex::lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); auto status = simplex::solve_linear_program_with_barrier( - user_problem, barrier_settings, timer.get_tic_start(), solution, handle_ptr); + user_problem, barrier_settings, timer.get_tic_start(), solution, cache, handle_ptr); if (status == simplex::lp_status_t::OPTIMAL) { barrier::project_barrier_solution_to_model_variables(user_problem, solution); @@ -557,12 +588,14 @@ template optimization_problem_solution_t run_barrier( mip::problem_t& problem, pdlp_solver_settings_t const& settings, - const timer_t& timer) + const timer_t& timer, + cuopt::cython::barrier_cache_t* cache = nullptr) { // Convert data structures to dual simplex format and back simplex::user_problem_t dual_simplex_problem = cuopt_problem_to_user_problem(problem.handle_ptr, problem, false); - auto sol_dual_simplex = run_barrier(dual_simplex_problem, settings, timer, problem.handle_ptr); + auto sol_dual_simplex = + run_barrier(dual_simplex_problem, settings, timer, problem.handle_ptr, cache); return convert_dual_simplex_sol(problem, std::get<0>(sol_dual_simplex), std::get<1>(sol_dual_simplex), @@ -1813,7 +1846,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.barrier_cache); } else if (settings.method == method_t::Concurrent) { return run_concurrent(problem, settings, timer, is_batch_mode); } else { @@ -1845,7 +1878,18 @@ optimization_problem_solution_t solve_qcqp( auto qcqp_timer = cuopt::timer_t(settings.time_limit); - if (problem_checking) { + auto* cache = settings.barrier_cache; + auto const* xf = (cache != nullptr && cache->c_dirty()) ? cache->transform() : nullptr; + const bool reuse_from_cache = + settings.user_problem_file.empty() && xf != nullptr && xf->barrier_lp != nullptr && + settings.barrier_presolve_bound_free_variables == 0 && + op_problem.has_quadratic_objective() && !op_problem.has_quadratic_constraints() && + xf->second_order_cone_dims.empty() && + static_cast(xf->row_sense.size()) == xf->user_num_rows && + op_problem.get_n_variables() == xf->user_num_cols && + op_problem.get_n_constraints() == xf->user_num_rows; + + if (problem_checking && !reuse_from_cache) { problem_checking_t::check_problem_representation(op_problem); if (problem_checking_t::has_crossing_bounds(op_problem)) { return optimization_problem_solution_t( @@ -1873,12 +1917,20 @@ 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); } - // Convert data structures to dual simplex format and back - 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, op_problem.get_handle_ptr()); - auto solution = convert_dual_simplex_sol(op_problem, + simplex::user_problem_t dual_simplex_problem(op_problem.get_handle_ptr()); + if (reuse_from_cache) { + dual_simplex_problem = user_problem_from_transform( + op_problem.get_handle_ptr(), op_problem, *settings.barrier_cache->transform()); + } else { + 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, + op_problem.get_handle_ptr(), + settings.barrier_cache); + auto solution = convert_dual_simplex_sol(op_problem, std::get<0>(sol_dual_simplex), std::get<1>(sol_dual_simplex), std::get<2>(sol_dual_simplex), diff --git a/cpp/src/pdlp/utilities/barrier_cache.cu b/cpp/src/pdlp/utilities/barrier_cache.cu new file mode 100644 index 0000000000..9690d8a59e --- /dev/null +++ b/cpp/src/pdlp/utilities/barrier_cache.cu @@ -0,0 +1,134 @@ +/* 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 +#include +#include + +namespace cuopt::cython { + +using barrier_iteration_data_t = + mathematical_optimization::barrier::iteration_data_t; +using barrier_iteration_data_ptr = + std::unique_ptr; + +struct barrier_cache_t::impl { + impl(std::unique_ptr stream_in, std::unique_ptr handle_in) + : stream(std::move(stream_in)), + handle(std::move(handle_in)), + iteration_data(nullptr, &mathematical_optimization::barrier::destroy_iteration_data) + { + } + + std::unique_ptr stream; + std::unique_ptr handle; + barrier_iteration_data_ptr iteration_data; + std::unique_ptr transform; + bool c_dirty{false}; +}; + +barrier_cache_t::barrier_cache_t(std::unique_ptr stream, + std::unique_ptr handle) + : impl_(std::make_unique(std::move(stream), std::move(handle))) +{ +} + +barrier_cache_t::~barrier_cache_t() = default; + +barrier_cache_t::barrier_cache_t(barrier_cache_t&&) noexcept = default; +barrier_cache_t& barrier_cache_t::operator=(barrier_cache_t&&) noexcept = default; + +std::unique_ptr barrier_cache_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 barrier_cache_t(std::move(stream), std::move(handle))); +} + +raft::handle_t* barrier_cache_t::handle_ptr() +{ + return impl_->handle.get(); +} + +raft::handle_t const* barrier_cache_t::handle_ptr() const +{ + return impl_->handle.get(); +} + +rmm::cuda_stream_view barrier_cache_t::stream_view() const +{ + return impl_->stream->view(); +} + +void barrier_cache_t::clear() +{ + impl_->iteration_data.reset(); + impl_->transform.reset(); + impl_->c_dirty = false; +} + +void barrier_cache_t::store_iteration_data(barrier_iteration_data_t* data) +{ + impl_->iteration_data.reset(data); +} + +barrier_iteration_data_t* barrier_cache_t::release_iteration_data() +{ + return impl_->iteration_data.release(); +} + +void barrier_cache_t::store_transform(std::unique_ptr transform) +{ + impl_->transform = std::move(transform); +} + +barrier_transform_t* barrier_cache_t::transform() { return impl_->transform.get(); } + +barrier_transform_t const* barrier_cache_t::transform() const { return impl_->transform.get(); } + +void barrier_cache_t::set_c_dirty(bool dirty) { impl_->c_dirty = dirty; } + +bool barrier_cache_t::c_dirty() const +{ + return impl_->c_dirty && impl_->transform != nullptr && impl_->iteration_data.get() != nullptr; +} + +void barrier_cache_t::update_linear_objective(double const* c, int n) +{ + cuopt_expects(impl_->transform != nullptr, + error_type_t::ValidationError, + "update_q: no barrier transform; Solve with sequence_solve first."); + cuopt_expects(impl_->iteration_data.get() != nullptr, + error_type_t::ValidationError, + "update_q: no cached iteration_data; Solve a QP to Optimal first."); + std::vector crushed; + try { + crushed = crush_user_linear_objective(*impl_->transform, c, n); + } catch (std::invalid_argument const& e) { + cuopt_expects(false, error_type_t::ValidationError, "%s", e.what()); + } + if (impl_->transform->linear_obj_shift.size() == crushed.size()) { + for (std::size_t j = 0; j < crushed.size(); ++j) { + crushed[j] += impl_->transform->linear_obj_shift[j]; + } + } + try { + mathematical_optimization::barrier::apply_barrier_linear_objective( + *impl_->iteration_data, crushed.data(), static_cast(crushed.size())); + } catch (std::invalid_argument const& e) { + cuopt_expects(false, error_type_t::ValidationError, "%s", e.what()); + } + impl_->c_dirty = true; +} + +} // namespace cuopt::cython diff --git a/cpp/src/pdlp/utilities/barrier_transform.hpp b/cpp/src/pdlp/utilities/barrier_transform.hpp new file mode 100644 index 0000000000..f4c6185498 --- /dev/null +++ b/cpp/src/pdlp/utilities/barrier_transform.hpp @@ -0,0 +1,105 @@ +/* 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::cython { + +/** + * User ↔ barrier transform retained on barrier_cache_t after Optimal: + * convert / presolve / scaling, plus the scaled LP. + * Enough to crush a new user-space linear objective into barrier space and to + * uncrush a solution without rerunning those algorithms. + */ +struct barrier_transform_t { + int user_num_cols{0}; + int user_num_rows{0}; + int original_num_cols{0}; + int original_num_rows{0}; + double obj_scale{1.0}; + double obj_constant{0.0}; + + // Enough of the user problem for reuse uncrush without rebuilding A. + std::vector row_sense; + int cone_var_start{0}; + std::vector second_order_cone_dims; + int expanded_original_num_cols{0}; + std::vector original_col_to_expanded_col; + + cuopt::mathematical_optimization::simplex::presolve_info_t presolve_info; + std::vector column_scales; + std::vector row_scales; + // Barrier linear objective minus crush(user c) from the first solve (Q·ℓ shift, etc.). + std::vector linear_obj_shift; + std::unique_ptr> barrier_lp; +}; + +inline std::vector crush_user_linear_objective(barrier_transform_t const& xf, + double const* c, + int n) +{ + if (c == nullptr || n != xf.user_num_cols) { + throw std::invalid_argument( + "update_q: linear objective length must match the cached user column count."); + } + if (xf.original_num_cols < xf.user_num_cols) { + throw std::invalid_argument("update_q: cached original column count is smaller than user n."); + } + if (xf.barrier_lp == nullptr) { + throw std::invalid_argument("update_q: cached barrier LP is missing."); + } + + std::vector orig(static_cast(xf.original_num_cols), 0.0); + for (int j = 0; j < n; ++j) { + orig[static_cast(j)] = c[j]; + } + for (int j : xf.presolve_info.negated_variables) { + orig[static_cast(j)] *= -1.0; + } + + std::vector presolved; + if (!xf.presolve_info.remaining_variables.empty()) { + presolved.resize(xf.presolve_info.remaining_variables.size()); + for (std::size_t k = 0; k < xf.presolve_info.remaining_variables.size(); ++k) { + presolved[k] = orig[static_cast(xf.presolve_info.remaining_variables[k])]; + } + } else { + presolved = std::move(orig); + } + + auto const& pairs = xf.presolve_info.free_variable_pairs; + if (!pairs.empty()) { + if (pairs.size() % 2 != 0) { + throw std::invalid_argument("update_q: free_variable_pairs size is not even."); + } + std::size_t extra = pairs.size() / 2; + presolved.resize(presolved.size() + extra); + for (std::size_t k = 0; k < extra; ++k) { + int u = pairs[2 * k]; + int v = pairs[2 * k + 1]; + presolved[static_cast(v)] = -presolved[static_cast(u)]; + } + } + + if (static_cast(presolved.size()) != xf.barrier_lp->num_cols || + xf.column_scales.size() != presolved.size()) { + throw std::invalid_argument( + "update_q: crushed objective size does not match barrier columns / column_scales."); + } + for (std::size_t j = 0; j < presolved.size(); ++j) { + presolved[j] /= xf.column_scales[j]; + } + return presolved; +} + +} // namespace cuopt::cython diff --git a/cpp/src/pdlp/utilities/cython_solve.cu b/cpp/src/pdlp/utilities/cython_solve.cu index ed77c4f722..78c2c5c192 100644 --- a/cpp/src/pdlp/utilities/cython_solve.cu +++ b/cpp/src/pdlp/utilities/cython_solve.cu @@ -6,6 +6,7 @@ /* clang-format on */ #include + #include #include #include @@ -18,6 +19,8 @@ #include #include #include +#include + #include #include #include @@ -30,6 +33,8 @@ #include #include +#include + #include namespace cuopt { @@ -96,24 +101,54 @@ std::unique_ptr call_solve( cuopt::mathematical_optimization::io::data_model_view_t* data_model, cuopt::mathematical_optimization::solver_settings_t* solver_settings, unsigned int flags, - bool is_batch_mode) + bool is_batch_mode, + barrier_cache_t* cache_in) { raft::common::nvtx::range fun_scope("Call Solve"); + 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::mathematical_optimization::get_memory_backend_type(); solver_ret_t response; + auto& pdlp_settings = solver_settings->get_pdlp_settings(); + const bool sequence_solve = pdlp_settings.sequence_solve; + const bool barrier_path = + data_model->has_quadratic_objective() || data_model->has_quadratic_constraints() || + pdlp_settings.method == cuopt::mathematical_optimization::method_t::Barrier; + const bool want_cache = (cache_in != nullptr || sequence_solve) && barrier_path && + memory_backend == cuopt::mathematical_optimization::memory_backend_t::GPU && + !is_batch_mode; + + std::unique_ptr owned_cache; + barrier_cache_t* active_cache = cache_in; + pdlp_settings.barrier_cache = 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::mathematical_optimization::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_cache) { + if (active_cache == nullptr) { + owned_cache = barrier_cache_t::create(flags); + active_cache = owned_cache.get(); + } + solve_handle = active_cache->handle_ptr(); + pdlp_settings.barrier_cache = active_cache; + } - auto problem = cuopt::mathematical_optimization::optimization_problem_t(&handle_); + auto problem = cuopt::mathematical_optimization::optimization_problem_t(solve_handle); cuopt::mathematical_optimization::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() == mathematical_optimization::problem_category_t::LP) { @@ -142,6 +177,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_cache) { response.lp_ret.barrier_cache = std::move(owned_cache); } + } else { // MIP solve auto mip_solution_ptr = @@ -200,6 +237,8 @@ std::unique_ptr call_solve( } } + pdlp_settings.barrier_cache = nullptr; + return std::make_unique(std::move(response)); } @@ -288,7 +327,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/python/cuopt/cuopt/linear_programming/data_model/data_model.py b/python/cuopt/cuopt/linear_programming/data_model/data_model.py index 7bcdfaea9b..b995815108 100644 --- a/python/cuopt/cuopt/linear_programming/data_model/data_model.py +++ b/python/cuopt/cuopt/linear_programming/data_model/data_model.py @@ -228,6 +228,24 @@ def set_objective_coefficients(self, c): """ super().set_objective_coefficients(c) + @catch_cuopt_exception + def update_q(self, c): + """ + Update the linear objective coefficients (c) for a sequence re-solve. + + Writes user-space ``c`` onto this DataModel. If a barrier cache is + present, also maps ``c`` into the cached barrier workspace and marks + it dirty (quadratic ``Q``, ``A``, and bounds must stay unchanged). + Cache reuse is QP-only: quadratic constraints take a full solve. + + Parameters + ---------- + c : array-like of float64 + Linear objective coefficients, length equal to the number of + variables on the first ``sequence_solve``. + """ + super().update_q(c) + @catch_cuopt_exception def set_objective_scaling_factor(self, objective_scaling_factor): """ diff --git a/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pxd b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pxd index 6c401b59f5..aae627a927 100644 --- a/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pxd +++ b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pxd @@ -12,6 +12,7 @@ from libcpp.memory cimport unique_ptr cdef class DataModel: cdef unique_ptr[data_model_view_t[int, double]] c_data_model_view + cdef object barrier_cache_capsule cdef void _set_cpp_quadratic_constraints( self, data_model_view_t[int, double]* c_data_model_view diff --git a/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx index b298cf57a5..de72f10267 100644 --- a/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx +++ b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx @@ -21,6 +21,14 @@ from libcpp.string cimport string from libcpp.utility cimport move from libcpp.vector cimport vector +cdef extern from "Python.h": + bint PyCapsule_IsValid(object cap, const char* name) + void* PyCapsule_GetPointer(object cap, const char* name) + +cdef extern from "cuopt/mathematical_optimization/utilities/barrier_cache.hpp" namespace "cuopt::cython": # noqa + cdef cppclass barrier_cache_t: + void update_linear_objective(const double* c, int n) except + + def type_cast(np_obj, np_type, name): if not isinstance(np_obj, np.ndarray): @@ -41,6 +49,7 @@ cdef class DataModel: def __init__(self): self.c_data_model_view.reset(new data_model_view_t[int, double]()) + self.barrier_cache_capsule = None self.maximize = False self.A_values = np.array([]) @@ -68,6 +77,14 @@ cdef class DataModel: self.row_names = np.array([]) self.quadratic_constraints = [] + def has_barrier_cache(self): + """Return whether this data model owns a reusable solver cache.""" + return self.barrier_cache_capsule is not None + + def clear_barrier_cache(self): + """Release this data model's reusable barrier cache.""" + self.barrier_cache_capsule = None + def clear_quadratic_constraints(self): self.quadratic_constraints = [] @@ -158,6 +175,34 @@ cdef class DataModel: def set_objective_coefficients(self, c): self.c = type_cast(c, np.float64, "c") + def update_q(self, c): + """Update linear objective coefficients (user-space ``c``). + + Always writes the DataModel objective. If this model owns a solver + cache from a prior Barrier solve, also crushes ``c`` into the cached + ``iteration_data_t`` and sets ``c_dirty`` so a later reuse can + skip convert/presolve. Crush runs first so a length error + leaves the DataModel coefficients unchanged. + """ + cdef barrier_cache_t* cache + cdef double[::1] c_view + new_c = type_cast(c, np.float64, "c") + if self.barrier_cache_capsule is not None: + if not PyCapsule_IsValid( + self.barrier_cache_capsule, b"cuopt.barrier_cache" + ): + raise ValueError("Invalid barrier cache stored on DataModel.") + cache = PyCapsule_GetPointer( + self.barrier_cache_capsule, + b"cuopt.barrier_cache", + ) + c_view = np.ascontiguousarray(new_c, dtype=np.float64) + if c_view.shape[0] == 0: + cache.update_linear_objective(NULL, 0) + else: + cache.update_linear_objective(&c_view[0], c_view.shape[0]) + self.c = new_c + def set_objective_scaling_factor(self, objective_scaling_factor): self.objective_scaling_factor = objective_scaling_factor diff --git a/python/cuopt/cuopt/linear_programming/solver/solver.pxd b/python/cuopt/cuopt/linear_programming/solver/solver.pxd index 04b75ce4e5..6fa11d874c 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/mathematical_optimization/utilities/cython_types.hpp" na vector[double] last_restart_duality_gap_primal_solution_ vector[double] last_restart_duality_gap_dual_solution_ +cdef extern from "cuopt/mathematical_optimization/utilities/barrier_cache.hpp" namespace "cuopt::cython": # noqa + cdef cppclass barrier_cache_t: + pass + cdef extern from "cuopt/mathematical_optimization/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/mathematical_optimization/utilities/cython_solve.hpp" na int nb_iterations_ double solve_time_ method_t solved_by_ + unique_ptr[barrier_cache_t] barrier_cache bool is_gpu() # Unified MIP solution struct — solution_ variant accessed via helpers @@ -144,6 +149,9 @@ cdef extern from "cuopt/mathematical_optimization/utilities/cython_solve.hpp" na 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, + barrier_cache_t* cache_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_wrapper.pyx b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx index 7a27a87140..a6c2a77aae 100644 --- a/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx +++ b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx @@ -18,7 +18,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 @@ -28,6 +28,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 @@ -43,6 +50,7 @@ from cuopt.linear_programming.solver.solver cimport ( linear_programming_ret_t, lp_cpu_solutions_t, lp_gpu_solutions_t, + barrier_cache_t, mip_ret_t, mip_termination_status_t, pdlp_solver_mode_t, @@ -79,6 +87,25 @@ cdef extern from "cuopt/mathematical_optimization/utilities/internals.hpp" names cdef cppclass base_solution_callback_t +cdef extern from *: + """ + #include + + static void cuopt_barrier_cache_capsule_dtor(PyObject *cap) noexcept + { + void *p = PyCapsule_GetPointer(cap, "cuopt.barrier_cache"); + if (p != nullptr) { + delete reinterpret_cast(p); + } + } + """ + void cuopt_barrier_cache_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 @@ -526,6 +553,19 @@ def prepare_solver_settings(SolverSettings settings, data_model=None, mip=False) def Solve(py_data_model_obj, SolverSettings settings, mip=False): cdef DataModel data_model_obj = py_data_model_obj + cdef barrier_cache_t* cache_in = NULL + cdef solver_ret_t* sol_ret + + if settings.sequence_solve and data_model_obj.barrier_cache_capsule is not None: + if not PyCapsule_IsValid( + data_model_obj.barrier_cache_capsule, + b"cuopt.barrier_cache", + ): + raise ValueError("Invalid barrier cache stored on DataModel.") + cache_in = PyCapsule_GetPointer( + data_model_obj.barrier_cache_capsule, + b"cuopt.barrier_cache", + ) data_model_obj.variable_types = type_cast( data_model_obj.variable_types, "S1", "variable_types" @@ -541,7 +581,22 @@ 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, + cache_in, )) + + sol_ret = sol_ret_ptr.get() + if ( + sol_ret.problem_type == ProblemCategory.LP + and sol_ret.lp_ret.barrier_cache.get() != NULL + ): + data_model_obj.barrier_cache_capsule = PyCapsule_New( + sol_ret.lp_ret.barrier_cache.release(), + b"cuopt.barrier_cache", + cuopt_barrier_cache_capsule_dtor, + ) + 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 03958d2286..02d1ec241f 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/mathematical_optimization/pdlp/solver_settings.hpp" name Barrier "cuopt::mathematical_optimization::method_t::Barrier" # noqa Unset "cuopt::mathematical_optimization::method_t::Unset" # noqa + cdef cppclass pdlp_solver_settings_t[i_t, f_t]: + bool sequence_solve + barrier_cache_t* barrier_cache + +cdef extern from "cuopt/mathematical_optimization/utilities/barrier_cache.hpp" namespace "cuopt::cython": # noqa + cdef cppclass barrier_cache_t: + pass + cdef extern from "cuopt/mathematical_optimization/solver_settings.hpp" namespace "cuopt::mathematical_optimization": # noqa cdef cppclass solver_settings_t[i_t, f_t]: @@ -90,9 +98,12 @@ cdef extern from "cuopt/mathematical_optimization/solver_settings.hpp" namespace 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 sequence_solve 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..17c3046105 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.sequence_solve = 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().sequence_solve = self.sequence_solve + + def set_sequence_solve(self, enabled): + """Enable barrier cache reuse across a sequence of solves (same sparsity).""" + self.sequence_solve = True if enabled else False + + def get_sequence_solve(self): + """Return whether sequence-solve barrier cache reuse is enabled.""" + return self.sequence_solve + def dump_parameters_to_file(self, path, hyperparameters_only=True): """Apply ``settings_dict`` / warm start to C++, then dump parameters to *path*.