diff --git a/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp b/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp index 2b12f680b5..853f7dec0f 100644 --- a/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp +++ b/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp @@ -23,6 +23,8 @@ void destroy_iteration_data(iteration_data_t* data); void apply_barrier_linear_objective(iteration_data_t& data, double const* barrier_c, int n); + +void apply_barrier_rhs(iteration_data_t& data, double const* barrier_b, int m); } // namespace cuopt::mathematical_optimization::barrier namespace cuopt { @@ -34,8 +36,8 @@ struct barrier_transform_t; * @brief GPU solve cache owned by DataModel when CUOPT_SEQUENCE_SOLVE is enabled. * * After an Optimal full solve, holds iteration_data_t and the user-barrier transform. - * update_linear_objective crushes the new linear objective and sets c_dirty so the next Solve - * reuses that workspace (skip convert/presolve/scaling). + * The update APIs crush new user data into that workspace and mark the cache dirty so the + * next Solve reuses it (skip convert/presolve/scaling). */ class barrier_cache_t { public: @@ -64,15 +66,25 @@ class barrier_cache_t { 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; + /** True when an update API has staged new data that the next Solve should reuse. */ + [[nodiscard]] bool dirty() const; + void mark_clean(); + + /** True when the last update_rhs made a row presolve dropped as empty infeasible. */ + [[nodiscard]] bool rhs_infeasible() const; /** - * Crush the input linear objective into cached iteration_data_t.c / d_c_ and set c_dirty. + * Crush the input linear objective into cached iteration_data_t.c / d_c_ and mark dirty. * Requires a stored transform and iteration_data from an Optimal solve. */ void update_linear_objective(double const* c, int n); + /** + * Crush the input constraint RHS into cached iteration_data_t.b / d_b_ and mark dirty. + * Requires a stored transform and iteration_data from an Optimal solve. + */ + void update_rhs(double const* b, int m); + private: barrier_cache_t(std::unique_ptr stream, std::unique_ptr handle); diff --git a/cpp/src/barrier/CMakeLists.txt b/cpp/src/barrier/CMakeLists.txt index 70a1f180cd..852a8fd64d 100644 --- a/cpp/src/barrier/CMakeLists.txt +++ b/cpp/src/barrier/CMakeLists.txt @@ -9,6 +9,7 @@ set(BARRIER_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/barrier_cache.cu ${CMAKE_CURRENT_SOURCE_DIR}/device_sparse_matrix.cu ${CMAKE_CURRENT_SOURCE_DIR}/pinned_host_allocator.cu + ${CMAKE_CURRENT_SOURCE_DIR}/scaling_gpu.cu ) set(BARRIER_SRC_FILES ${BARRIER_SRC_FILES} PARENT_SCOPE) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index 24dfd3c411..56bd242391 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -445,11 +445,44 @@ class barrier_reduce_helper_t { template class iteration_data_t { public: + /** The device Q to factorize: adopted from the scaling when already there, uploaded otherwise. */ + static device_csc_matrix_t make_device_Q( + std::shared_ptr> scaled_device_Q, + const lp_problem_t& lp, + const csc_matrix_t& Qin, + cuda::stream_ref stream) + { + const bool has_entries = Qin.n > 0 && Qin.col_start[Qin.n] > 0; + if (scaled_device_Q && has_entries) { + device_csc_matrix_t dQ(std::move(*scaled_device_Q)); + // All that is missing is the slack padding create_Q applies on host, which only extends + // col_start with the final nz. + const i_t old_n = dQ.n; + cuopt_assert(old_n <= Qin.n, "device Q has more columns than the host Q"); + dQ.m = dQ.n = Qin.n; + if (old_n < Qin.n) { + dQ.col_start.resize(Qin.n + 1, stream); + thrust::fill(rmm::exec_policy(stream), + dQ.col_start.begin() + old_n + 1, + dQ.col_start.end(), + dQ.nz_max); + } + return dQ; + } + if (has_entries) { return device_csc_matrix_t(Qin, stream); } + // Keep an empty but correctly shaped Q so device views are never zero-sized/uninitialized. + device_csc_matrix_t empty(stream); + empty.reset_empty(lp.num_cols, lp.num_cols, stream); + return empty; + } + iteration_data_t(const lp_problem_t& lp, i_t num_upper_bounds, const std::vector& direct_free_variables, const csc_matrix_t& Qin, - const simplex_solver_settings_t& settings) + const simplex_solver_settings_t& settings, + std::shared_ptr> scaled_device_A, + std::shared_ptr> scaled_device_Q) : upper_bounds(num_upper_bounds), c(lp.objective), b(lp.rhs), @@ -473,7 +506,6 @@ class iteration_data_t { inv_diag(lp.num_cols), inv_sqrt_diag(lp.num_cols), AD(lp.num_cols, lp.num_rows, 0), - AT(lp.num_rows, lp.num_cols, 0), ADAT(lp.num_rows, lp.num_rows, 0), // augmented(lp.num_cols + lp.num_rows, lp.num_cols + lp.num_rows, 0), A_dense(lp.num_rows, 0), @@ -482,17 +514,28 @@ class iteration_data_t { Hchol(0, 0), A(lp.A), Q(Qin), - cusparse_Q_view_(lp.handle_ptr, Q), - cusparse_view_(lp.handle_ptr, lp.A), + // Q is stored fully symmetric, so CSC(Q) is CSR(Q) and both descriptors borrow the one + // device_Q_csc_, which is declared earlier and so is already built. + cusparse_Q_view_(lp.handle_ptr, device_Q_csc_, device_Q_csc_), cusparse_info_(nullptr), + // Borrows device_A_csc_ / device_AT_csc_, both of which are declared before it and so are + // already built. + cusparse_view_(lp.handle_ptr, device_A_csc_, device_AT_csc_), 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()), device_augmented( lp.num_cols + lp.num_rows, lp.num_cols + lp.num_rows, 0, lp.handle_ptr->get_stream()), - device_A_csc_(lp.handle_ptr->get_stream()), - device_Q_csc_(lp.handle_ptr->get_stream()), - device_AT_csc_(lp.handle_ptr->get_stream()), + // Take over the scaled A when the scaling already left it on device, so it is neither + // downloaded there nor uploaded again here. + device_A_csc_(scaled_device_A + ? std::move(*scaled_device_A) + : device_csc_matrix_t(lp.A, lp.handle_ptr->get_stream())), + device_Q_csc_( + make_device_Q(std::move(scaled_device_Q), lp, Qin, lp.handle_ptr->get_stream())), + device_AT_csc_(typename device_csc_matrix_t::transposed_t{}, + device_A_csc_, + lp.handle_ptr->get_stream()), d_original_A_values(0, lp.handle_ptr->get_stream()), d_inv_diag_prime(0, lp.handle_ptr->get_stream()), d_flag_buffer(0, lp.handle_ptr->get_stream()), @@ -800,11 +843,17 @@ class iteration_data_t { if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; } - { + // AD only feeds ADAT; the augmented path derives A^T on device instead. + if (!use_augmented) { raft::common::nvtx::range scope("Barrier: LP Data: AD matrix setup"); - // Copy A into AD - AD = lp.A; - if (!use_augmented && n_dense_columns > 0) { + if (n_dense_columns == 0) { + // AD is A, which is already on device; only AD's dimensions are read on the host, so its + // entries are never materialised here. + AD.m = lp.A.m; + AD.n = lp.A.n; + } else { + // Copy A into AD + AD = lp.A; cols_to_remove.resize(lp.num_cols, 0); for (i_t k : dense_columns_unordered) { cols_to_remove[k] = 1; @@ -830,34 +879,29 @@ class iteration_data_t { A_dense.from_sparse(lp.A, j, k++); } } - - AD.transpose(AT); - } - - if (use_augmented) { - raft::common::nvtx::range scope("Barrier: augmented: device CSC upload"); - device_A_csc_.copy(A, handle_ptr->get_stream()); - device_AT_csc_.copy(AT, handle_ptr->get_stream()); - if (Q.n > 0 && Q.col_start[Q.n] > 0) { - device_Q_csc_.copy(Q, handle_ptr->get_stream()); - } else { - // Keep an empty but correctly shaped Q so device views are never zero-sized/uninitialized. - device_Q_csc_.reset_empty(A.n, A.n, handle_ptr->get_stream()); - } } // device_AD / device_A / ADAT path is only used when forming ADAT (!use_augmented). if (!use_augmented) { raft::common::nvtx::range scope("Barrier: LP Data: device AD path"); - device_AD.copy(AD, handle_ptr->get_stream()); - d_original_A_values.resize(device_AD.x.size(), handle_ptr->get_stream()); - raft::copy(d_original_A_values.data(), - device_AD.x.data(), - device_AD.x.size(), - handle_ptr->get_stream()); + if (n_dense_columns > 0) { + device_AD.copy(AD, handle_ptr->get_stream()); + // AD differs from A once dense columns are dropped, so form_adat needs its own snapshot + // of the unscaled values to restore from. + d_original_A_values.resize(device_AD.x.size(), handle_ptr->get_stream()); + raft::copy(d_original_A_values.data(), + device_AD.x.data(), + device_AD.x.size(), + handle_ptr->get_stream()); + device_AD.to_compressed_row(device_A, handle_ptr->get_stream()); + } else { + // AD == A, so device_AD is seeded straight from device_A_csc_, which also doubles as + // form_adat's restore source, and device_AT_csc_ (already CSR(A)) serves as the SpGEMM's + // left operand -- neither needs a second copy. Both stay read-only for the whole solve. + device_AD.copy(device_A_csc_, handle_ptr->get_stream()); + } // For efficient scaling of AD col we form the col index array device_AD.form_col_index(handle_ptr->get_stream()); - device_AD.to_compressed_row(device_A, handle_ptr->get_stream()); RAFT_CHECK_CUDA(handle_ptr->get_stream().get()); } @@ -895,7 +939,7 @@ class iteration_data_t { } // Attach this solve's settings and rewind iterate-dependent state so barrier can - // start with the new c. A and Q are unchanged; the previous solve + // start with the new c / b. 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 reset_iterate_state(const simplex_solver_settings_t& settings) @@ -1066,7 +1110,6 @@ class iteration_data_t { stream_view_); settings_.log.debug("augmented nz %d (gpu build)\n", total_nnz); - cuopt_assert(A.col_start[n] == AT.col_start[m], "A nz != AT nz"); handle_ptr->sync_stream(); #ifdef CHECK_SYMMETRY @@ -1154,10 +1197,11 @@ 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(), - handle_ptr->get_stream()); + // device_A_csc_ holds A's unscaled values and is never written, so when AD == A it is the + // snapshot; with dense columns removed AD differs and carries its own. + const f_t* original_values = + n_dense_columns > 0 ? d_original_A_values.data() : device_A_csc_.x.data(); + raft::copy(device_AD.x.data(), original_values, device_AD.x.size(), handle_ptr->get_stream()); } { raft::common::nvtx::range scope("Barrier: Form ADAT: inv_diag prime"); @@ -1199,12 +1243,29 @@ class iteration_data_t { if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { return; } if (first_call) { raft::common::nvtx::range scope("Barrier: Form ADAT: cusparse init"); + // With no dense columns AD == A, so CSC(A^T) is the CSR(A) the SpGEMM needs and + // device_AT_csc_ is used directly instead of keeping a second copy in device_A. + const bool own_csr = n_dense_columns > 0; + const i_t A_rows = own_csr ? device_A.m : device_AT_csc_.n; + const i_t A_cols = own_csr ? device_A.n : device_AT_csc_.m; + const i_t A_nnz = own_csr ? device_A.nz_max : device_AT_csc_.nz_max; + i_t* A_offsets = own_csr ? device_A.row_start.data() : device_AT_csc_.col_start.data(); + i_t* A_indices = own_csr ? device_A.j.data() : device_AT_csc_.i.data(); + f_t* A_values = own_csr ? device_A.x.data() : device_AT_csc_.x.data(); try { if (!cusparse_info_) { cusparse_info_ = std::make_unique>(handle_ptr); } - initialize_cusparse_data( - handle_ptr, device_A, device_AD, device_ADAT, spgemm_info()); + initialize_cusparse_data(handle_ptr, + A_rows, + A_cols, + A_nnz, + A_offsets, + A_indices, + A_values, + device_AD, + device_ADAT, + spgemm_info()); } catch (const raft::cuda_error& e) { settings_.log.printf("Error in initialize_cusparse_data: %s\n", e.what()); return; @@ -1214,7 +1275,7 @@ class iteration_data_t { { raft::common::nvtx::range scope("Barrier: Form ADAT: ADAT multiply"); - multiply_kernels(handle_ptr, device_A, device_AD, device_ADAT, spgemm_info()); + multiply_kernels(handle_ptr, device_ADAT, spgemm_info()); handle_ptr->sync_stream(); } @@ -2210,7 +2271,6 @@ class iteration_data_t { rmm::device_uvector d_original_A_values; csc_matrix_t AD; - csc_matrix_t AT; csc_matrix_t ADAT; // csc_matrix_t augmented; device_csr_matrix_t device_augmented; @@ -2431,10 +2491,18 @@ void cholesky_debug_check(const iteration_data_t& data, } template -barrier_solver_t::barrier_solver_t(const lp_problem_t& lp, - const simplex::presolve_info_t& presolve, - const simplex_solver_settings_t& settings) - : lp(lp), settings(settings), presolve_info(presolve), stream_view_(lp.handle_ptr->get_stream()) +barrier_solver_t::barrier_solver_t( + const lp_problem_t& lp, + const simplex::presolve_info_t& presolve, + const simplex_solver_settings_t& settings, + std::shared_ptr> device_A, + std::shared_ptr> device_Q) + : lp(lp), + settings(settings), + presolve_info(presolve), + stream_view_(lp.handle_ptr->get_stream()), + device_A_(std::move(device_A)), + device_Q_(std::move(device_Q)) { } @@ -4863,8 +4931,13 @@ lp_status_t barrier_solver_t::solve( Qin = xf->barrier_Q.get(); } if (lp.Q.n > 0) { create_Q(lp, *Qin); } - owned_data = std::make_unique>( - lp, num_upper_bounds, presolve_info.direct_free_variables, *Qin, settings); + owned_data = std::make_unique>(lp, + num_upper_bounds, + presolve_info.direct_free_variables, + *Qin, + settings, + std::move(device_A_), + std::move(device_Q_)); lp_status_t status = barrier_advanced_solve(start_time, solution, *owned_data); return store_or_clear_cache(cache, owned_data, status); } catch (const raft::cuda_error& e) { @@ -4911,6 +4984,18 @@ void apply_barrier_linear_objective(iteration_data_t& data, data.d_c_.data(), data.c.data(), static_cast(n), data.handle_ptr->get_stream()); } +void apply_barrier_rhs(iteration_data_t& data, double const* barrier_b, int m) +{ + cuopt_expects( + barrier_b != nullptr && static_cast(data.b.size()) == m && + static_cast(data.d_b_.size()) == m, + error_type_t::ValidationError, + "update_rhs: barrier RHS size does not match cached iteration_data_t."); + std::copy(barrier_b, barrier_b + m, data.b.data()); + raft::copy( + data.d_b_.data(), data.b.data(), static_cast(m), data.handle_ptr->get_stream()); +} + #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE template bool validate_barrier_cone_layout( const lp_problem_t& problem, const simplex_solver_settings_t& settings); diff --git a/cpp/src/barrier/barrier.hpp b/cpp/src/barrier/barrier.hpp index 58fb6d63f5..50d80fbfd5 100644 --- a/cpp/src/barrier/barrier.hpp +++ b/cpp/src/barrier/barrier.hpp @@ -20,6 +20,7 @@ #include #include +#include #include namespace cuopt::mathematical_optimization { @@ -36,12 +37,21 @@ bool validate_barrier_cone_layout(const simplex::lp_problem_t& problem template class iteration_data_t; // Forward declare +template +class device_csc_matrix_t; // Forward declare + template class barrier_solver_t { public: + // `device_A` / `device_Q` are the scaled matrices the GPU scaling already left on device, taken + // over here so they are neither downloaded there nor uploaded again. Null means the solver + // uploads them from `lp`. Only solve() consumes them; solve_with_cache() reuses the cached + // iteration_data_t and never looks at them. barrier_solver_t(const simplex::lp_problem_t& lp, const simplex::presolve_info_t& presolve, - const simplex::simplex_solver_settings_t& settings); + const simplex::simplex_solver_settings_t& settings, + std::shared_ptr> device_A = nullptr, + std::shared_ptr> device_Q = nullptr); simplex::lp_status_t solve(f_t start_time, simplex::lp_solution_t& solution, cuopt::mathematical_optimization::barrier_cache_t* cache = nullptr); @@ -117,6 +127,9 @@ class barrier_solver_t { const simplex::simplex_solver_settings_t& settings; const simplex::presolve_info_t& presolve_info; cuda::stream_ref stream_view_; + // Handed over to iteration_data_t by solve(), which empties them. + std::shared_ptr> device_A_; + std::shared_ptr> device_Q_; }; } // namespace cuopt::mathematical_optimization::barrier diff --git a/cpp/src/barrier/barrier_cache.cu b/cpp/src/barrier/barrier_cache.cu index f4577bc839..7af14bf8c4 100644 --- a/cpp/src/barrier/barrier_cache.cu +++ b/cpp/src/barrier/barrier_cache.cu @@ -22,6 +22,29 @@ using barrier_iteration_data_t = barrier::iteration_data_t; using barrier_iteration_data_ptr = std::unique_ptr; +static void require_warm_cache(barrier_transform_t const* transform, + barrier_iteration_data_t const* data, + char const* api) +{ + cuopt_expects(transform != nullptr, + error_type_t::ValidationError, + "%s: no barrier transform; Solve with CUOPT_SEQUENCE_SOLVE enabled first.", + api); + cuopt_expects(data != nullptr, + error_type_t::ValidationError, + "%s: no cached iteration_data; Solve a QP to Optimal first.", + api); +} + +// Re-adds the first solve's barrier-minus-crush shift so the update lands in cached coordinates. +static void add_shift(std::vector& crushed, std::vector const& shift) +{ + if (shift.size() != crushed.size()) { return; } + for (std::size_t i = 0; i < crushed.size(); ++i) { + crushed[i] += shift[i]; + } +} + struct barrier_cache_t::impl { impl(std::unique_ptr stream_in, std::unique_ptr handle_in) : stream(std::move(stream_in)), @@ -36,6 +59,8 @@ struct barrier_cache_t::impl { std::unique_ptr transform; barrier_iteration_data_ptr iteration_data; bool c_dirty{false}; + bool b_dirty{false}; + bool rhs_infeasible{false}; }; barrier_cache_t::barrier_cache_t(std::unique_ptr stream, @@ -66,7 +91,7 @@ void barrier_cache_t::clear() { impl_->iteration_data.reset(); impl_->transform.reset(); - impl_->c_dirty = false; + mark_clean(); } void barrier_cache_t::store_iteration_data(barrier_iteration_data_t* data) @@ -88,22 +113,25 @@ 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::dirty() const +{ + return (impl_->c_dirty || impl_->b_dirty) && impl_->transform != nullptr && + impl_->iteration_data.get() != nullptr; +} -bool barrier_cache_t::c_dirty() const +void barrier_cache_t::mark_clean() { - return impl_->c_dirty && impl_->transform != nullptr && impl_->iteration_data.get() != nullptr; + impl_->c_dirty = false; + impl_->b_dirty = false; + impl_->rhs_infeasible = false; } +bool barrier_cache_t::rhs_infeasible() const { return impl_->rhs_infeasible; } + void barrier_cache_t::update_linear_objective(double const* c, int n) { - cuopt_expects(impl_->transform != nullptr, - error_type_t::ValidationError, - "update_linear_objective: no barrier transform; Solve with CUOPT_SEQUENCE_SOLVE " - "enabled first."); - cuopt_expects(impl_->iteration_data.get() != nullptr, - error_type_t::ValidationError, - "update_linear_objective: no cached iteration_data; Solve a QP to Optimal first."); + require_warm_cache( + impl_->transform.get(), impl_->iteration_data.get(), "update_linear_objective"); // Cached Q and c are in minimization space. std::vector user_objective; if (impl_->transform->maximize && c != nullptr && n > 0) { @@ -135,11 +163,7 @@ void barrier_cache_t::update_linear_objective(double const* c, int n) } barrier_lp.obj_constant += obj_constant_delta; } - if (linear_obj_shift.size() == crushed.size()) { - for (std::size_t j = 0; j < crushed.size(); ++j) { - crushed[j] += linear_obj_shift[j]; - } - } + add_shift(crushed, linear_obj_shift); // The next solve builds its solver from barrier_lp, so keep its objective and the cached // iteration workspace on the same c. auto& barrier_objective = barrier_lp.objective; @@ -153,4 +177,33 @@ void barrier_cache_t::update_linear_objective(double const* c, int n) impl_->c_dirty = true; } +void barrier_cache_t::update_rhs(double const* b, int m) +{ + require_warm_cache(impl_->transform.get(), impl_->iteration_data.get(), "update_rhs"); + std::vector crushed; + try { + crushed = crush_user_rhs(*impl_->transform, b, m); + } catch (update_rhs_infeasible_error const&) { + // Cache stays usable for a later feasible update; the next Solve reports INFEASIBLE from + // this flag without running IPM. + impl_->rhs_infeasible = true; + impl_->b_dirty = true; + return; + } catch (std::invalid_argument const& e) { + cuopt_expects(false, error_type_t::ValidationError, "%s", e.what()); + } + impl_->rhs_infeasible = false; + add_shift(crushed, impl_->transform->rhs_shift); + // barrier_lp->rhs also seeds the next solve's Mehrotra start, so keep it and the cached + // workspace on the same b. + auto& barrier_rhs = impl_->transform->barrier_lp->rhs; + cuopt_expects(barrier_rhs.size() == crushed.size(), + error_type_t::ValidationError, + "update_rhs: crushed RHS size does not match the cached barrier LP."); + barrier_rhs = crushed; + barrier::apply_barrier_rhs( + *impl_->iteration_data, crushed.data(), static_cast(crushed.size())); + impl_->b_dirty = true; +} + } // namespace cuopt::mathematical_optimization diff --git a/cpp/src/barrier/barrier_transform.hpp b/cpp/src/barrier/barrier_transform.hpp index c4c96c9f1b..73cefa86a0 100644 --- a/cpp/src/barrier/barrier_transform.hpp +++ b/cpp/src/barrier/barrier_transform.hpp @@ -8,19 +8,36 @@ #pragma once #include +#include #include +#include +#include +#include #include #include +#include +#include #include namespace cuopt::mathematical_optimization { +/** + * Singleton rows that force a cone head nonnegative. The SOC expansion requires every head to be + * provably >= 0, and proves it from these rows when the head has no explicit bound, so a new RHS + * can invalidate a cached expansion. Only heads that need the proof are recorded. + */ +struct cone_head_bound_t { + int head_col{0}; + // (row, coefficient) pairs, each implying head >= rhs[row] / coefficient. + std::vector> rows; +}; + /** * User-to-barrier transform retained on barrier_cache_t after Optimal: * convert / presolve / scaling, plus the scaled LP. - * Enough to crush a new linear objective from the original problem into barrier - * coordinates and to uncrush a solution without rerunning those algorithms. + * Enough to crush new linear objective or RHS data from the original problem into + * barrier coordinates and to uncrush a solution without rerunning those algorithms. */ struct barrier_transform_t { int user_num_cols{0}; @@ -34,28 +51,168 @@ struct barrier_transform_t { // Enough of the user problem for reuse uncrush without rebuilding A. std::vector row_sense; int cone_var_start{0}; + // cone_var_start after convert, which inserts inequality slacks ahead of the cone block. + int converted_cone_var_start{0}; std::vector second_order_cone_dims; - int expanded_original_num_cols{0}; + // Quadratic constraint count the cached expansion was built from. + int num_quadratic_constraints{0}; + // Dimensions before the QCMATRIX->SOC expansion, which permutes columns and appends rows. + // Updates arrive in these coordinates, not the expanded ones. Zero when no expansion ran. + int pre_expansion_num_cols{0}; + int pre_expansion_num_rows{0}; std::vector original_col_to_expanded_col; + // RHS of the rows the expansion appended. Fixed by the quadratic constraints, so an RHS + // update keeps them and only overwrites the model's own rows. + std::vector cone_row_rhs; + std::vector cone_head_bounds; 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*ell shift, etc.). std::vector linear_obj_shift; + // Barrier RHS minus crush(user b) from the first solve (fixed/lower-bound shifts). + std::vector rhs_shift; + // False when range rows or folding put the user RHS somewhere other than barrier_lp->rhs. + bool rhs_update_supported{false}; + // Absolute primal tolerance of the first solve, used to test rows presolve dropped as empty. + double primal_tol{1e-6}; std::unique_ptr> barrier_lp; // CSC Q with slack columns, as consumed by iteration_data_t. Not the same object as // barrier_lp->Q. std::unique_ptr> barrier_Q; }; +// Dimensions an update is sized in: the cached user counts, or the smaller pre-expansion counts +// when the QCMATRIX->SOC expansion grew the problem. +inline int model_num_cols(barrier_transform_t const& xf) +{ + return xf.pre_expansion_num_cols > 0 ? xf.pre_expansion_num_cols : xf.user_num_cols; +} + +inline int model_num_rows(barrier_transform_t const& xf) +{ + return xf.pre_expansion_num_rows > 0 ? xf.pre_expansion_num_rows : xf.user_num_rows; +} + +// The expansion permutes model columns into a [linear | cone] layout. An empty map means no +// expansion ran, so the layouts coincide. +inline int model_col_to_expanded_col(barrier_transform_t const& xf, int model_col) +{ + return xf.original_col_to_expanded_col.empty() + ? model_col + : xf.original_col_to_expanded_col[static_cast(model_col)]; +} + +// convert inserts the inequality slacks ahead of the cone block, pushing every cone column +// right. Mirrors user_col_to_problem_col in presolve.cpp. +inline int expanded_col_to_converted_col(barrier_transform_t const& xf, int expanded_col) +{ + if (xf.second_order_cone_dims.empty() || xf.converted_cone_var_start <= xf.cone_var_start || + expanded_col < xf.cone_var_start) { + return expanded_col; + } + return xf.converted_cone_var_start + (expanded_col - xf.cone_var_start); +} + +// Move the model's own coefficients into the expanded layout the cached problem is sized for. +template +std::vector scatter_model_objective(barrier_transform_t const& xf, + std::vector const& model_objective) +{ + std::vector expanded(static_cast(xf.user_num_cols), f_t(0)); + for (int j = 0; j < static_cast(model_objective.size()); ++j) { + expanded[static_cast(model_col_to_expanded_col(xf, j))] = + model_objective[static_cast(j)]; + } + return expanded; +} + +// Inverse of scatter_model_objective, giving crush_user_linear_objective the model-sized input +// it expects. The expansion adds no objective coefficients, so nothing is lost. +template +std::vector gather_model_objective(barrier_transform_t const& xf, + std::vector const& expanded_objective) +{ + std::vector model_objective(static_cast(model_num_cols(xf))); + for (int j = 0; j < static_cast(model_objective.size()); ++j) { + model_objective[static_cast(j)] = static_cast( + expanded_objective[static_cast(model_col_to_expanded_col(xf, j))]); + } + return model_objective; +} + +// Reuse never re-runs the expansion, so the cone block must be laid out exactly as the cached +// one left it. +template +bool cone_layout_matches(barrier_transform_t const& xf, + simplex::user_problem_t const& user_problem) +{ + return user_problem.cone_var_start == static_cast(xf.cone_var_start) && + user_problem.second_order_cone_dims.size() == xf.second_order_cone_dims.size() && + std::equal(user_problem.second_order_cone_dims.begin(), + user_problem.second_order_cone_dims.end(), + xf.second_order_cone_dims.begin(), + [](i_t dim, int cached) { return dim == static_cast(cached); }); +} + +// A cone head with no explicit nonnegative bound is only admissible because some singleton row +// forces it nonnegative. The expansion checks that once; collect the rows it relied on so an RHS +// update can re-check them against the new RHS. +template +std::vector record_cone_head_bounds( + simplex::user_problem_t const& user_problem) +{ + std::vector bounds; + if (user_problem.second_order_cone_dims.empty()) { return bounds; } + + // Cones only reach here via the expansion, which always sets original_num_rows. + const auto& A = user_problem.A; + const i_t model_rows = user_problem.original_num_rows; + std::vector row_nz(model_rows, 0); + for (i_t j = 0; j < user_problem.num_cols; ++j) { + for (i_t p = A.col_start[j]; p < A.col_start[j + 1]; ++p) { + if (A.i[p] < model_rows) { ++row_nz[A.i[p]]; } + } + } + + // Only heads that were already model variables carry the precondition. A head the expansion + // created is nonnegative by cone membership, so no row has to prove it. + std::vector is_model_col(user_problem.num_cols, 0); + for (i_t expanded : user_problem.original_col_to_expanded_col) { + if (expanded >= 0 && expanded < user_problem.num_cols) { is_model_col[expanded] = 1; } + } + + i_t head = user_problem.cone_var_start; + for (i_t q_k : user_problem.second_order_cone_dims) { + if (head < 0 || head >= user_problem.num_cols) { break; } + if (is_model_col[head] && !(user_problem.lower[head] >= f_t(0))) { + cone_head_bound_t bound; + bound.head_col = static_cast(head); + // A is CSC, so the head's own column already lists every row it appears in. + for (i_t p = A.col_start[head]; p < A.col_start[head + 1]; ++p) { + const i_t i = A.i[p]; + if (i >= model_rows || row_nz[i] != 1) { continue; } + const f_t a = A.x[p]; + const char sense = user_problem.row_sense[i]; + if ((sense == 'G' && a > f_t(0)) || (sense == 'L' && a < f_t(0))) { + bound.rows.emplace_back(static_cast(i), static_cast(a)); + } + } + bounds.push_back(std::move(bound)); + } + head += q_k; + } + return bounds; +} + 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) { + if (c == nullptr || n != model_num_cols(xf)) { throw std::invalid_argument( - "update_linear_objective: linear objective length must match the cached user column count."); + "update_linear_objective: linear objective length must match the cached model column count."); } if (xf.original_num_cols < xf.user_num_cols) { throw std::invalid_argument( @@ -64,10 +221,22 @@ inline std::vector crush_user_linear_objective(barrier_transform_t const if (xf.barrier_lp == nullptr) { throw std::invalid_argument("update_linear_objective: cached barrier LP is missing."); } + if (!xf.original_col_to_expanded_col.empty() && + static_cast(xf.original_col_to_expanded_col.size()) != n) { + throw std::invalid_argument( + "update_linear_objective: cached column map does not cover the model columns."); + } + // The expansion leaves the variables it adds out of the objective, so only the positions of + // the model's own coefficients move. std::vector orig(static_cast(xf.original_num_cols), 0.0); for (int j = 0; j < n; ++j) { - orig[static_cast(j)] = c[j]; + int const converted_col = expanded_col_to_converted_col(xf, model_col_to_expanded_col(xf, j)); + if (converted_col < 0 || converted_col >= xf.original_num_cols) { + throw std::invalid_argument( + "update_linear_objective: cached column map points outside the converted problem."); + } + orig[static_cast(converted_col)] = c[j]; } for (int j : xf.presolve_info.negated_variables) { orig[static_cast(j)] *= -1.0; @@ -109,4 +278,99 @@ inline std::vector crush_user_linear_objective(barrier_transform_t const return presolved; } +// Distinct from the invalid_argument cases so the caller can report INFEASIBLE rather than a +// validation failure. +struct update_rhs_infeasible_error : std::runtime_error { + explicit update_rhs_infeasible_error(std::string const& message) : std::runtime_error(message) {} +}; + +inline std::vector crush_user_rhs(barrier_transform_t const& xf, double const* b, int m) +{ + if (b == nullptr || m != model_num_rows(xf)) { + throw std::invalid_argument("update_rhs: RHS length must match the cached model row count."); + } + if (!xf.rhs_update_supported) { + throw std::invalid_argument( + "update_rhs: cached convert used range rows or folding; run a full Solve."); + } + if (xf.original_num_rows != xf.user_num_rows) { + throw std::invalid_argument( + "update_rhs: cached original row count does not match the user row count."); + } + if (static_cast(xf.row_sense.size()) != xf.user_num_rows) { + throw std::invalid_argument( + "update_rhs: cached row-sense count does not match the user row count."); + } + if (xf.barrier_lp == nullptr) { + throw std::invalid_argument("update_rhs: cached barrier LP is missing."); + } + + // The quadratic constraints fix the RHS of the appended rows, so an update overwrites the + // model's own rows and keeps the cached tail. The tail is empty without an expansion. + std::vector expanded(b, b + m); + expanded.insert(expanded.end(), xf.cone_row_rhs.begin(), xf.cone_row_rhs.end()); + if (static_cast(expanded.size()) != xf.user_num_rows) { + throw std::invalid_argument("update_rhs: cached cone-row RHS does not span the expanded rows."); + } + + // The expansion proved these heads nonnegative from the old RHS. A full solve rejects the + // model once that no longer holds, so re-prove it here rather than trust the cached verdict. + for (cone_head_bound_t const& bound : xf.cone_head_bounds) { + double implied = -std::numeric_limits::infinity(); + for (auto const& [row, coefficient] : bound.rows) { + implied = std::max(implied, expanded[static_cast(row)] / coefficient); + } + if (!(implied >= 0.0)) { + throw std::invalid_argument( + "update_rhs: new RHS no longer implies second-order cone head variable " + + std::to_string(bound.head_col) + " is nonnegative."); + } + } + + // convert turns 'G' rows into 'L' rows by negating the row and its RHS. + std::vector original(static_cast(xf.original_num_rows)); + for (int i = 0; i < xf.user_num_rows; ++i) { + original[static_cast(i)] = + xf.row_sense[static_cast(i)] == 'G' ? -expanded[i] : expanded[i]; + } + + // Dropped rows were empty, so the new RHS never reaches the barrier: 'E' needs 0 == b_i and + // the rest need 0 <= b_i. + for (int i : xf.presolve_info.removed_constraints) { + if (i < 0 || i >= xf.user_num_rows) { + throw std::invalid_argument("update_rhs: removed constraint index is out of range."); + } + double const converted_rhs = original[static_cast(i)]; + bool const infeasible = xf.row_sense[static_cast(i)] == 'E' + ? std::abs(converted_rhs) > xf.primal_tol + : converted_rhs < -xf.primal_tol; + if (infeasible) { + throw update_rhs_infeasible_error("update_rhs: empty constraint row " + std::to_string(i) + + " is infeasible with the new RHS."); + } + } + + // Empty remaining_constraints means either no empty-row pass ran, or every row was dropped + // and accepted above. + std::vector presolved; + if (!xf.presolve_info.remaining_constraints.empty()) { + presolved.resize(xf.presolve_info.remaining_constraints.size()); + for (std::size_t k = 0; k < xf.presolve_info.remaining_constraints.size(); ++k) { + presolved[k] = original[static_cast(xf.presolve_info.remaining_constraints[k])]; + } + } else if (xf.presolve_info.removed_constraints.empty()) { + presolved = std::move(original); + } + + if (static_cast(presolved.size()) != xf.barrier_lp->num_rows || + xf.row_scales.size() != presolved.size()) { + throw std::invalid_argument( + "update_rhs: crushed RHS size does not match barrier rows / row_scales."); + } + for (std::size_t i = 0; i < presolved.size(); ++i) { + presolved[i] /= xf.row_scales[i]; + } + return presolved; +} + } // namespace cuopt::mathematical_optimization diff --git a/cpp/src/barrier/cusparse_view.cu b/cpp/src/barrier/cusparse_view.cu index c7a9cac067..66d59d0e1d 100644 --- a/cpp/src/barrier/cusparse_view.cu +++ b/cpp/src/barrier/cusparse_view.cu @@ -23,6 +23,8 @@ #include +#include + namespace cuopt::mathematical_optimization::barrier { #define CUDA_VER_12_4_UP (CUDART_VERSION >= 12040) @@ -181,25 +183,26 @@ cusparse_view_t::cusparse_view_t(raft::handle_t const* handle_ptr, RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsesetpointermode(handle_ptr->get_cusparse_handle(), CUSPARSE_POINTER_MODE_DEVICE, handle_ptr->get_stream().get())); - // TMP matrix data should already be on the GPU constexpr bool debug = false; if (debug) { printf("A hash: %zu\n", A.hash()); } - csr_matrix_t A_csr(A.m, A.n, 1); - A.to_compressed_row(A_csr); - rows_ = A_csr.m; - i_t cols = A_csr.n; - i_t nnz = A_csr.x.size(); - const std::vector& offsets = A_csr.row_start; - const std::vector& indices = A_csr.j; - const std::vector& data = A_csr.x; - - A_offsets_ = device_copy(offsets, handle_ptr->get_stream()); - A_indices_ = device_copy(indices, handle_ptr->get_stream()); - A_data_ = device_copy(data, handle_ptr->get_stream()); - - A_T_offsets_ = device_copy(A.col_start, handle_ptr->get_stream()); - A_T_indices_ = device_copy(A.i, handle_ptr->get_stream()); - A_T_data_ = device_copy(A.x, handle_ptr->get_stream()); + + // A^T's CSR is A's CSC verbatim, so one upload serves the transpose view; the forward CSR is + // then derived from it on device instead of being converted on the host and uploaded again. + device_csc_matrix_t d_A(A, handle_ptr->get_stream()); + device_csr_matrix_t d_A_csr(handle_ptr->get_stream()); + d_A.to_compressed_row(d_A_csr, handle_ptr->get_stream()); + + rows_ = A.m; + const i_t cols = A.n; + const i_t nnz = A.col_start[A.n]; + + A_offsets_ = std::move(d_A_csr.row_start); + A_indices_ = std::move(d_A_csr.j); + A_data_ = std::move(d_A_csr.x); + + A_T_offsets_ = std::move(d_A.col_start); + A_T_indices_ = std::move(d_A.i); + A_T_data_ = std::move(d_A.x); A_ = pdlp::make_csr( rows_, cols, nnz, A_offsets_.data(), A_indices_.data(), A_data_.data()); @@ -217,6 +220,49 @@ cusparse_view_t::cusparse_view_t(raft::handle_t const* handle_ptr, A_T_.get(), y.get(), x.get(), spmv_buffer_transpose_, A_T_offsets_.size() - 1); } +template +cusparse_view_t::cusparse_view_t(raft::handle_t const* handle_ptr, + device_csc_matrix_t& A_csc, + device_csc_matrix_t& AT_csc) + : handle_ptr_(handle_ptr), + A_offsets_(0, handle_ptr->get_stream()), + A_indices_(0, handle_ptr->get_stream()), + A_data_(0, handle_ptr->get_stream()), + A_T_offsets_(0, handle_ptr->get_stream()), + A_T_indices_(0, handle_ptr->get_stream()), + A_T_data_(0, handle_ptr->get_stream()), + spmv_buffer_(0, handle_ptr->get_stream()), + spmv_buffer_transpose_(0, handle_ptr->get_stream()), + d_one_(one_v, handle_ptr->get_stream()), + d_minus_one_(neg_one_v, handle_ptr->get_stream()), + d_zero_(zero_v, handle_ptr->get_stream()) +{ + RAFT_CUBLAS_TRY(raft::linalg::detail::cublassetpointermode( + handle_ptr->get_cublas_handle(), CUBLAS_POINTER_MODE_DEVICE, handle_ptr->get_stream().get())); + RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsesetpointermode(handle_ptr->get_cusparse_handle(), + CUSPARSE_POINTER_MODE_DEVICE, + handle_ptr->get_stream().get())); + rows_ = A_csc.m; + const i_t cols = A_csc.n; + const i_t nnz = A_csc.nz_max; + + // Both descriptors are relabellings of the caller's buffers: CSC(A^T) is CSR(A), and CSC(A) + // is CSR(A^T). + A_ = pdlp::make_csr( + rows_, cols, nnz, AT_csc.col_start.data(), AT_csc.i.data(), AT_csc.x.data()); + A_T_ = pdlp::make_csr( + cols, rows_, nnz, A_csc.col_start.data(), A_csc.i.data(), A_csc.x.data()); + + // Temporary vectors used to initialize the SpMV buffers and preprocessing data. + rmm::device_uvector d_x(cols, handle_ptr_->get_stream()); + rmm::device_uvector d_y(rows_, handle_ptr_->get_stream()); + auto x = pdlp::make_dnvec(d_x.size(), d_x.data()); + auto y = pdlp::make_dnvec(d_y.size(), d_y.data()); + + init_spmv_buffer_and_preprocess(A_.get(), x.get(), y.get(), spmv_buffer_, rows_); + init_spmv_buffer_and_preprocess(A_T_.get(), y.get(), x.get(), spmv_buffer_transpose_, cols); +} + template pdlp::cusparse_dn_vec_uptr 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 8a3c38e33c..573b86efe8 100644 --- a/cpp/src/barrier/cusparse_view.hpp +++ b/cpp/src/barrier/cusparse_view.hpp @@ -30,6 +30,13 @@ class cusparse_view_t { // TMP matrix data should already be on the GPU and in CSR not CSC cusparse_view_t(raft::handle_t const* handle_ptr, const csc_matrix_t& A); + // Borrowing overload: the descriptors point at caller-owned device buffers, which must outlive + // this view, and the A_* / A_T_* members below stay empty. A_csc supplies the transpose view + // (CSC(A) is CSR(A^T)) and AT_csc the forward view (CSC(A^T) is CSR(A)). + cusparse_view_t(raft::handle_t const* handle_ptr, + device_csc_matrix_t& A_csc, + device_csc_matrix_t& AT_csc); + pdlp::cusparse_dn_vec_uptr create_vector(rmm::device_uvector const& vec); template diff --git a/cpp/src/barrier/device_sparse_matrix.cuh b/cpp/src/barrier/device_sparse_matrix.cuh index 7595e85fa2..f647dec0ad 100644 --- a/cpp/src/barrier/device_sparse_matrix.cuh +++ b/cpp/src/barrier/device_sparse_matrix.cuh @@ -194,6 +194,7 @@ class device_csc_matrix_t { { } + /** Move leaves the source empty; needed to hand a matrix over without a device-to-device copy. */ 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; @@ -242,6 +243,20 @@ class device_csc_matrix_t { raft::copy(x.data(), A.x.data(), A.x.size(), stream); } + /** Copy from another device CSC matrix, without going through the host. */ + void copy(const device_csc_matrix_t& A, cuda::stream_ref stream) + { + m = A.m; + n = A.n; + nz_max = A.nz_max; + col_start.resize(A.col_start.size(), stream); + raft::copy(col_start.data(), A.col_start.data(), A.col_start.size(), stream); + i.resize(A.i.size(), stream); + raft::copy(i.data(), A.i.data(), A.i.size(), stream); + x.resize(A.x.size(), stream); + raft::copy(x.data(), A.x.data(), A.x.size(), stream); + } + /** Reset to an empty (all-zero col_start, no nonzeros) matrix of the given shape. */ void reset_empty(i_t rows, i_t cols, cuda::stream_ref stream) { @@ -256,6 +271,19 @@ class device_csc_matrix_t { * device. */ void to_compressed_row(device_csr_matrix_t& Arow, cuda::stream_ref stream) const; + /** Same semantics as csc_matrix_t::transpose, entirely on device. */ + void transpose(device_csc_matrix_t& AT, cuda::stream_ref stream) const; + + /** Tag selecting the transpose constructor below. */ + struct transposed_t {}; + + /** Construct as A^T, entirely on device. */ + device_csc_matrix_t(transposed_t, const device_csc_matrix_t& A, cuda::stream_ref stream) + : col_start(0, stream), i(0, stream), x(0, stream), col_index(0, stream) + { + A.transpose(*this, stream); + } + void form_col_index(cuda::stream_ref stream) { col_index.resize(x.size(), stream); @@ -386,6 +414,20 @@ class device_csr_matrix_t { raft::copy(x.data(), A.x.data(), A.x.size(), stream); } + /** Copy from a device CSC matrix holding this matrix's transpose; the arrays are identical. */ + void copy_transposed(const device_csc_matrix_t& AT, cuda::stream_ref stream) + { + m = AT.n; + n = AT.m; + nz_max = AT.nz_max; + row_start.resize(AT.col_start.size(), stream); + raft::copy(row_start.data(), AT.col_start.data(), AT.col_start.size(), stream); + j.resize(AT.i.size(), stream); + raft::copy(j.data(), AT.i.data(), AT.i.size(), stream); + x.resize(AT.x.size(), stream); + raft::copy(x.data(), AT.x.data(), AT.x.size(), stream); + } + i_t nz_max; // maximum number of entries i_t m; // number of rows i_t n; // number of columns @@ -397,86 +439,164 @@ class device_csr_matrix_t { // to avoid extra space / computation) }; +// One block per CSC column; each nonzero claims its CSR slot through its row's atomic cursor. template -void device_csc_matrix_t::to_compressed_row(device_csr_matrix_t& Arow, - cuda::stream_ref stream) const +__global__ void csc_to_csr_scatter_kernel(i_t n_cols, + const i_t* __restrict__ col_start, + const i_t* __restrict__ row_ind, + const f_t* __restrict__ csc_val, + i_t* __restrict__ next_pos, + i_t* __restrict__ col_ind_out, + f_t* __restrict__ val_out) { - static_assert(std::is_signed_v); - - // Device CSC -> CSR: col_start[], i[], x[] (this) -> Arow.row_start[], j[], x[]. - // Nonzeros are reordered by sorting (row, col) so each CSR row segment is contiguous. - - i_t const nz = nz_max; - - Arow.m = m; - Arow.n = n; - Arow.nz_max = nz_max; - Arow.row_start.resize(m + 1, stream); - Arow.j.resize(nz, stream); - Arow.x.resize(nz, stream); + const i_t col = static_cast(blockIdx.x); + if (col >= n_cols) { return; } + const i_t col_end = col_start[col + 1]; + for (i_t p = col_start[col] + static_cast(threadIdx.x); p < col_end; + p += static_cast(blockDim.x)) { + const i_t q = atomicAdd(next_pos + row_ind[p], i_t(1)); + col_ind_out[q] = col; + val_out[q] = csc_val[p]; + } +} - auto exec = rmm::exec_policy(stream); +// Device CSC -> CSR on raw arrays. Doubles as a CSC transpose: CSR(A) and CSC(A^T) hold the +// same three arrays, so only the dimensions the caller records differ. +template +void csc_to_csr_on_device(i_t m, + i_t n, + i_t nz, + const i_t* col_start, + const i_t* row_ind, + const f_t* csc_val, + i_t* out_offsets, + i_t* out_indices, + f_t* out_values, + cuda::stream_ref stream) +{ + static_assert(std::is_signed_v); if (nz == 0) { - // Empty matrix: row_start all zero; j/x unused. - RAFT_CUDA_TRY(cudaMemsetAsync(Arow.row_start.data(), 0, sizeof(i_t) * (m + 1), stream.get())); + // Empty matrix: offsets all zero; indices/values unused. + RAFT_CUDA_TRY(cudaMemsetAsync(out_offsets, 0, sizeof(i_t) * (m + 1), stream.get())); return; } - // Per-row nnz from CSC row indices i[] (one atomic add per nonzero). + auto exec = rmm::exec_policy(stream); + + // Per-row nnz from the CSC row indices (one atomic add per nonzero). rmm::device_uvector row_counts(m, stream); RAFT_CUDA_TRY(cudaMemsetAsync(row_counts.data(), 0, sizeof(i_t) * m, stream.get())); thrust::for_each(exec, thrust::make_counting_iterator(0), thrust::make_counting_iterator(nz), - [row_ind = i.data(), counts = row_counts.data()] __device__(i_t p) { + [row_ind, counts = row_counts.data()] __device__(i_t p) { atomicAdd(counts + row_ind[p], i_t(1)); }); - // CSR row pointers: exclusive prefix sum of row_counts; Arow.row_start[m] = nz. + // Row pointers: exclusive prefix sum of row_counts; out_offsets[m] = nz. rmm::device_buffer scan_tmp; std::size_t scan_bytes = 0; cub::DeviceScan::ExclusiveSum( - nullptr, scan_bytes, row_counts.data(), Arow.row_start.data(), m, stream.get()); + nullptr, scan_bytes, row_counts.data(), out_offsets, m, stream.get()); scan_tmp.resize(scan_bytes, stream); cub::DeviceScan::ExclusiveSum( - scan_tmp.data(), scan_bytes, row_counts.data(), Arow.row_start.data(), m, stream.get()); - - RAFT_CUDA_TRY(cudaMemcpyAsync( - Arow.row_start.data() + m, &nz, sizeof(i_t), cudaMemcpyHostToDevice, stream.get())); - - // rows[]: CSC row indices (sort key). Arow.j / Arow.x hold (col, val) per flat CSC index, - // then sort_by_key permutes j and x in place into CSR (row, col) order. - rmm::device_uvector rows(nz, stream); - raft::copy(rows.data(), i.data(), nz, stream); - raft::copy(Arow.x.data(), x.data(), nz, stream); - - // Global CSC position p lies in column c iff col_start[c] <= p < col_start[c+1]. - thrust::tabulate(exec, - thrust::device_pointer_cast(Arow.j.data()), - thrust::device_pointer_cast(Arow.j.data() + nz), - [cs = col_start.data(), nn_c = n] __device__(i_t p) { - i_t lo = 0; - i_t hi = nn_c; - while (lo < hi) { - i_t mid = lo + (hi - lo) / 2; - if (cs[mid] <= p) { - lo = mid + 1; - } else { - hi = mid; - } - } - return lo - 1; - }); + scan_tmp.data(), scan_bytes, row_counts.data(), out_offsets, m, stream.get()); + + RAFT_CUDA_TRY( + cudaMemcpyAsync(out_offsets + m, &nz, sizeof(i_t), cudaMemcpyHostToDevice, stream.get())); + + // Scatter every nonzero into its row's segment. + rmm::device_uvector next_pos(m, stream); + raft::copy(next_pos.data(), out_offsets, m, stream); + + rmm::device_uvector indices_unsorted(nz, stream); + rmm::device_uvector values_unsorted(nz, stream); + constexpr int scatter_block_size = 256; + csc_to_csr_scatter_kernel + <<(n), scatter_block_size, 0, stream.get()>>>(n, + col_start, + row_ind, + csc_val, + next_pos.data(), + indices_unsorted.data(), + values_unsorted.data()); + RAFT_CUDA_TRY(cudaPeekAtLastError()); + + // Sort each segment by index; column ids are unique per row, so the result is deterministic. + rmm::device_buffer sort_tmp; + std::size_t sort_bytes = 0; + cub::DeviceSegmentedSort::SortPairs(nullptr, + sort_bytes, + indices_unsorted.data(), + out_indices, + values_unsorted.data(), + out_values, + nz, + m, + out_offsets, + out_offsets + 1, + stream.get()); + sort_tmp.resize(sort_bytes, stream); + cub::DeviceSegmentedSort::SortPairs(sort_tmp.data(), + sort_bytes, + indices_unsorted.data(), + out_indices, + values_unsorted.data(), + out_values, + nz, + m, + out_offsets, + out_offsets + 1, + stream.get()); +} - // CSR column order: sort (row, col) lexicographically; values follow the same permutation. - auto row_iter = thrust::device_pointer_cast(rows.data()); - auto col_iter = thrust::device_pointer_cast(Arow.j.data()); - thrust::sort_by_key(exec, - thrust::make_zip_iterator(thrust::make_tuple(row_iter, col_iter)), - thrust::make_zip_iterator(thrust::make_tuple(row_iter + nz, col_iter + nz)), - thrust::device_pointer_cast(Arow.x.data())); +template +void device_csc_matrix_t::to_compressed_row(device_csr_matrix_t& Arow, + cuda::stream_ref stream) const +{ + Arow.m = m; + Arow.n = n; + Arow.nz_max = nz_max; + Arow.row_start.resize(m + 1, stream); + Arow.j.resize(nz_max, stream); + Arow.x.resize(nz_max, stream); + + csc_to_csr_on_device(m, + n, + nz_max, + col_start.data(), + i.data(), + x.data(), + Arow.row_start.data(), + Arow.j.data(), + Arow.x.data(), + stream); +} + +template +void device_csc_matrix_t::transpose(device_csc_matrix_t& AT, + cuda::stream_ref stream) const +{ + // A^T is n x m, and its CSC arrays are exactly the CSR arrays of A. + AT.m = n; + AT.n = m; + AT.nz_max = nz_max; + AT.col_start.resize(m + 1, stream); + AT.i.resize(nz_max, stream); + AT.x.resize(nz_max, stream); + + csc_to_csr_on_device(m, + n, + nz_max, + col_start.data(), + i.data(), + x.data(), + AT.col_start.data(), + AT.i.data(), + AT.x.data(), + stream); } } // namespace cuopt::mathematical_optimization::barrier diff --git a/cpp/src/barrier/scaling_gpu.cu b/cpp/src/barrier/scaling_gpu.cu new file mode 100644 index 0000000000..57e76920a5 --- /dev/null +++ b/cpp/src/barrier/scaling_gpu.cu @@ -0,0 +1,465 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +// GPU port of Ruiz equilibration for the barrier/SOCP-QP path (see scaling.cpp's +// `scaling()` for the CPU reference implementation this mirrors step-for-step). Only +// the Ruiz-equilibration branch is implemented here; callers must only invoke this +// for problems with second-order cones or a quadratic objective (the same condition +// `scaling()` uses to select the Ruiz branch over the plain geometric-mean scaling). + +#include + +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace cuopt::mathematical_optimization::simplex { + +namespace { + +using cuopt::mathematical_optimization::barrier::device_csc_matrix_t; +using cuopt::mathematical_optimization::barrier::device_csr_matrix_t; + +// row_norm[i] = max_j |A(i,j)|, computed straight off CSC: A.i[p] is the row of nonzero p, +// so the per-row maxima need no row-contiguous (CSR) copy of the matrix. Mirrors +// compute_row_inf_norms in scaling.cpp. +template +void compute_row_inf_norms(const device_csc_matrix_t& A, + rmm::device_uvector& row_norm, + rmm::cuda_stream_view stream) +{ + row_norm.resize(A.m, stream); + thrust::fill(rmm::exec_policy(stream), row_norm.begin(), row_norm.end(), f_t(0)); + thrust::for_each(rmm::exec_policy(stream), + thrust::make_counting_iterator(i_t(0)), + thrust::make_counting_iterator(A.nz_max), + [x = A.x.data(), row = A.i.data(), rmax = row_norm.data()] __device__(i_t p) { + const f_t a = raft::abs(x[p]); + // Read first: most nonzeros lose to the running max, so the CAS loop + // inside myAtomicMax is usually skipped entirely. Collisions are rare + // anyway -- a CSC column cannot repeat a row, so consecutive threads + // within a column target distinct accumulators. + if (a > rmax[row[p]]) { raft::myAtomicMax(rmax + row[p], a); } + }); +} + +// max_i |transform(input[i])| over the whole array (single segment); used for the +// one-shot imbalance-ratio heuristic, not the per-iteration row/column reduces. +template +f_t whole_array_abs_max(InputIt input, size_t n, rmm::cuda_stream_view stream) +{ + if (n == 0) return f_t(0); + auto abs_it = thrust::make_transform_iterator(input, cuopt::abs_value_transform_t{}); + return thrust::reduce( + rmm::exec_policy(stream), abs_it, abs_it + n, f_t(0), cuopt::max_op_t{}); +} + +// min_i (|transform(input[i])| > 0 ? |value| : sentinel) over the whole array, sentinel == +// std::numeric_limits::max() when nothing is nonzero -- mirrors the host code's +// "ignore exact zeros" min-norm loop (scaling.cpp:45-46,53-56) bit for bit, including the +// degenerate all-zero case where the ratio check below reduces to `min_row_norm > 0`. +template +f_t whole_array_nonzero_abs_min(InputIt input, size_t n, rmm::cuda_stream_view stream) +{ + const f_t sentinel = std::numeric_limits::max(); + if (n == 0) return sentinel; + auto nz_it = thrust::make_transform_iterator(input, [sentinel] __device__(f_t value) -> f_t { + const f_t abs_value = raft::abs(value); + return abs_value > f_t(0) ? abs_value : sentinel; + }); + return thrust::reduce( + rmm::exec_policy(stream), nz_it, nz_it + n, sentinel, cuopt::min_op_t{}); +} + +// max/min per-segment |value| via cub::DeviceSegmentedReduce over an arbitrary offsets +// pair (row_start for CSR, col_start for CSC); segments with no nonzeros reduce to 0. +template +void segmented_abs_max(const f_t* values, + OffsetBeginIt begin_offsets, + OffsetEndIt end_offsets, + i_t num_segments, + f_t* out, + rmm::device_buffer& temp_storage, + rmm::cuda_stream_view stream) +{ + if (num_segments == 0) return; + auto abs_it = thrust::make_transform_iterator(values, cuopt::abs_value_transform_t{}); + size_t bytes = 0; + cub::DeviceSegmentedReduce::Reduce(nullptr, + bytes, + abs_it, + out, + num_segments, + begin_offsets, + end_offsets, + cuopt::max_op_t{}, + f_t(0), + stream); + temp_storage.resize(bytes, stream); + cub::DeviceSegmentedReduce::Reduce(temp_storage.data(), + bytes, + abs_it, + out, + num_segments, + begin_offsets, + end_offsets, + cuopt::max_op_t{}, + f_t(0), + stream); +} + +} // namespace + +template +i_t scaling_ruiz_gpu(const lp_problem_t& unscaled, + const simplex_solver_settings_t& settings, + lp_problem_t& scaled, + std::vector& column_scaling, + std::vector& row_scaling, + std::shared_ptr>& device_A, + std::shared_ptr>& device_Q) +{ + scaled = unscaled; + i_t m = scaled.num_rows; + i_t n = scaled.num_cols; + bool has_q = unscaled.Q.n > 0; + + rmm::cuda_stream_view stream = unscaled.handle_ptr->get_stream(); + + // Unconditional, so the early return below cannot leave a stale matrix in the caller's hands. + device_A.reset(); + device_Q.reset(); + row_scaling.assign(m, 1.0); + + // --- Upload only what the skip heuristic needs; the rest of the setup is deferred + // until after the decision, so a skipped problem pays almost nothing. --- + device_csc_matrix_t dA(scaled.A, stream); + device_csr_matrix_t dQ(scaled.Q, stream); + + // --- One-shot imbalance heuristic (mirrors scaling.cpp:43-116) --- + rmm::device_buffer scratch; + // Holds the raw row inf-norms, both for the heuristic here and in each Ruiz iteration + // below, where it is then converted in place into that iteration's row scale factors. + rmm::device_uvector r(0, stream); + compute_row_inf_norms(dA, r, stream); + f_t max_row_norm = whole_array_abs_max(r.data(), m, stream); + f_t min_row_norm = whole_array_nonzero_abs_min(r.data(), m, stream); + f_t row_norm_ratio = (min_row_norm > 0) ? max_row_norm / min_row_norm : f_t(1.0); + + rmm::device_uvector col_max_full(n, stream); + segmented_abs_max(dA.x.data(), + dA.col_start.data(), + dA.col_start.data() + 1, + n, + col_max_full.data(), + scratch, + stream); + f_t max_col_norm = whole_array_abs_max(col_max_full.data(), n, stream); + f_t min_col_norm = whole_array_nonzero_abs_min(col_max_full.data(), n, stream); + f_t col_norm_ratio = (min_col_norm > 0) ? max_col_norm / min_col_norm : f_t(1.0); + + f_t q_ratio = f_t(1.0); + if (has_q) { + f_t max_q = whole_array_abs_max(dQ.x.data(), dQ.nz_max, stream); + f_t min_q = whole_array_nonzero_abs_min(dQ.x.data(), dQ.nz_max, stream); + if (min_q <= max_q) { q_ratio = max_q / min_q; } + } + + const i_t ruiz_mode = settings.qcqp_ruiz_equilibration; + const bool balanced = row_norm_ratio < 100.0 && col_norm_ratio < 5e4 && q_ratio < 100.0; + const bool skip_ruiz = (ruiz_mode == 0) || (ruiz_mode < 0 && balanced); + + if (skip_ruiz) { + if (ruiz_mode == 0) { + settings.log.printf("Skipping Ruiz equilibration (qcqp_hyper_ruiz_equilibration = 0)\n"); + } else { + settings.log.printf( + "Skipping Ruiz equilibration (row norm ratio %.1f, column norm ratio %.1f < 5e4, Q coeff " + "ratio %.1f < 100)\n", + row_norm_ratio, + col_norm_ratio, + q_ratio); + } + column_scaling.assign(n, 1.0); + return 0; + } + if (ruiz_mode > 0) { + settings.log.printf( + "Applying Ruiz equilibration (qcqp_hyper_ruiz_equilibration = 1, row norm ratio %.1f, " + "column norm ratio %.1f, Q coeff ratio %.1f) [GPU]\n", + row_norm_ratio, + col_norm_ratio, + q_ratio); + } + + // --- Ruiz is actually going to run: upload the rest and build the index arrays. --- + std::vector col_scale_host(n, 1.0); + rmm::device_uvector d_rhs = cuopt::device_copy(scaled.rhs, stream); + rmm::device_uvector d_objective = cuopt::device_copy(scaled.objective, stream); + rmm::device_uvector d_lower = cuopt::device_copy(scaled.lower, stream); + rmm::device_uvector d_upper = cuopt::device_copy(scaled.upper, stream); + rmm::device_uvector d_row_scale = cuopt::device_copy(row_scaling, stream); + rmm::device_uvector d_col_scale = cuopt::device_copy(col_scale_host, stream); + + // Per-nonzero column ids for A, built once (sparsity pattern is fixed across iterations). + dA.form_col_index(stream); // dA.col_index[p] = column of nonzero p + + const i_t cone_start = unscaled.second_order_cone_dims.empty() ? n : unscaled.cone_var_start; + const i_t num_cones = static_cast(unscaled.second_order_cone_dims.size()); + // Column boundaries of each cone (in the global column index space) and, for every + // cone column, which cone it belongs to -- both built once, used every iteration. + std::vector cone_col_offsets_host(num_cones + 1, cone_start); + for (i_t k = 0; k < num_cones; ++k) { + cone_col_offsets_host[k + 1] = cone_col_offsets_host[k] + unscaled.second_order_cone_dims[k]; + } + std::vector col_cone_id_host(n - cone_start); + for (i_t k = 0; k < num_cones; ++k) { + for (i_t j = cone_col_offsets_host[k]; j < cone_col_offsets_host[k + 1]; ++j) { + col_cone_id_host[j - cone_start] = k; + } + } + rmm::device_uvector d_cone_col_offsets = cuopt::device_copy(cone_col_offsets_host, stream); + rmm::device_uvector d_col_cone_id = cuopt::device_copy(col_cone_id_host, stream); + + // --- Ruiz iteration loop (mirrors scaling.cpp:123-224) --- + constexpr i_t max_ruiz_iterations = 10; + rmm::device_uvector c(n, stream); + rmm::device_uvector col_max_linear(cone_start, stream); + rmm::device_uvector qrow_max_linear(has_q ? cone_start : 0, stream); + rmm::device_uvector cone_max(num_cones, stream); + + for (i_t iter = 0; iter < max_ruiz_iterations; ++iter) { + f_t max_deviation = 0.0; + + // --- Row scaling: r[i] = 1/sqrt(max_j |A(i,j)|) --- + // On the first pass r still holds the row inf-norms computed for the skip heuristic + // above, and A has not been touched since, so only recompute once A has been scaled. + if (iter > 0) { compute_row_inf_norms(dA, r, stream); } + max_deviation = std::max( + max_deviation, + whole_array_abs_max(thrust::make_transform_iterator( + r.data(), [] __device__(f_t v) -> f_t { return v - f_t(1); }), + m, + stream)); + thrust::transform( + rmm::exec_policy(stream), r.data(), r.data() + m, r.data(), [] __device__(f_t rm) { + return rm > 0 ? f_t(1) / std::sqrt(rm) : f_t(1); + }); + + thrust::for_each( + rmm::exec_policy(stream), + thrust::make_counting_iterator(i_t(0)), + thrust::make_counting_iterator(dA.nz_max), + [x = dA.x.data(), row = dA.i.data(), r = r.data()] __device__(i_t p) { x[p] *= r[row[p]]; }); + thrust::transform(rmm::exec_policy(stream), + d_rhs.data(), + d_rhs.data() + m, + r.data(), + d_rhs.data(), + cuda::std::multiplies{}); + thrust::transform(rmm::exec_policy(stream), + d_row_scale.data(), + d_row_scale.data() + m, + r.data(), + d_row_scale.data(), + cuda::std::multiplies{}); + + // --- Column scaling: linear columns [0, cone_start) combine A and Q; cone columns + // use one uniform scale per cone. --- + if (cone_start > 0) { + segmented_abs_max(dA.x.data(), + dA.col_start.data(), + dA.col_start.data() + 1, + cone_start, + col_max_linear.data(), + scratch, + stream); + if (has_q) { + segmented_abs_max(dQ.x.data(), + dQ.row_start.data(), + dQ.row_start.data() + 1, + cone_start, + qrow_max_linear.data(), + scratch, + stream); + thrust::transform(rmm::exec_policy(stream), + col_max_linear.data(), + col_max_linear.data() + cone_start, + qrow_max_linear.data(), + col_max_linear.data(), + cuopt::max_op_t{}); + } + max_deviation = + std::max(max_deviation, + whole_array_abs_max( + thrust::make_transform_iterator( + col_max_linear.data(), [] __device__(f_t v) -> f_t { return v - f_t(1); }), + cone_start, + stream)); + thrust::transform(rmm::exec_policy(stream), + col_max_linear.data(), + col_max_linear.data() + cone_start, + c.data(), + [] __device__(f_t cm) { return cm > 0 ? f_t(1) / std::sqrt(cm) : f_t(1); }); + } + if (num_cones > 0) { + auto begin_it = + thrust::make_permutation_iterator(dA.col_start.data(), d_cone_col_offsets.data()); + auto end_it = + thrust::make_permutation_iterator(dA.col_start.data(), d_cone_col_offsets.data() + 1); + segmented_abs_max( + dA.x.data(), begin_it, end_it, num_cones, cone_max.data(), scratch, stream); + max_deviation = + std::max(max_deviation, + whole_array_abs_max( + thrust::make_transform_iterator( + cone_max.data(), [] __device__(f_t v) -> f_t { return v - f_t(1); }), + num_cones, + stream)); + thrust::transform(rmm::exec_policy(stream), + cone_max.data(), + cone_max.data() + num_cones, + cone_max.data(), + [] __device__(f_t cm) { return cm > 0 ? f_t(1) / std::sqrt(cm) : f_t(1); }); + thrust::gather(rmm::exec_policy(stream), + d_col_cone_id.data(), + d_col_cone_id.data() + (n - cone_start), + cone_max.data(), + c.data() + cone_start); + } + + thrust::for_each(rmm::exec_policy(stream), + thrust::make_counting_iterator(i_t(0)), + thrust::make_counting_iterator(dA.nz_max), + [x = dA.x.data(), col = dA.col_index.data(), c = c.data()] __device__(i_t p) { + x[p] *= c[col[p]]; + }); + thrust::transform(rmm::exec_policy(stream), + d_objective.data(), + d_objective.data() + n, + c.data(), + d_objective.data(), + cuda::std::multiplies{}); + thrust::transform(rmm::exec_policy(stream), + d_col_scale.data(), + d_col_scale.data() + n, + c.data(), + d_col_scale.data(), + cuda::std::multiplies{}); + thrust::for_each( + rmm::exec_policy(stream), + thrust::make_counting_iterator(i_t(0)), + thrust::make_counting_iterator(n), + [lower = d_lower.data(), upper = d_upper.data(), c = c.data()] __device__(i_t j) { + if (lower[j] > f_t(-1e20)) lower[j] /= c[j]; + if (upper[j] < f_t(1e20)) upper[j] /= c[j]; + }); + if (has_q) { + // Row-parallel so the row index comes from the iteration variable, as in scaling.cpp. + // Deriving it per nonzero from row_start instead would mis-attribute every nonzero + // after an empty row, and Q has an empty row for every variable with no quadratic term. + thrust::for_each( + rmm::exec_policy(stream), + thrust::make_counting_iterator(i_t(0)), + thrust::make_counting_iterator(dQ.m), + [x = dQ.x.data(), rs = dQ.row_start.data(), col = dQ.j.data(), c = c.data()] __device__( + i_t i) { + for (i_t p = rs[i]; p < rs[i + 1]; ++p) { + x[p] *= c[i] * c[col[p]]; + } + }); + } + + if (max_deviation < 0.1) break; + } + + // --- Finalize: invert accumulated reciprocal scales (mirrors scaling.cpp:226-235) --- + thrust::transform(rmm::exec_policy(stream), + d_col_scale.data(), + d_col_scale.data() + n, + d_col_scale.data(), + [] __device__(f_t v) { return f_t(1) / v; }); + thrust::transform(rmm::exec_policy(stream), + d_row_scale.data(), + d_row_scale.data() + m, + d_row_scale.data(), + [] __device__(f_t v) { return f_t(1) / v; }); + + const f_t a_max = whole_array_abs_max(dA.x.data(), dA.nz_max, stream); + const f_t a_min = whole_array_nonzero_abs_min(dA.x.data(), dA.nz_max, stream); + + // --- Download scaled problem and scale vectors back to host --- + // SOCP goes straight to the barrier's augmented path, where no host code reads A's values, so + // A stays on device rather than being downloaded and immediately uploaded again. Ruiz rescales + // values but never changes the sparsity pattern, so scaled.A already carries the right + // col_start/i from the host copy above; only x differs, and clearing it keeps anyone from + // reading the stale unscaled values that would otherwise be left behind. + if (!unscaled.second_order_cone_dims.empty()) { + scaled.A.x.clear(); + scaled.A.x.shrink_to_fit(); + device_A = std::make_shared>(std::move(dA)); + } else { + scaled.A = dA.to_host(stream); + } + scaled.Q = dQ.to_host(stream); + // Q stays on host as well, so this is purely so the barrier need not upload it again. Q is + // symmetric, so its CSR arrays are also its CSC arrays and the handover is a relabel. + if (dQ.nz_max > 0) { + auto dQ_csc = std::make_shared>(stream); + dQ_csc->m = dQ.m; + dQ_csc->n = dQ.m; + dQ_csc->nz_max = dQ.nz_max; + dQ_csc->col_start = std::move(dQ.row_start); + dQ_csc->i = std::move(dQ.j); + dQ_csc->x = std::move(dQ.x); + device_Q = std::move(dQ_csc); + } + scaled.rhs = cuopt::host_copy(d_rhs, stream); + scaled.objective = cuopt::host_copy(d_objective, stream); + scaled.lower = cuopt::host_copy(d_lower, stream); + scaled.upper = cuopt::host_copy(d_upper, stream); + column_scaling = cuopt::host_copy(d_col_scale, stream); + row_scaling = cuopt::host_copy(d_row_scale, stream); + + settings.log.printf("Ruiz equilibration: coefficient range [%e, %e] [GPU]\n", a_min, a_max); + return 0; +} + +#ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE + +template int scaling_ruiz_gpu( + const lp_problem_t& unscaled, + const simplex_solver_settings_t& settings, + lp_problem_t& scaled, + std::vector& column_scaling, + std::vector& row_scaling, + std::shared_ptr>& device_A, + std::shared_ptr>& device_Q); + +#endif + +} // namespace cuopt::mathematical_optimization::simplex diff --git a/cpp/src/barrier/scaling_gpu.cuh b/cpp/src/barrier/scaling_gpu.cuh new file mode 100644 index 0000000000..0bdad182b7 --- /dev/null +++ b/cpp/src/barrier/scaling_gpu.cuh @@ -0,0 +1,39 @@ +/* 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 { +// Only ever held by shared_ptr below, so the definition (which needs nvcc) stays out of here. +template +class device_csc_matrix_t; +} // namespace cuopt::mathematical_optimization::barrier + +namespace cuopt::mathematical_optimization::simplex { + +// GPU-based Ruiz scaling. +// +// `device_A` receives the scaled A when it is left on device (second-order cones only, where +// `scaled.A` then keeps its sparsity pattern but has an empty `x` and this is the only copy of the +// values); `device_Q` receives the scaled Q. Both are null when there is nothing to hand over, and +// the barrier uploads from `scaled` instead. +template +i_t scaling_ruiz_gpu(const lp_problem_t& unscaled, + const simplex_solver_settings_t& settings, + lp_problem_t& scaled, + std::vector& column_scaling, + std::vector& row_scaling, + std::shared_ptr>& device_A, + std::shared_ptr>& device_Q); + +} // namespace cuopt::mathematical_optimization::simplex diff --git a/cpp/src/barrier/sparse_matrix_kernels.cuh b/cpp/src/barrier/sparse_matrix_kernels.cuh index c736e67aa4..3fbfa703ee 100644 --- a/cpp/src/barrier/sparse_matrix_kernels.cuh +++ b/cpp/src/barrier/sparse_matrix_kernels.cuh @@ -18,19 +18,25 @@ namespace cuopt::mathematical_optimization::barrier { template +// A is passed as its CSR pieces rather than as a matrix so callers can hand over a matrix that +// only happens to be CSR(A), such as a CSC holding A^T. void initialize_cusparse_data(raft::handle_t const* handle, - device_csr_matrix_t& A, + i_t A_rows, + i_t A_cols, + i_t A_nnz, + i_t* A_offsets, + i_t* A_indices, + f_t* A_values, device_csc_matrix_t& DAT, device_csr_matrix_t& ADAT, cusparse_info_t& cusparse_data) { - auto A_nnz = A.nz_max; auto DAT_nnz = DAT.nz_max; f_t chunk_fraction = 0.15; // Create matrix descriptors cusparse_data.matA_descr = - pdlp::make_csr(A.m, A.n, A_nnz, A.row_start.data(), A.j.data(), A.x.data()); + pdlp::make_csr(A_rows, A_cols, A_nnz, A_offsets, A_indices, A_values); cusparse_data.matDAT_descr = pdlp::make_csr( DAT.n, DAT.m, DAT_nnz, DAT.col_start.data(), DAT.i.data(), DAT.x.data()); cusparse_data.matADAT_descr = pdlp::make_csr( @@ -114,9 +120,9 @@ void initialize_cusparse_data(raft::handle_t const* handle, } template +// Operates purely on the descriptors built by initialize_cusparse_data; ADAT is resized to the +// nnz cuSPARSE reports. void multiply_kernels(raft::handle_t const* handle, - device_csr_matrix_t& A, - device_csc_matrix_t& DAT, device_csr_matrix_t& ADAT, cusparse_info_t& cusparse_data) { diff --git a/cpp/src/barrier/translate_soc.hpp b/cpp/src/barrier/translate_soc.hpp index f89994b949..e388cce026 100644 --- a/cpp/src/barrier/translate_soc.hpp +++ b/cpp/src/barrier/translate_soc.hpp @@ -69,6 +69,9 @@ void convert_quadratic_constraints_to_second_order_cones( // Use a practical tolerance for text-parsed MPS numeric values. const f_t tol = std::numeric_limits::epsilon() * 2; + // Rows appended below all land after these, so the model's own keep their indices. + user_problem.original_num_rows = csr_A.m; + // Derive implied lower bounds from singleton inequality rows. // Used to check if SOC head variables have implied non-negativity from the constraint system // without actually modifying the variable bounds (which would add barrier terms). @@ -838,6 +841,8 @@ void convert_quadratic_constraints_to_second_order_cones( const i_t m_old = csr_A.m; const i_t m_new = static_cast(m_old + cone_alias_pairs.size()); + user_problem.cone_variables_aliased = true; + user_problem.objective.resize(n_new, 0); user_problem.lower.resize(n_new, -std::numeric_limits::infinity()); user_problem.upper.resize(n_new, std::numeric_limits::infinity()); diff --git a/cpp/src/dual_simplex/scaling.cpp b/cpp/src/dual_simplex/scaling.cpp index 98c409a630..86db8512a7 100644 --- a/cpp/src/dual_simplex/scaling.cpp +++ b/cpp/src/dual_simplex/scaling.cpp @@ -12,6 +12,21 @@ namespace cuopt::mathematical_optimization::simplex { +namespace { + +// row_norm[i] = max_j |A(i,j)|, the infinity norm of row i of A. +template +void compute_row_inf_norms(const csc_matrix_t& A, std::vector& row_norm) +{ + row_norm.assign(A.m, 0.0); + const i_t nz = A.col_start[A.n]; + for (i_t p = 0; p < nz; ++p) { + row_norm[A.i[p]] = std::max(row_norm[A.i[p]], std::abs(A.x[p])); + } +} + +} // namespace + template i_t scaling(const lp_problem_t& unscaled, const simplex_solver_settings_t& settings, @@ -35,21 +50,18 @@ i_t scaling(const lp_problem_t& unscaled, if (!unscaled.second_order_cone_dims.empty() || unscaled.Q.n > 0) { // col_scale and row_scale accumulate reciprocal scale factors during Ruiz iterations. std::vector col_scale(n, 1.0); + // row inf-norms, used for both the skip heuristic and the Ruiz iterations. + std::vector r; // Decide whether Ruiz scaling is needed by checking row- and column-norm // imbalance. If both max_norm / min_norm ratios are small, the matrix is // already well-conditioned and scaling can hurt (e.g. by amplifying tiny // noise coefficients). - csr_matrix_t Arow_check(0, 0, 0); - scaled.A.to_compressed_row(Arow_check); + compute_row_inf_norms(scaled.A, r); f_t max_row_norm = 0; f_t min_row_norm = std::numeric_limits::max(); for (i_t i = 0; i < m; ++i) { - f_t row_norm = 0; - for (i_t p = Arow_check.row_start[i]; p < Arow_check.row_start[i + 1]; ++p) { - f_t a = std::abs(Arow_check.x[p]); - if (a > row_norm) row_norm = a; - } + const f_t row_norm = r[i]; if (row_norm > 0) { max_row_norm = std::max(max_row_norm, row_norm); min_row_norm = std::min(min_row_norm, row_norm); @@ -116,21 +128,16 @@ i_t scaling(const lp_problem_t& unscaled, } // Apply Ruiz equilibration - csr_matrix_t Arow(0, 0, 0); - scaled.A.to_compressed_row(Arow); - constexpr i_t max_ruiz_iterations = 10; for (i_t iter = 0; iter < max_ruiz_iterations; ++iter) { f_t max_deviation = 0.0; // --- Row scaling: scale each row by 1/sqrt(max|a_ij|) --- - std::vector r(m); + // On the first pass r still holds the row inf-norms computed for the skip heuristic + // above, and A has not been touched since, so only recompute once A has been scaled. + if (iter > 0) { compute_row_inf_norms(scaled.A, r); } for (i_t i = 0; i < m; ++i) { - f_t rm = 0.0; - for (i_t p = Arow.row_start[i]; p < Arow.row_start[i + 1]; ++p) { - f_t a = std::abs(Arow.x[p]); - if (a > rm) rm = a; - } + const f_t rm = r[i]; r[i] = rm > 0 ? 1.0 / std::sqrt(rm) : 1.0; max_deviation = std::max(max_deviation, std::abs(rm - 1.0)); } @@ -140,9 +147,6 @@ i_t scaling(const lp_problem_t& unscaled, } } for (i_t i = 0; i < m; ++i) { - for (i_t p = Arow.row_start[i]; p < Arow.row_start[i + 1]; ++p) { - Arow.x[p] *= r[i]; - } scaled.rhs[i] *= r[i]; row_scaling[i] *= r[i]; } @@ -196,11 +200,6 @@ i_t scaling(const lp_problem_t& unscaled, scaled.A.x[p] *= c[j]; } } - for (i_t i = 0; i < m; ++i) { - for (i_t p = Arow.row_start[i]; p < Arow.row_start[i + 1]; ++p) { - Arow.x[p] *= c[Arow.j[p]]; - } - } for (i_t j = 0; j < n; ++j) { scaled.objective[j] *= c[j]; col_scale[j] *= c[j]; diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 8b3eba56d3..df2ff95c09 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -85,6 +85,7 @@ struct simplex_solver_settings_t { postsolve_info(-1), barrier_presolve_bound_free_variables(-1), qcqp_ruiz_equilibration(-1), + gpu_ruiz_nnz_threshold(500000), barrier_initial_point_safeguard(10.0), check_Q(false), crossover(false), @@ -194,6 +195,8 @@ struct simplex_solver_settings_t { i_t postsolve_info; // -1 automatic (disabled), 0 disabled, 1 enabled i_t barrier_presolve_bound_free_variables; // -1 automatic, 0 disabled, 1 enabled i_t qcqp_ruiz_equilibration; // -1 automatic (imbalance heuristic), 0 disabled, 1 enabled + i_t gpu_ruiz_nnz_threshold; // nnz(A)+nnz(Q) above which the barrier path's Ruiz + // equilibration runs on GPU instead of CPU. f_t barrier_initial_point_safeguard; // margin pushing the barrier initial iterate into // the interior of the nonnegative orthant / SOC bool check_Q; // true to check if Q is positive semidefinite diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 6f6ecc88f9..3a92254b87 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -8,6 +8,7 @@ #include #include +#include #include @@ -36,6 +37,7 @@ #include #include #include +#include #include namespace cuopt::mathematical_optimization::simplex { @@ -47,6 +49,7 @@ void unscale_uncrush_barrier_to_user(const user_problem_t& user_proble const raft::handle_t* handle_ptr, i_t original_num_rows, i_t original_num_cols, + i_t converted_cone_var_start, const lp_problem_t& barrier_lp, const presolve_info_t& presolve_info, const std::vector& column_scales, @@ -68,7 +71,10 @@ void unscale_uncrush_barrier_to_user(const user_problem_t& user_proble unscaled_z); // Dummy converted LP: sizes only. Bound-free=0 so uncrush_solution never reads A. + // cone_var_start is the exception: uncrush_primal_solution needs the real value to undo the + // shift convert applied to the cone block. lp_problem_t converted(handle_ptr, original_num_rows, original_num_cols, 0); + converted.cone_var_start = converted_cone_var_start; lp_solution_t lp_solution(original_num_rows, original_num_cols); uncrush_solution(presolve_info, barrier_settings, @@ -91,6 +97,22 @@ void unscale_uncrush_barrier_to_user(const user_problem_t& user_proble solution.iterations = barrier_solution.iterations; } +// Presolve and scaling offset the data by a constant the maps alone cannot recover, so record +// what the barrier values are worth beyond crush(user values); updates re-add it. +template +std::vector shift_from(const std::vector& barrier_values, + const std::vector& crushed) +{ + if (crushed.size() != barrier_values.size()) { + throw std::runtime_error("crushed length disagrees with the cached barrier LP"); + } + std::vector shift(barrier_values.size()); + for (std::size_t k = 0; k < shift.size(); ++k) { + shift[k] = static_cast(barrier_values[k]) - crushed[k]; + } + return shift; +} + template void write_matlab(const std::string& filename, const simplex::lp_problem_t& lp) { @@ -419,15 +441,23 @@ lp_status_t solve_linear_program_with_barrier( 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() && + auto const* xf = (cache != nullptr && cache->dirty()) ? cache->transform() : nullptr; + const bool reuse_cached_data = + xf != nullptr && xf->barrier_lp != nullptr && + // Only quadratic-objective and cone models reach the barrier at all. + (!user_problem.Q_values.empty() || !user_problem.second_order_cone_dims.empty()) && + cuopt::mathematical_optimization::cone_layout_matches(*xf, user_problem) && + // run_barrier already resolved -1 to 0. The second check covers caches built by an earlier + // solve that did bound free variables, whose presolve state the reuse path cannot replay. settings.barrier_presolve_bound_free_variables == 0 && + xf->presolve_info.bounded_free_variables.empty() && user_problem.num_cols == xf->user_num_cols && user_problem.num_rows == xf->user_num_rows; - if (reuse_c_only) { + if (reuse_cached_data) { + if (cache->rhs_infeasible()) { + settings.log.printf("Barrier: update_rhs made an empty constraint row infeasible\n"); + return lp_status_t::INFEASIBLE; + } 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( @@ -439,6 +469,7 @@ lp_status_t solve_linear_program_with_barrier( cache->handle_ptr(), xf->original_num_rows, xf->original_num_cols, + xf->converted_cone_var_start, *xf->barrier_lp, xf->presolve_info, xf->column_scales, @@ -446,7 +477,7 @@ lp_status_t solve_linear_program_with_barrier( barrier_settings, barrier_solution, solution); - cache->set_c_dirty(false); + cache->mark_clean(); } else { cache->clear(); } @@ -480,7 +511,19 @@ lp_status_t solve_linear_program_with_barrier( 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); + // Scaled A/Q that the GPU scaling may leave on device for the barrier to adopt; null otherwise. + std::shared_ptr> device_A; + std::shared_ptr> device_Q; + const bool is_ruiz_candidate = + !presolved_lp.second_order_cone_dims.empty() || presolved_lp.Q.n > 0; + const i_t presolved_nnz = presolved_lp.A.col_start[presolved_lp.num_cols] + + (presolved_lp.Q.n > 0 ? presolved_lp.Q.row_start[presolved_lp.Q.m] : 0); + if (is_ruiz_candidate && presolved_nnz >= barrier_settings.gpu_ruiz_nnz_threshold) { + scaling_ruiz_gpu( + presolved_lp, barrier_settings, barrier_lp, column_scales, row_scales, device_A, device_Q); + } else { + 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); @@ -498,35 +541,64 @@ lp_status_t solve_linear_program_with_barrier( 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->pre_expansion_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); - solver_lp = xf->barrier_lp.get(); + xf->pre_expansion_num_rows = user_problem.original_num_rows; + xf->converted_cone_var_start = original_lp.cone_var_start; + xf->cone_head_bounds = cuopt::mathematical_optimization::record_cone_head_bounds(user_problem); + // Rows the expansion appended past the model's own; none when no expansion ran. + if (user_problem.original_num_rows > 0) { + xf->cone_row_rhs.assign(user_problem.rhs.begin() + user_problem.original_num_rows, + user_problem.rhs.end()); + } + xf->presolve_info = presolve_info; + xf->column_scales = column_scales; + xf->row_scales = row_scales; + xf->primal_tol = static_cast(barrier_settings.primal_tol); + // convert_range_rows zeroes rhs[i] onto the slack bounds and folding aggregates rows, so + // neither leaves the user RHS in barrier_lp->rhs. Plain inequality/equality slacks do. + // Aliased cone variables hide which model variable a head stands for, so the head bounds + // cannot be re-proved against a new RHS. + xf->rhs_update_supported = user_problem.num_range_rows == 0 && + !presolve_info.folding_info.is_folded && + !user_problem.cone_variables_aliased; + xf->barrier_lp = std::make_unique>(barrier_lp); + solver_lp = xf->barrier_lp.get(); cache->store_transform(std::move(xf)); } - barrier::barrier_solver_t barrier_solver(*solver_lp, presolve_info, barrier_settings); + barrier::barrier_solver_t barrier_solver( + *solver_lp, presolve_info, barrier_settings, std::move(device_A), std::move(device_Q)); 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 = cache->transform(); + // Crushing this solve's own data also checks the maps still describe it: presolve may + // have dualized an LP, in which case the crush throws. + // Both crushes take model coordinates. The expansion only appends rows, so the RHS prefix + // is already model-sized, but it permutes columns, so gather the objective back. + auto const model_objective = + cuopt::mathematical_optimization::gather_model_objective(*xf, user_problem.objective); + const int model_m = cuopt::mathematical_optimization::model_num_rows(*xf); try { - auto crushed = cuopt::mathematical_optimization::crush_user_linear_objective( - *xf, user_problem.objective.data(), user_problem.num_cols); - xf->linear_obj_shift.resize(static_cast(solver_lp->num_cols), 0.0); - if (static_cast(crushed.size()) == solver_lp->num_cols) { - for (int j = 0; j < solver_lp->num_cols; ++j) { - xf->linear_obj_shift[static_cast(j)] = - solver_lp->objective[static_cast(j)] - - crushed[static_cast(j)]; - } - } + auto const crushed = cuopt::mathematical_optimization::crush_user_linear_objective( + *xf, model_objective.data(), static_cast(model_objective.size())); + xf->linear_obj_shift = shift_from(solver_lp->objective, crushed); } catch (std::exception const&) { - xf->linear_obj_shift.assign(static_cast(solver_lp->num_cols), 0.0); + // A zero shift still lets an update run; it just contributes nothing. + xf->linear_obj_shift.assign(solver_lp->objective.size(), 0.0); + } + if (xf->rhs_update_supported) { + try { + auto const crushed = + cuopt::mathematical_optimization::crush_user_rhs(*xf, user_problem.rhs.data(), model_m); + xf->rhs_shift = shift_from(solver_lp->rhs, crushed); + } catch (std::exception const&) { + // Maps cannot reproduce this RHS, so refuse later updates. + xf->rhs_update_supported = false; + xf->rhs_shift.clear(); + } } } else { cache->clear(); diff --git a/cpp/src/dual_simplex/user_problem.hpp b/cpp/src/dual_simplex/user_problem.hpp index 419b33e1e1..879856fabe 100644 --- a/cpp/src/dual_simplex/user_problem.hpp +++ b/cpp/src/dual_simplex/user_problem.hpp @@ -69,6 +69,12 @@ struct user_problem_t { // expanded layout (num_cols) and must be projected back via original_col_to_expanded_col. i_t original_num_cols{0}; std::vector original_col_to_expanded_col; + // Row count before QCMATRIX->SOC expansion. The expansion only appends rows, so rows + // [0, original_num_rows) still hold the model's own constraints at their original indices. + i_t original_num_rows{0}; + // Set when a variable shared by several cones was given an alias column. The cone head a + // cache sees is then not the variable whose bounds the expansion checked. + bool cone_variables_aliased{false}; }; } // namespace cuopt::mathematical_optimization::simplex diff --git a/cpp/src/mip_heuristics/mip_scaling_strategy.cu b/cpp/src/mip_heuristics/mip_scaling_strategy.cu index 41525a9c65..67537ea242 100644 --- a/cpp/src/mip_heuristics/mip_scaling_strategy.cu +++ b/cpp/src/mip_heuristics/mip_scaling_strategy.cu @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -68,9 +69,13 @@ constexpr double big_m_abs_threshold = 1.0e4; constexpr double big_m_ratio_threshold = 1.0e4; template -struct abs_value_transform_t { - __device__ f_t operator()(f_t value) const { return raft::abs(value); } -}; +using abs_value_transform_t = cuopt::abs_value_transform_t; + +template +using max_op_t = cuopt::max_op_t; + +template +using min_op_t = cuopt::min_op_t; template struct nonzero_abs_or_inf_transform_t { @@ -86,22 +91,6 @@ struct nonzero_count_transform_t { __device__ i_t operator()(f_t value) const { return raft::abs(value) > f_t(0) ? i_t(1) : i_t(0); } }; -template -struct max_op_t { - __host__ __device__ item_t operator()(const item_t& lhs, const item_t& rhs) const - { - return lhs > rhs ? lhs : rhs; - } -}; - -template -struct min_op_t { - __host__ __device__ item_t operator()(const item_t& lhs, const item_t& rhs) const - { - return lhs < rhs ? lhs : rhs; - } -}; - struct gcd_op_t { __host__ __device__ std::int64_t operator()(std::int64_t lhs, std::int64_t rhs) const { diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 38202c6b51..7ab036a880 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -87,7 +87,7 @@ simplex::user_problem_t user_problem_from_transform( 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.objective = scatter_model_objective(xf, 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); @@ -96,8 +96,9 @@ simplex::user_problem_t user_problem_from_transform( 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_num_cols = xf.pre_expansion_num_cols; user_problem.original_col_to_expanded_col = xf.original_col_to_expanded_col; + user_problem.original_num_rows = xf.pre_expansion_num_rows; return user_problem; } @@ -519,6 +520,16 @@ optimization_problem_solution_t convert_dual_simplex_sol( method); } +// Bounding free variables writes presolve state the reuse path cannot replay, so the automatic +// (-1) choice resolves to 0 for a sequence solve. An explicit 1 is honored and forgoes reuse. +template +i_t effective_bound_free_variables(pdlp_solver_settings_t const& settings) +{ + return (settings.sequence_solve && settings.barrier_presolve_bound_free_variables < 0) + ? 0 + : settings.barrier_presolve_bound_free_variables; +} + template std::tuple, simplex::lp_status_t, f_t, f_t, f_t> run_barrier( const simplex::user_problem_t& user_problem, @@ -531,31 +542,30 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t f_t norm_rhs = vector_norm2(user_problem.rhs); simplex::simplex_solver_settings_t barrier_settings; - barrier_settings.num_gpus = settings.num_gpus; - barrier_settings.time_limit = settings.time_limit; - barrier_settings.iteration_limit = settings.iteration_limit; - barrier_settings.concurrent_halt = settings.concurrent_halt; - barrier_settings.folding = settings.folding; - barrier_settings.augmented = settings.augmented; - barrier_settings.dualize = settings.dualize; - barrier_settings.ordering = settings.ordering; - barrier_settings.barrier_dual_initial_point = settings.barrier_dual_initial_point; - barrier_settings.postsolve_info = settings.postsolve_info; - barrier_settings.barrier_presolve_bound_free_variables = - settings.barrier_presolve_bound_free_variables; - barrier_settings.barrier_initial_point_safeguard = settings.barrier_initial_point_safeguard; - barrier_settings.barrier = true; - barrier_settings.barrier_presolve = true; - barrier_settings.crossover = settings.crossover; - barrier_settings.eliminate_dense_columns = settings.eliminate_dense_columns; - barrier_settings.barrier_iterative_refinement = settings.barrier_iterative_refinement; - barrier_settings.barrier_adaptive_regularization = settings.barrier_adaptive_regularization; - barrier_settings.barrier_primal_regularization = settings.barrier_primal_regularization; - barrier_settings.barrier_dual_regularization = settings.barrier_dual_regularization; - barrier_settings.barrier_soc_threshold = settings.barrier_soc_threshold; - barrier_settings.barrier_step_scale = settings.barrier_step_scale; - barrier_settings.qcqp_ruiz_equilibration = settings.qcqp_ruiz_equilibration; - barrier_settings.cudss_deterministic = settings.cudss_deterministic; + barrier_settings.num_gpus = settings.num_gpus; + barrier_settings.time_limit = settings.time_limit; + barrier_settings.iteration_limit = settings.iteration_limit; + barrier_settings.concurrent_halt = settings.concurrent_halt; + barrier_settings.folding = settings.folding; + barrier_settings.augmented = settings.augmented; + barrier_settings.dualize = settings.dualize; + barrier_settings.ordering = settings.ordering; + barrier_settings.barrier_dual_initial_point = settings.barrier_dual_initial_point; + barrier_settings.postsolve_info = settings.postsolve_info; + barrier_settings.barrier_presolve_bound_free_variables = effective_bound_free_variables(settings); + barrier_settings.barrier_initial_point_safeguard = settings.barrier_initial_point_safeguard; + barrier_settings.barrier = true; + barrier_settings.barrier_presolve = true; + barrier_settings.crossover = settings.crossover; + barrier_settings.eliminate_dense_columns = settings.eliminate_dense_columns; + barrier_settings.barrier_iterative_refinement = settings.barrier_iterative_refinement; + barrier_settings.barrier_adaptive_regularization = settings.barrier_adaptive_regularization; + barrier_settings.barrier_primal_regularization = settings.barrier_primal_regularization; + barrier_settings.barrier_dual_regularization = settings.barrier_dual_regularization; + barrier_settings.barrier_soc_threshold = settings.barrier_soc_threshold; + barrier_settings.barrier_step_scale = settings.barrier_step_scale; + barrier_settings.qcqp_ruiz_equilibration = settings.qcqp_ruiz_equilibration; + barrier_settings.cudss_deterministic = settings.cudss_deterministic; barrier_settings.barrier_relaxed_feasibility_tol = settings.tolerances.relative_primal_tolerance; barrier_settings.barrier_relaxed_optimality_tol = settings.tolerances.relative_dual_tolerance; barrier_settings.barrier_relaxed_complementarity_tol = settings.tolerances.relative_gap_tolerance; @@ -1887,14 +1897,21 @@ optimization_problem_solution_t solve_qcqp( auto qcqp_timer = cuopt::timer_t(settings.time_limit); auto* cache = settings.barrier_cache; - auto const* xf = (cache != nullptr && cache->c_dirty()) ? cache->transform() : nullptr; + auto const* xf = (cache != nullptr && cache->dirty()) ? cache->transform() : nullptr; + // Must stay in lockstep with the gate in solve_linear_program_with_barrier: this path swaps + // in the slim user_problem_from_transform, so disagreement runs presolve on a fabricated + // problem. + // Cone models are compared in model coordinates: the cached counts are post-expansion. 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() && + effective_bound_free_variables(settings) == 0 && + xf->presolve_info.bounded_free_variables.empty() && + (op_problem.has_quadratic_objective() || op_problem.has_quadratic_constraints()) && + static_cast(op_problem.get_quadratic_constraints().size()) == + xf->num_quadratic_constraints && 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; + op_problem.get_n_variables() == model_num_cols(*xf) && + op_problem.get_n_constraints() == model_num_rows(*xf); if (problem_checking && !reuse_from_cache) { problem_checking_t::check_problem_representation(op_problem); @@ -1937,8 +1954,15 @@ optimization_problem_solution_t solve_qcqp( qcqp_timer, op_problem.get_handle_ptr(), settings.barrier_cache); + // A full solve creates the transform inside run_barrier. Record the model sense afterward, + // before a later update_linear_objective needs to map raw user c into barrier minimization + // space. Reuse keeps the sense of the workspace (and Q) built by that full solve. if (!reuse_from_cache && cache != nullptr && cache->transform() != nullptr) { cache->transform()->maximize = op_problem.get_sense(); + // Reuse never re-runs the cone expansion, so the gate above rejects a model that has + // gained or lost a quadratic constraint since the cache was built. + cache->transform()->num_quadratic_constraints = + static_cast(op_problem.get_quadratic_constraints().size()); } auto solution = convert_dual_simplex_sol(op_problem, std::get<0>(sol_dual_simplex), diff --git a/cpp/src/utilities/reduce_ops.cuh b/cpp/src/utilities/reduce_ops.cuh new file mode 100644 index 0000000000..e39d71d5a2 --- /dev/null +++ b/cpp/src/utilities/reduce_ops.cuh @@ -0,0 +1,36 @@ +/* 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 + +namespace cuopt { + +// abs(value); used as a cub/thrust transform-iterator input op ahead of an inf-norm reduce. +template +struct abs_value_transform_t { + __device__ f_t operator()(f_t value) const { return raft::abs(value); } +}; + +template +struct max_op_t { + __host__ __device__ item_t operator()(const item_t& lhs, const item_t& rhs) const + { + return lhs > rhs ? lhs : rhs; + } +}; + +template +struct min_op_t { + __host__ __device__ item_t operator()(const item_t& lhs, const item_t& rhs) const + { + return lhs < rhs ? lhs : rhs; + } +}; + +} // namespace cuopt diff --git a/cpp/tests/dual_simplex/unit_tests/device_sparse_matrix_test.cu b/cpp/tests/dual_simplex/unit_tests/device_sparse_matrix_test.cu new file mode 100644 index 0000000000..c76865e8f8 --- /dev/null +++ b/cpp/tests/dual_simplex/unit_tests/device_sparse_matrix_test.cu @@ -0,0 +1,157 @@ +/* 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::mathematical_optimization::barrier::test { + +namespace { + +// A distinct value encoding the entry's position, so a mis-permutation names the entry it came +// from in the failure output. Integral, so host and device copies compare bit-exact. +double value_at(int row, int col) +{ + const double magnitude = 1000.0 * (row + 1) + (col + 1); + return col % 2 == 0 ? magnitude : -magnitude; +} + +// Host CSC with the given row indices per column, in the order given (deliberately unsorted). +csc_matrix_t make_csc(int m, const std::vector>& rows_per_col) +{ + const int n = static_cast(rows_per_col.size()); + int nz = 0; + for (const auto& rows : rows_per_col) { + nz += static_cast(rows.size()); + } + + csc_matrix_t A(m, n, nz); + int p = 0; + A.col_start[0] = 0; + for (int j = 0; j < n; ++j) { + for (int r : rows_per_col[j]) { + A.i[p] = r; + A.x[p] = value_at(r, j); + ++p; + } + A.col_start[j + 1] = p; + } + return A; +} + +// The device conversion must reproduce the host reference exactly. +void expect_device_matches_host(const csc_matrix_t& A) +{ + auto stream = cuda::stream_ref{cudaStream_t{cudaStreamDefault}}; + + csr_matrix_t expected(A.m, A.n, A.col_start[A.n]); + A.to_compressed_row(expected); + + device_csc_matrix_t d_A(A, stream); + device_csr_matrix_t d_Arow(stream); + d_A.to_compressed_row(d_Arow, stream); + auto got = d_Arow.to_host(stream); + + ASSERT_EQ(got.m, expected.m); + ASSERT_EQ(got.n, expected.n); + EXPECT_EQ(got.row_start, expected.row_start); + EXPECT_EQ(got.j, expected.j); + EXPECT_EQ(got.x, expected.x); + + // The transpose shares the conversion, and CSC(A^T) holds the same arrays as CSR(A). + csc_matrix_t expected_t(1, 1, 1); + A.transpose(expected_t); + + device_csc_matrix_t d_AT(stream); + d_A.transpose(d_AT, stream); + auto got_t = d_AT.to_host(stream); + + ASSERT_EQ(got_t.m, expected_t.m); + ASSERT_EQ(got_t.n, expected_t.n); + EXPECT_EQ(got_t.col_start, expected_t.col_start); + EXPECT_EQ(got_t.i, expected_t.i); + EXPECT_EQ(got_t.x, expected_t.x); +} + +} // namespace + +TEST(device_sparse_matrix, csc_to_csr_empty_rows_and_columns) +{ + // 9 x 6 sparsity pattern: + // c0 c1 c2 c3 c4 c5 + // r0 . . . . . . + // r1 x . . x . x + // r2 . . . . . . + // r3 . . x . . x + // r4 . . . . . . + // r5 x . . . . x + // r6 . . . . . . + // r7 x . x . . x + // r8 . . . . . . + // + // Rows 0, 2, 4, 6 and 8 hold no entries, giving leading, interior and trailing empty CSR rows, + // i.e. zero-length sort segments. Columns 1 and 4 are empty, so their scatter blocks do no work. + // Every row indices list is out of order, and rows 1, 3, 5 and 7 each hold more than one entry, + // so the segmented sort has to restore column order rather than inherit it from the input. + const std::vector> rows_per_col = { + {5, 1, 7}, // c0 + {}, // c1, empty + {7, 3}, // c2 + {1}, // c3 + {}, // c4, empty + {3, 7, 5, 1}, // c5 + }; + + expect_device_matches_host(make_csc(9, rows_per_col)); +} + +TEST(device_sparse_matrix, csc_to_csr_dense_column) +{ + constexpr int m = 1000; + + // Column 2 holds every row, in descending order. The scatter kernel gives each column one + // 256-thread block, so this column makes its strided loop wrap several times. + std::vector all_rows_descending; + all_rows_descending.reserve(m); + for (int r = m - 1; r >= 0; --r) { + all_rows_descending.push_back(r); + } + + // The short columns repeat rows 0 and 999 so the first and last CSR rows hold several entries + // rather than just the one the dense column contributes. + const std::vector> rows_per_col = { + {900, 4, 500}, // c0 + {0, 999}, // c1 + all_rows_descending, // c2 + {999, 0, 7}, // c3 + {251, 250}, // c4 + }; + + expect_device_matches_host(make_csc(m, rows_per_col)); +} + +TEST(device_sparse_matrix, csc_to_csr_empty_matrix) +{ + // No nonzeros at all: the conversion takes its early return, which zeroes the offsets and + // launches no kernel. + expect_device_matches_host(make_csc(4, {{}, {}, {}})); +} + +TEST(device_sparse_matrix, csc_to_csr_single_entry) +{ + expect_device_matches_host(make_csc(1, {{0}})); +} + +} // namespace cuopt::mathematical_optimization::barrier::test diff --git a/cpp/tests/internal/CMakeLists.txt b/cpp/tests/internal/CMakeLists.txt index 9beef78df3..8375376344 100644 --- a/cpp/tests/internal/CMakeLists.txt +++ b/cpp/tests/internal/CMakeLists.txt @@ -10,6 +10,7 @@ ConfigureTest(NUMOPT_INTERNAL_TEST # dual_simplex ${CUOPT_TEST_DIR}/dual_simplex/unit_tests/solve.cpp ${CUOPT_TEST_DIR}/dual_simplex/unit_tests/solve_barrier.cu + ${CUOPT_TEST_DIR}/dual_simplex/unit_tests/device_sparse_matrix_test.cu ${CUOPT_TEST_DIR}/dual_simplex/unit_tests/right_looking_ldlt.cpp # linear_programming ${CUOPT_TEST_DIR}/linear_programming/pdlp_test.cu 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 4358471c37..03cc47f87b 100644 --- a/python/cuopt/cuopt/linear_programming/data_model/data_model.py +++ b/python/cuopt/cuopt/linear_programming/data_model/data_model.py @@ -231,12 +231,12 @@ def set_objective_coefficients(self, c): @catch_cuopt_exception def update_linear_objective(self, coefficients): """ - Cache reuse is QP-only: quadratic constraints take a full solve. - Update the linear objective coefficients for a sequence re-solve. + Writes ``coefficients`` onto this DataModel. If a barrier cache is present, also maps them into the cached barrier workspace and marks - it dirty (quadratic ``Q``, ``A``, and bounds must stay unchanged). + it dirty (quadratic ``Q``, ``A``, bounds, and the quadratic + constraints must stay unchanged). Parameters ---------- @@ -246,6 +246,29 @@ def update_linear_objective(self, coefficients): """ super().update_linear_objective(coefficients) + @catch_cuopt_exception + def update_rhs(self, b): + """ + Update the constraint right-hand sides (b) for a sequence re-solve. + + Writes ``b`` onto this DataModel. If a barrier cache is present, also + maps ``b`` into the cached barrier workspace and marks it dirty + (quadratic ``Q``, ``A``, row senses, bounds, and the quadratic + constraints must stay unchanged). + + Range rows and folding in the first solve are not supported and raise; + run a full solve for those models. Rows that presolve dropped as empty + are allowed: if the new ``b`` makes one infeasible, the next solve + reports infeasible without rerunning the interior point method. + + Parameters + ---------- + b : array-like of float64 + Constraint right-hand sides, length equal to the number of + constraints on the first ``sequence_solve``. + """ + super().update_rhs(b) + @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.pyx b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx index 017a81eeb7..3cb37f33b0 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 @@ -28,6 +28,7 @@ cdef extern from "Python.h": cdef extern from "cuopt/mathematical_optimization/utilities/barrier_cache.hpp" namespace "cuopt::mathematical_optimization": # noqa cdef cppclass barrier_cache_t: void update_linear_objective(const double* c, int n) except + + void update_rhs(const double* b, int m) except + def type_cast(np_obj, np_type, name): @@ -195,6 +196,32 @@ cdef class DataModel: cache.update_linear_objective(&c_view[0], c_view.shape[0]) self.c = new_c + def update_rhs(self, b): + """Update constraint right-hand sides (user-space ``b``). + + Always writes the DataModel RHS, and additionally crushes ``b`` into + the barrier cache when this model owns one. Crush runs first so a + length error leaves the DataModel RHS unchanged. + """ + cdef barrier_cache_t* cache + cdef double[::1] b_view + new_b = type_cast(b, np.float64, "b") + 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", + ) + b_view = np.ascontiguousarray(new_b, dtype=np.float64) + if b_view.shape[0] == 0: + cache.update_rhs(NULL, 0) + else: + cache.update_rhs(&b_view[0], b_view.shape[0]) + self.b = new_b + def set_objective_scaling_factor(self, objective_scaling_factor): self.objective_scaling_factor = objective_scaling_factor diff --git a/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py b/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py new file mode 100644 index 0000000000..9b1c0421e3 --- /dev/null +++ b/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py @@ -0,0 +1,513 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 + +"""Barrier cache reuse (``sequence_solve``) for the DataModel update APIs. + +Every re-solve through the cache is compared against a fresh full solve of the +same model, so no assertion depends on a hand-derived optimum. The models are +picked to exercise the parts of the RHS crush that a one-row QP leaves as +no-ops: mixed row senses, presolve-dropped empty rows, non-unit row scaling, +and the shift that comes from translating nonzero variable lower bounds. + +Each test also asserts the reuse path actually ran. Without that, a test passes +just as happily when the gate rejects the model and the solver quietly falls +back to a full solve, which returns the same answer. +""" + +import numpy as np +import pytest + +from cuopt.linear_programming import ( + data_model, + solver, + solver_settings, +) + +REUSE_LOG = "reusing cache" +RHS_INFEASIBLE_LOG = "update_rhs made an empty constraint row infeasible" +IPM_LOG = "Optimal solution found" + + +def _sequence_settings(): + settings = solver_settings.SolverSettings() + settings.set_parameter("sequence_solve", True) + # barrier_presolve_bound_free_variables stays at its -1 default on purpose: + # sequence_solve resolves automatic to 0, which the reuse assertions cover. + return settings + + +def _build( + values, indices, offsets, rhs, senses, lower, upper, objective=None +): + """A QP with quadratic term x^T x, so the model is barrier-eligible.""" + n = len(lower) + model = data_model.DataModel() + model.set_csr_constraint_matrix( + np.asarray(values, dtype=np.float64), + np.asarray(indices, dtype=np.int32), + np.asarray(offsets, dtype=np.int32), + ) + model.set_constraint_bounds(np.asarray(rhs, dtype=np.float64)) + model.set_row_types(senses) + model.set_objective_coefficients( + np.zeros(n) + if objective is None + else np.asarray(objective, dtype=np.float64) + ) + model.set_quadratic_objective_matrix( + np.ones(n), + np.arange(n, dtype=np.int32), + np.arange(n + 1, dtype=np.int32), + ) + model.set_variable_lower_bounds(np.asarray(lower, dtype=np.float64)) + model.set_variable_upper_bounds(np.asarray(upper, dtype=np.float64)) + return model + + +def _solve(model, settings, capfd): + """Solve and return the log this solve alone produced.""" + capfd.readouterr() + solution = solver.Solve(model, settings) + return solution, capfd.readouterr().out + + +def _full_solve(model_args, rhs, **overrides): + """Oracle: a fresh model with default settings, no cache in play.""" + args = dict(model_args, rhs=rhs, **overrides) + return solver.Solve(_build(**args), solver_settings.SolverSettings()) + + +def _assert_matches_oracle(reused, oracle): + assert reused.get_termination_reason() == "Optimal" + assert reused.get_termination_reason() == oracle.get_termination_reason() + assert reused.get_primal_objective() == pytest.approx( + oracle.get_primal_objective(), rel=1e-6, abs=1e-6 + ) + np.testing.assert_allclose( + np.asarray(reused.get_primal_solution()), + np.asarray(oracle.get_primal_solution()), + rtol=1e-5, + atol=1e-5, + ) + + +# x0 + x1 == b0 ; x1 + x2 <= b1 ; x0 + x2 >= b2. The equality and the greater +# than row are both active at the optimum, so the 'E' and 'G' crush paths +# (the latter negates the row) both affect the answer. +MIXED_SENSES = dict( + values=[1.0, 1.0, 1.0, 1.0, 1.0, 1.0], + indices=[0, 1, 1, 2, 0, 2], + offsets=[0, 2, 4, 6], + senses="ELG", + lower=[0.0, 0.0, 0.0], + upper=[10.0, 10.0, 10.0], +) + +# Row norms differ by seven orders of magnitude, so equilibration cannot leave +# row_scales at 1 and the crush has to divide the new RHS by the right ones. +BADLY_SCALED = dict( + values=[1e4, 1e4, 1e-3, -1e-3], + indices=[0, 1, 0, 1], + offsets=[0, 2, 4], + senses="GL", + lower=[0.0, 0.0], + upper=[10.0, 10.0], +) + +# Nonzero lower bounds make presolve translate x = x' + l, subtracting +# sum_j a_ij * l_j from each row's RHS. That is rhs_shift, -7 on the first row +# here, so dropping it would solve x0 + x1 >= b0 + 7 and miss the oracle. +LOWER_BOUNDED = dict( + values=[1.0, 1.0, 1.0, -1.0], + indices=[0, 1, 0, 1], + offsets=[0, 2, 4], + senses="GL", + lower=[3.0, 4.0], + upper=[20.0, 20.0], +) + +# Row 0 has no coefficients, so presolve drops it and it never reaches the +# barrier problem. Its RHS can only be checked for feasibility, not applied. +EMPTY_ROW = dict( + values=[1.0, 1.0], + indices=[0, 1], + offsets=[0, 0, 2], + senses="EG", + lower=[0.0, 0.0], + upper=[10.0, 10.0], +) + + +@pytest.mark.parametrize( + "model_args,first_rhs,updates", + [ + (MIXED_SENSES, [5.0, 8.0, 3.0], [[6.0, 7.0, 4.0], [4.0, 9.0, 2.5]]), + (BADLY_SCALED, [2e4, 5e-3], [[3e4, 2e-3], [1e4, 8e-3]]), + (LOWER_BOUNDED, [12.0, 1.0], [[15.0, 2.0], [9.0, 0.5]]), + ], + ids=["mixed_senses", "row_scaling", "rhs_shift"], +) +def test_update_rhs_matches_full_solve(model_args, first_rhs, updates, capfd): + """Reused solves must agree with a fresh full solve of the same model.""" + settings = _sequence_settings() + model = _build(**dict(model_args, rhs=first_rhs)) + + first, _ = _solve(model, settings, capfd) + assert first.get_termination_reason() == "Optimal" + + for rhs in updates: + model.update_rhs(np.asarray(rhs, dtype=np.float64)) + reused, log = _solve(model, settings, capfd) + assert REUSE_LOG in log, "update_rhs fell back to a full solve" + _assert_matches_oracle(reused, _full_solve(model_args, rhs)) + + +def test_update_rhs_keeps_dropped_empty_row(capfd): + """An empty row still satisfied by the new RHS must not block reuse.""" + settings = _sequence_settings() + model = _build(**dict(EMPTY_ROW, rhs=[0.0, 2.0])) + + first, log = _solve(model, settings, capfd) + assert first.get_termination_reason() == "Optimal" + assert "empty rows" in log, "presolve did not drop the empty row" + + rhs = [0.0, 4.0] + model.update_rhs(np.asarray(rhs, dtype=np.float64)) + reused, log = _solve(model, settings, capfd) + assert REUSE_LOG in log + _assert_matches_oracle(reused, _full_solve(EMPTY_ROW, rhs)) + + +def test_update_rhs_infeasible_empty_row_short_circuits(capfd): + """A dropped 'E' row needs 0 == b_i; violating it is infeasible. + + The row has no variables, so it is either satisfied for every x or for + none. That makes the verdict exact and lets the next Solve answer without + running IPM. The cache is kept, so a later feasible RHS still reuses it. + """ + settings = _sequence_settings() + model = _build(**dict(EMPTY_ROW, rhs=[0.0, 2.0])) + assert _solve(model, settings, capfd)[0].get_termination_reason() == ( + "Optimal" + ) + + model.update_rhs(np.array([1.0, 4.0])) + infeasible, log = _solve(model, settings, capfd) + assert infeasible.get_termination_reason() == "PrimalInfeasible" + assert RHS_INFEASIBLE_LOG in log + assert IPM_LOG not in log, "IPM ran despite a provably infeasible row" + + # Same cache, feasible RHS again. + rhs = [0.0, 4.0] + model.update_rhs(np.asarray(rhs, dtype=np.float64)) + recovered, log = _solve(model, settings, capfd) + assert REUSE_LOG in log, "cache was discarded by the infeasible update" + _assert_matches_oracle(recovered, _full_solve(EMPTY_ROW, rhs)) + + +# Lower bounds away from zero make presolve translate x = x' + l, which folds +# sum_j c_j * l_j into obj_constant. An RHS update must leave that alone. +TRANSLATED = dict( + values=[1.0, 1.0], + indices=[0, 1], + offsets=[0, 2], + senses="G", + lower=[3.0, 4.0], + upper=[20.0, 20.0], +) + + +def test_update_rhs_leaves_objective_constant_alone(capfd): + """An RHS update must not disturb the objective constant. + + obj_constant depends on c and on the translated lower bounds, not on b, so + the reused objective has to stay exact for a lower-bounded model. + """ + settings = _sequence_settings() + model = _build(**dict(TRANSLATED, rhs=[12.0], objective=[2.0, 0.0])) + + first, _ = _solve(model, settings, capfd) + assert first.get_termination_reason() == "Optimal" + + rhs = [15.0] + model.update_rhs(np.asarray(rhs, dtype=np.float64)) + reused, log = _solve(model, settings, capfd) + assert REUSE_LOG in log + _assert_matches_oracle( + reused, _full_solve(TRANSLATED, rhs, objective=[2.0, 0.0]) + ) + + +# x1 is free but the rows imply bounds on it, so presolve's free-variable +# bounding leaves state in presolve_info that the reuse path cannot replay. +FREE_VARIABLE = dict( + values=[1.0, 1.0, 1.0, 1.0], + indices=[0, 1, 0, 1], + offsets=[0, 2, 4], + senses="GL", + lower=[0.0, -np.inf], + upper=[10.0, np.inf], +) +FREE_VARIABLE_RHS = [4.0, 8.0] + + +@pytest.mark.parametrize( + "first_solve_bfv,expect_reuse", + [(None, True), (1, False)], + ids=["automatic_reuses", "explicit_bounding_refuses"], +) +def test_bounded_free_variables_block_reuse( + first_solve_bfv, expect_reuse, capfd +): + """A cache is only reusable if its own presolve left free variables alone. + + The two cases share a model, so the refusal in the second can only come + from the first solve having bounded a free variable. Without the pair, an + assertion that reuse did not happen would pass for any unrelated reason. + """ + settings = _sequence_settings() + if first_solve_bfv is not None: + settings.set_parameter( + "barrier_presolve_bound_free_variables", first_solve_bfv + ) + model = _build( + **dict(FREE_VARIABLE, rhs=FREE_VARIABLE_RHS, objective=[1.0, 0.0]) + ) + first, _ = _solve(model, settings, capfd) + assert first.get_termination_reason() == "Optimal" + + new_objective = [3.0, -2.0] + model.update_linear_objective(np.asarray(new_objective, dtype=np.float64)) + # Whatever the first solve asked for, this one asks for 0, which is what + # the gate used to key on all by itself. + settings.set_parameter("barrier_presolve_bound_free_variables", 0) + second, log = _solve(model, settings, capfd) + assert (REUSE_LOG in log) == expect_reuse + + _assert_matches_oracle( + second, + _full_solve(FREE_VARIABLE, FREE_VARIABLE_RHS, objective=new_objective), + ) + + +# Quadratic constraints reach the barrier as second-order cones. The conversion appends rows and +# permutes columns, so an update has to be mapped from model coordinates into the expanded layout +# the cache holds, and updates stay model-sized: the appended rows belong to the conversion. +# The conversion also rejects cone variables that carry an explicit upper bound or a nonzero +# lower bound, so the builders below leave the bounds of anything a cone touches open. + + +def _build_lorentz(rhs, objective, cone_head_free=False): + """``||(x1, x2)|| <= t`` plus two linear rows, no quadratic objective. + + Written as ``-t^2 + x1^2 + x2^2 <= 0``, which the conversion recognizes + and lifts by permuting ``(t, x1, x2)`` into a cone block, so all three are + conic variables. With ``cone_head_free``, ``t`` loses its lower bound and + row 1 becomes the singleton ``t >= b1`` that proves the head nonnegative. + """ + n = 3 + model = data_model.DataModel() + # row 0: x1 + x2 >= b0 ; row 1: t <= b1, or t >= b1 when the head is free + model.set_csr_constraint_matrix( + np.array([1.0, 1.0, 1.0], dtype=np.float64), + np.array([1, 2, 0], dtype=np.int32), + np.array([0, 2, 3], dtype=np.int32), + ) + model.set_constraint_bounds(np.asarray(rhs, dtype=np.float64)) + model.set_row_types("GG" if cone_head_free else "GL") + model.set_objective_coefficients(np.asarray(objective, dtype=np.float64)) + lower = np.zeros(n) + if cone_head_free: + lower[0] = -np.inf + model.set_variable_lower_bounds(lower) + model.set_variable_upper_bounds(np.full(n, np.inf)) + model.add_quadratic_constraint( + vals=np.array([-1.0, 1.0, 1.0]), + rows=np.array([0, 1, 2], dtype=np.int32), + cols=np.array([0, 1, 2], dtype=np.int32), + rhs_value=0.0, + sense="L", + ) + return model + + +def _build_general_cone(rhs, objective): + """``x1^2 + x2^2 + 2*x1 <= 8``, a shifted disk, plus two linear rows. + + The linear part and the nonzero RHS keep this off the recognized Lorentz + pattern, so the conversion takes its general path: it factors Q, adds cone + variables of its own, and appends four rows instead of two, which is a + different expanded shape to map an update through. The disk bounds + ``min x1``, so the model stays bounded. + """ + n = 2 + model = data_model.DataModel() + # row 0: x1 + x2 >= b0 ; row 1: x1 - x2 <= b1 + model.set_csr_constraint_matrix( + np.array([1.0, 1.0, 1.0, -1.0], dtype=np.float64), + np.array([0, 1, 0, 1], dtype=np.int32), + np.array([0, 2, 4], dtype=np.int32), + ) + model.set_constraint_bounds(np.asarray(rhs, dtype=np.float64)) + model.set_row_types("GL") + model.set_objective_coefficients(np.asarray(objective, dtype=np.float64)) + model.set_variable_lower_bounds(np.full(n, -np.inf)) + model.set_variable_upper_bounds(np.full(n, np.inf)) + model.add_quadratic_constraint( + vals=np.array([1.0, 1.0]), + rows=np.array([0, 1], dtype=np.int32), + cols=np.array([0, 1], dtype=np.int32), + linear_values=np.array([2.0]), + linear_indices=np.array([0], dtype=np.int32), + rhs_value=8.0, + sense="L", + ) + return model + + +def _full_cone_solve(build, rhs, objective, **kwargs): + """Oracle: a fresh cone model with default settings, no cache in play.""" + return solver.Solve( + build(rhs, objective, **kwargs), solver_settings.SolverSettings() + ) + + +# Builder, objective, first RHS, then the RHSs to update to. The updates stay +# feasible: on the shifted disk, x1 + x2 tops out just above 2.16. +CONE_CASES = { + "lorentz": ( + _build_lorentz, + [1.0, 0.0, 0.0], + [2.0, 9.0], + [[3.0, 9.0], [1.5, 9.0]], + ), + "general": ( + _build_general_cone, + [1.0, 0.0], + [1.0, 5.0], + [[1.5, 5.0], [0.5, 5.0]], + ), +} + + +@pytest.mark.parametrize("case", list(CONE_CASES)) +def test_cone_update_rhs_matches_full_solve(case, capfd): + """An RHS update on a cone model must agree with a fresh full solve. + + The model RHS has two entries while the converted problem has more rows, + so this fails outright if the update is not mapped into the expanded + layout, and gives a wrong answer if the appended rows lose their own RHS. + """ + build, objective, first_rhs, later_rhs = CONE_CASES[case] + settings = _sequence_settings() + model = build(first_rhs, objective) + + first, _ = _solve(model, settings, capfd) + assert first.get_termination_reason() == "Optimal" + + for rhs in later_rhs: + model.update_rhs(np.asarray(rhs, dtype=np.float64)) + reused, log = _solve(model, settings, capfd) + assert REUSE_LOG in log, "update_rhs fell back to a full solve" + _assert_matches_oracle(reused, _full_cone_solve(build, rhs, objective)) + + +@pytest.mark.parametrize("case", list(CONE_CASES)) +def test_cone_update_rhs_moves_the_optimum(case): + """Guard the oracles: the RHS being updated has to matter. + + Without this, a cached solve that ignored the new RHS would still match an + oracle that ignored it too, and every comparison above would pass. + """ + build, objective, first_rhs, later_rhs = CONE_CASES[case] + objectives = { + _full_cone_solve(build, rhs, objective).get_primal_objective() + for rhs in [first_rhs] + later_rhs + } + assert len(objectives) == 1 + len(later_rhs) + + +def test_lorentz_optimum_is_the_expected_one(): + """Pin one closed form, so the oracles are not the only reference. + + ``||(x1, x2)|| <= t`` with ``x1 + x2 >= b`` and ``x >= 0`` splits the sum + evenly, putting the optimum of ``min t`` at ``b / sqrt(2)``. + """ + for b in (2.0, 3.0): + solution = _full_cone_solve(_build_lorentz, [b, 9.0], [1.0, 0.0, 0.0]) + assert solution.get_termination_reason() == "Optimal" + assert solution.get_primal_objective() == pytest.approx( + b / np.sqrt(2.0), rel=1e-5, abs=1e-5 + ) + + +def test_cone_update_linear_objective_matches_full_solve(capfd): + """The conversion permutes columns, so the objective needs remapping too.""" + settings = _sequence_settings() + rhs = [2.0, 9.0] + model = _build_lorentz(rhs, [1.0, 0.0, 0.0]) + + first, _ = _solve(model, settings, capfd) + assert first.get_termination_reason() == "Optimal" + + new_objective = [1.0, 0.25, 0.0] + model.update_linear_objective(np.asarray(new_objective, dtype=np.float64)) + reused, log = _solve(model, settings, capfd) + assert REUSE_LOG in log, ( + "update_linear_objective fell back to a full solve" + ) + _assert_matches_oracle( + reused, _full_cone_solve(_build_lorentz, rhs, new_objective) + ) + + +def test_cone_update_rhs_rejects_lost_cone_head_bound(capfd): + """A cone head proved nonnegative by a row cannot lose that proof. + + The head here is free, so the conversion only accepts the model because + row 1 forces t >= 0. An RHS that relaxes the row to t >= -1 makes this a + model a full solve refuses, and reuse has to refuse it the same way rather + than solving a stale cone formulation. + """ + settings = _sequence_settings() + objective = [1.0, 0.0, 0.0] + model = _build_lorentz([2.0, 0.0], objective, cone_head_free=True) + + first, _ = _solve(model, settings, capfd) + assert first.get_termination_reason() == "Optimal" + + with pytest.raises(Exception, match="nonnegative"): + model.update_rhs(np.array([2.0, -1.0])) + + # A full solve of the same model does not return an optimum either, + # whether it reports the rejection as an exception or a status. + try: + oracle = _full_cone_solve( + _build_lorentz, [2.0, -1.0], objective, cone_head_free=True + ) + except Exception: + pass + else: + assert oracle.get_termination_reason() != "Optimal" + + +def test_cone_update_rhs_rejects_wrong_length(capfd): + """Length is validated against the model rows, not the converted rows.""" + settings = _sequence_settings() + model = _build_lorentz([2.0, 9.0], [1.0, 0.0, 0.0]) + first, _ = _solve(model, settings, capfd) + assert first.get_termination_reason() == "Optimal" + + # Two model rows. Passing the converted row count must not be accepted. + with pytest.raises(Exception, match="match the cached model row count"): + model.update_rhs(np.array([2.0, 9.0, 0.0, 0.0])) + + +def test_update_rhs_rejects_wrong_length(): + """Length is validated against the cached user row count.""" + settings = _sequence_settings() + model = _build(**dict(MIXED_SENSES, rhs=[5.0, 8.0, 3.0])) + assert solver.Solve(model, settings).get_termination_reason() == "Optimal" + + with pytest.raises(Exception, match="match the cached model row count"): + model.update_rhs(np.array([1.0, 2.0]))