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/barrier.cu b/cpp/src/barrier/barrier.cu index 24dfd3c411..b33a1704ff 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -895,7 +895,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) @@ -4911,6 +4911,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_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..e33eb7d295 100644 --- a/cpp/src/barrier/barrier_transform.hpp +++ b/cpp/src/barrier/barrier_transform.hpp @@ -10,8 +10,10 @@ #include #include +#include #include #include +#include #include namespace cuopt::mathematical_optimization { @@ -19,8 +21,8 @@ namespace cuopt::mathematical_optimization { /** * 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}; @@ -43,6 +45,12 @@ struct barrier_transform_t { 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. @@ -109,4 +117,77 @@ 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 != xf.user_num_rows) { + throw std::invalid_argument("update_rhs: RHS length must match the cached user 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."); + } + + // 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 < m; ++i) { + original[static_cast(i)] = + xf.row_sense[static_cast(i)] == 'G' ? -b[i] : b[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 >= m) { + 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/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 6f6ecc88f9..532586cb12 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -36,6 +36,7 @@ #include #include #include +#include #include namespace cuopt::mathematical_optimization::simplex { @@ -91,6 +92,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 +436,22 @@ 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 = + auto const* xf = (cache != nullptr && cache->dirty()) ? cache->transform() : nullptr; + const bool reuse_cached_data = 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() && + // 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( @@ -446,7 +470,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(); } @@ -503,8 +527,13 @@ lp_status_t solve_linear_program_with_barrier( 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->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. + xf->rhs_update_supported = + user_problem.num_range_rows == 0 && !presolve_info.folding_info.is_folded; + xf->barrier_lp = std::make_unique>(barrier_lp); + solver_lp = xf->barrier_lp.get(); cache->store_transform(std::move(xf)); } @@ -514,19 +543,26 @@ lp_status_t solve_linear_program_with_barrier( 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. try { - auto crushed = cuopt::mathematical_optimization::crush_user_linear_objective( + auto const 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)]; - } - } + 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(), user_problem.num_rows); + 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/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 63444ab5bb..e5a18af8d7 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -518,6 +518,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, @@ -530,31 +540,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; @@ -1880,10 +1889,14 @@ 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. 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() && + effective_bound_free_variables(settings) == 0 && + xf->presolve_info.bounded_free_variables.empty() && op_problem.has_quadratic_objective() && !op_problem.has_quadratic_constraints() && xf->second_order_cone_dims.empty() && static_cast(xf->row_sense.size()) == xf->user_num_rows && op_problem.get_n_variables() == xf->user_num_cols && @@ -1930,6 +1943,9 @@ 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(); } 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..8f3fa0687a 100644 --- a/python/cuopt/cuopt/linear_programming/data_model/data_model.py +++ b/python/cuopt/cuopt/linear_programming/data_model/data_model.py @@ -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, and bounds must stay unchanged). + Cache reuse is QP-only: quadratic constraints take a full solve. + + 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..fdd4d22fef --- /dev/null +++ b/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py @@ -0,0 +1,300 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 + +"""Barrier cache reuse (``sequence_solve``) for ``DataModel.update_rhs``. + +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 + ), + ) + + +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 user row count"): + model.update_rhs(np.array([1.0, 2.0]))