diff --git a/cpp/src/branch_and_bound/CMakeLists.txt b/cpp/src/branch_and_bound/CMakeLists.txt index 1e40c1bbf1..a3c8fa064d 100644 --- a/cpp/src/branch_and_bound/CMakeLists.txt +++ b/cpp/src/branch_and_bound/CMakeLists.txt @@ -5,6 +5,7 @@ set(BRANCH_AND_BOUND_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/branch_and_bound.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/concurrent_root_solver.cu ${CMAKE_CURRENT_SOURCE_DIR}/pseudo_costs.cpp ${CMAKE_CURRENT_SOURCE_DIR}/diving_heuristics.cpp ) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index e2ea33bb6f..fe54aed873 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -6,6 +6,7 @@ /* clang-format on */ #include +#include #include #include #include @@ -3050,9 +3051,63 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( root_vstatus_, edge_norms_, nullptr); + // Dual simplex has finished; stop the GPU competitors if they are still running. + gpu_root_concurrent_halt_.store(1, std::memory_order_release); } - // Wait for the root relaxation solution to be sent by the diversity manager or dual simplex + // The diversity manager prepares the GPU problem while dual simplex starts on the CPU. + // Once the GPU problem is ready, launch PDLP and barrier from here so all root-LP + // competitors are owned by this function. + while (!concurrent_root_problem_ready_.load(std::memory_order_acquire) && + *get_root_concurrent_halt() == 0) { + std::this_thread::sleep_for(std::chrono::milliseconds(1)); +#pragma omp taskyield + } + + if (*get_root_concurrent_halt() == 0 && + concurrent_root_problem_ready_.load(std::memory_order_acquire)) { + cuopt_assert(concurrent_root_problem_ != nullptr, "Concurrent root problem is not configured"); + gpu_root_concurrent_halt_.store(0, std::memory_order_release); + try { + cuopt_assert(concurrent_root_settings_ != nullptr, + "Concurrent root settings are not configured"); + const f_t remaining_time = + std::max(settings_.time_limit - toc(exploration_stats_.start_time), 0); + const f_t root_time_limit = + std::min(concurrent_root_max_time_, remaining_time * concurrent_root_time_ratio_); + auto result = solve_concurrent_root_relaxation(concurrent_root_problem_, + *concurrent_root_settings_, + root_time_limit, + &gpu_root_concurrent_halt_); + if (result.usable) { + // Release the heuristics first: they only need the relaxation values, whereas + // crossover below can run until the time limit without ever producing a root. + if (root_lp_solution_callback_ != nullptr) { + root_lp_solution_callback_( + result.primal, result.dual, result.user_objective, result.optimal); + } + set_root_relaxation_solution(result.primal, + result.dual, + result.reduced_cost, + result.solver_objective, + result.user_objective, + result.iterations, + result.method); + // Same as the old diversity-manager path: an Optimal GPU root LP is a + // valid MIP dual bound even if dual simplex / crossover has not finished. + if (result.optimal) { update_user_bound(result.solver_objective); } + } else if (root_lp_solution_callback_ != nullptr) { + // No usable relaxation, but the heuristics must still be released so they can run + // without LP guidance rather than block on a dual simplex that may never finish. + root_lp_solution_callback_({}, {}, std::numeric_limits::infinity(), false); + } + } catch (const std::exception& e) { + settings_.log.printf("Concurrent GPU root LP failed: %s\n", e.what()); + } + } + + // Wait until either the GPU root solve supplies a crossover point or CPU dual + // simplex finishes. If dual simplex wins, stop and join the GPU solve. while (!root_crossover_solution_set_.load(std::memory_order_acquire) && *get_root_concurrent_halt() == 0) { if (received_halt_signal()) { diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 17ebdafa33..305c2c928c 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -83,6 +83,9 @@ struct clique_table_t; template struct mip_symmetry_t; +template +class problem_t; + template struct nondeterministic_policy_t; template @@ -146,7 +149,32 @@ class branch_and_bound_t { user_bound_callback_ = std::move(callback); } + // Hand the concurrent GPU root LP (PDLP/barrier) to the primal heuristics as soon as it + // finishes. The heuristics must not have to wait for crossover or dual simplex to complete + // the root, since on hard instances neither finishes within the time limit. + // Arguments are primal, dual, the user-space objective, and whether the LP proved optimality. + void set_root_lp_solution_callback( + std::function&, const std::vector&, f_t, bool)> callback) + { + root_lp_solution_callback_ = std::move(callback); + } + void set_concurrent_lp_root_solve(bool enable) { enable_concurrent_lp_root_solve_ = enable; } + void configure_concurrent_lp_root_solve(problem_t* problem, + const pdlp_solver_settings_t& settings, + f_t max_time, + f_t time_ratio) + { + concurrent_root_problem_ = problem; + concurrent_root_settings_ = std::make_unique>(settings); + concurrent_root_max_time_ = max_time; + concurrent_root_time_ratio_ = time_ratio; + enable_concurrent_lp_root_solve_ = true; + } + void notify_concurrent_root_problem_ready() + { + concurrent_root_problem_ready_.store(true, std::memory_order_release); + } // Seed the global upper bound from an external source (e.g., early FJ during presolve). // `bound` must be in B&B's internal objective space. @@ -257,6 +285,12 @@ class branch_and_bound_t { omp_atomic_t root_lp_current_lower_bound_; omp_atomic_t solving_root_relaxation_{false}; bool enable_concurrent_lp_root_solve_{false}; + problem_t* concurrent_root_problem_{nullptr}; + std::unique_ptr> concurrent_root_settings_; + f_t concurrent_root_max_time_{0}; + f_t concurrent_root_time_ratio_{0}; + std::atomic concurrent_root_problem_ready_{false}; + std::atomic gpu_root_concurrent_halt_{0}; std::atomic root_concurrent_halt_{0}; std::atomic node_concurrent_halt_{0}; bool is_root_solution_set{false}; @@ -293,6 +327,8 @@ class branch_and_bound_t { // corresponding subtree. omp_atomic_t lower_bound_numerical_; std::function user_bound_callback_; + std::function&, const std::vector&, f_t, bool)> + root_lp_solution_callback_; void print_table_header(); void report_heuristic(f_t obj, heuristics_origin_t origin); diff --git a/cpp/src/branch_and_bound/concurrent_root_solver.cu b/cpp/src/branch_and_bound/concurrent_root_solver.cu new file mode 100644 index 0000000000..88fa98c5a1 --- /dev/null +++ b/cpp/src/branch_and_bound/concurrent_root_solver.cu @@ -0,0 +1,61 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ + +#include +#include + +#include +#include + +#include + +namespace cuopt::mathematical_optimization::mip { + +template +concurrent_root_solution_t solve_concurrent_root_relaxation( + problem_t* problem, + const pdlp_solver_settings_t& settings, + f_t time_limit, + std::atomic* concurrent_halt) +{ + concurrent_root_solution_t result; + auto root_settings = settings; + root_settings.time_limit = time_limit; + root_settings.concurrent_halt = concurrent_halt; + + timer_t root_timer(time_limit); + auto lp_result = solve_lp_with_method(*problem, root_settings, root_timer); + const auto status = lp_result.get_termination_status(); + result.usable = + status != pdlp_termination_status_t::NumericalError && + status != pdlp_termination_status_t::ConcurrentLimit && + lp_result.get_primal_solution().size() == static_cast(problem->n_variables) && + lp_result.get_dual_solution().size() == static_cast(problem->n_constraints); + result.optimal = status == pdlp_termination_status_t::Optimal; + if (!result.usable) { return result; } + + auto& d_primal = lp_result.get_primal_solution(); + auto& d_dual = lp_result.get_dual_solution(); + auto& d_reduced_cost = lp_result.get_reduced_cost(); + result.primal.resize(d_primal.size()); + result.dual.resize(d_dual.size()); + result.reduced_cost.resize(d_reduced_cost.size()); + auto stream = problem->handle_ptr->get_stream(); + raft::copy(result.primal.data(), d_primal.data(), d_primal.size(), stream); + raft::copy(result.dual.data(), d_dual.data(), d_dual.size(), stream); + raft::copy(result.reduced_cost.data(), d_reduced_cost.data(), d_reduced_cost.size(), stream); + problem->handle_ptr->sync_stream(); + + result.user_objective = lp_result.get_objective_value(); + result.solver_objective = problem->get_solver_obj_from_user_obj(result.user_objective); + result.iterations = lp_result.get_additional_termination_information().number_of_steps_taken; + result.method = lp_result.get_additional_termination_information().solved_by; + return result; +} + +template concurrent_root_solution_t solve_concurrent_root_relaxation( + problem_t*, const pdlp_solver_settings_t&, double, std::atomic*); + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/branch_and_bound/concurrent_root_solver.hpp b/cpp/src/branch_and_bound/concurrent_root_solver.hpp new file mode 100644 index 0000000000..9591f7ecf7 --- /dev/null +++ b/cpp/src/branch_and_bound/concurrent_root_solver.hpp @@ -0,0 +1,37 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +#pragma once + +#include + +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +template +class problem_t; + +template +struct concurrent_root_solution_t { + bool usable{false}; + bool optimal{false}; + std::vector primal; + std::vector dual; + std::vector reduced_cost; + f_t solver_objective{0}; + f_t user_objective{0}; + i_t iterations{0}; + method_t method{method_t::Unset}; +}; + +template +concurrent_root_solution_t solve_concurrent_root_relaxation( + problem_t* problem, + const pdlp_solver_settings_t& settings, + f_t time_limit, + std::atomic* concurrent_halt); + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cu b/cpp/src/mip_heuristics/diversity/diversity_manager.cu index ec82c4b423..5d87bcde84 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -517,10 +517,7 @@ solution_t diversity_manager_t::run_solver() return population.best_feasible(); } - population.timer = timer; - const f_t time_limit = timer.remaining_time(); - const auto& hp = context.settings.heuristic_params; - const f_t lp_time_limit = std::min(hp.root_lp_max_time, time_limit * hp.root_lp_time_ratio); + population.timer = timer; // after every change to the problem, we should resize all the relevant vars // we need to encapsulate that to prevent repetitions recombine_stats.reset(); @@ -558,135 +555,164 @@ solution_t diversity_manager_t::run_solver() } else if (!fj_only_run) { convert_greater_to_less(*problem_ptr); - f_t absolute_tolerance = context.settings.tolerances.absolute_tolerance; - - pdlp_solver_settings_t pdlp_settings{}; - pdlp_settings.tolerances.absolute_dual_tolerance = absolute_tolerance; - pdlp_settings.tolerances.relative_dual_tolerance = - context.settings.tolerances.relative_tolerance; - pdlp_settings.tolerances.absolute_primal_tolerance = absolute_tolerance; - pdlp_settings.tolerances.relative_primal_tolerance = - context.settings.tolerances.relative_tolerance; - pdlp_settings.time_limit = lp_time_limit; - pdlp_settings.first_primal_feasible = false; - pdlp_settings.concurrent_halt = &global_concurrent_halt; - pdlp_settings.method = method_t::Concurrent; - pdlp_settings.inside_mip = true; - pdlp_settings.pdlp_solver_mode = pdlp_solver_mode_t::Stable2; - pdlp_settings.num_gpus = context.settings.num_gpus; - pdlp_settings.presolver = presolver_t::None; - pdlp_settings.per_constraint_residual = true; - set_pdlp_solver_mode(pdlp_settings); - timer_t lp_timer(lp_time_limit); - auto lp_result = solve_lp_with_method(*problem_ptr, pdlp_settings, lp_timer); - - // The concurrent root LP can fail to produce a usable solution -- e.g. the barrier - // hits a numerical error on an infeasible problem and PDLP returns NumericalError - // with empty primal/dual. In that case we must not copy or hand off the empty - // result (copying n elements from an empty buffer throws), and we must still - // release B&B's root-relaxation wait so it proceeds with its own dual-simplex root - // instead of spinning forever. - const bool root_lp_usable = - lp_result.get_termination_status() != pdlp_termination_status_t::NumericalError && - lp_result.get_primal_solution().size() == lp_optimal_solution.size() && - lp_result.get_dual_solution().size() == lp_dual_optimal_solution.size(); - - bool use_staged_simplex_solution = false; - { - std::lock_guard guard(relaxed_solution_mutex); - use_staged_simplex_solution = simplex_solution_exists.load(); - if (!use_staged_simplex_solution && root_lp_usable) { - raft::copy(lp_optimal_solution.data(), - lp_result.get_primal_solution().data(), - lp_optimal_solution.size(), - problem_ptr->handle_ptr->get_stream()); - raft::copy(lp_dual_optimal_solution.data(), - lp_result.get_dual_solution().data(), - lp_dual_optimal_solution.size(), - problem_ptr->handle_ptr->get_stream()); + if (context.branch_and_bound_ptr != nullptr && + context.branch_and_bound_ptr->enable_concurrent_lp_root_solve()) { + // B&B owns the CPU/GPU concurrent root solve. Signal that the GPU problem is ready, + // then wait for the first root relaxation B&B publishes: either the finished root + // (set_simplex_solution_callback) or the raw GPU LP values it hands over before + // crossover (set_root_lp_solution_callback). Waiting only for the finished root would + // starve the LP-guided heuristics on instances where neither dual simplex nor + // crossover completes within the time limit. + context.branch_and_bound_ptr->notify_concurrent_root_problem_ready(); + while (!simplex_solution_exists.load(std::memory_order_acquire) && + !root_lp_solution_exists.load(std::memory_order_acquire) && !check_b_b_preemption()) { + std::this_thread::sleep_for(std::chrono::milliseconds(1)); } - } - if (use_staged_simplex_solution) { consume_staged_simplex_solution(lp_state); } - if (use_staged_simplex_solution || root_lp_usable) { - cuopt_assert(thrust::all_of(problem_ptr->handle_ptr->get_thrust_policy(), - lp_optimal_solution.begin(), - lp_optimal_solution.end(), - [] __host__ __device__(f_t val) { return std::isfinite(val); }), - "LP optimal solution contains non-finite values"); - } - ls.lp_optimal_exists = true; - if (!use_staged_simplex_solution) { - if (!root_lp_usable) { - // The concurrent root LP produced no usable solution. Do not hand an empty - // solution to B&B; instead release its root-relaxation wait loop so it falls - // back to its own dual-simplex root rather than deadlocking. - CUOPT_LOG_DEBUG("Root LP produced no usable solution (status %d); releasing B&B root solve", - (int)lp_result.get_termination_status()); + if (simplex_solution_exists.load(std::memory_order_acquire)) { + consume_staged_simplex_solution(lp_state); + ls.lp_optimal_exists = true; + } else if (root_lp_solution_exists.load(std::memory_order_acquire)) { + ls.lp_optimal_exists = consume_staged_root_lp_solution(); + } else { ls.lp_optimal_exists = false; - if (context.branch_and_bound_ptr != nullptr) { - context.branch_and_bound_ptr->set_root_concurrent_halt(1); + } + } else { + // Heuristics-only fallback: no B&B object exists to own the root solve. + const f_t time_limit = timer.remaining_time(); + const auto& hp = context.settings.heuristic_params; + const f_t lp_time_limit = std::min(hp.root_lp_max_time, time_limit * hp.root_lp_time_ratio); + f_t absolute_tolerance = context.settings.tolerances.absolute_tolerance; + + pdlp_solver_settings_t pdlp_settings{}; + pdlp_settings.tolerances.absolute_dual_tolerance = absolute_tolerance; + pdlp_settings.tolerances.relative_dual_tolerance = + context.settings.tolerances.relative_tolerance; + pdlp_settings.tolerances.absolute_primal_tolerance = absolute_tolerance; + pdlp_settings.tolerances.relative_primal_tolerance = + context.settings.tolerances.relative_tolerance; + pdlp_settings.time_limit = lp_time_limit; + pdlp_settings.first_primal_feasible = false; + pdlp_settings.concurrent_halt = &global_concurrent_halt; + pdlp_settings.method = method_t::Concurrent; + pdlp_settings.inside_mip = true; + pdlp_settings.pdlp_solver_mode = pdlp_solver_mode_t::Stable2; + pdlp_settings.num_gpus = context.settings.num_gpus; + pdlp_settings.presolver = presolver_t::None; + pdlp_settings.per_constraint_residual = true; + set_pdlp_solver_mode(pdlp_settings); + timer_t lp_timer(lp_time_limit); + auto lp_result = solve_lp_with_method(*problem_ptr, pdlp_settings, lp_timer); + + // The concurrent root LP can fail to produce a usable solution -- e.g. the barrier + // hits a numerical error on an infeasible problem and PDLP returns NumericalError + // with empty primal/dual. In that case we must not copy or hand off the empty + // result (copying n elements from an empty buffer throws), and we must still + // release B&B's root-relaxation wait so it proceeds with its own dual-simplex root + // instead of spinning forever. + const bool root_lp_usable = + lp_result.get_termination_status() != pdlp_termination_status_t::NumericalError && + lp_result.get_primal_solution().size() == lp_optimal_solution.size() && + lp_result.get_dual_solution().size() == lp_dual_optimal_solution.size(); + + bool use_staged_simplex_solution = false; + { + std::lock_guard guard(relaxed_solution_mutex); + use_staged_simplex_solution = simplex_solution_exists.load(); + if (!use_staged_simplex_solution && root_lp_usable) { + raft::copy(lp_optimal_solution.data(), + lp_result.get_primal_solution().data(), + lp_optimal_solution.size(), + problem_ptr->handle_ptr->get_stream()); + raft::copy(lp_dual_optimal_solution.data(), + lp_result.get_dual_solution().data(), + lp_dual_optimal_solution.size(), + problem_ptr->handle_ptr->get_stream()); + } + } + if (use_staged_simplex_solution) { consume_staged_simplex_solution(lp_state); } + if (use_staged_simplex_solution || root_lp_usable) { + cuopt_assert(thrust::all_of(problem_ptr->handle_ptr->get_thrust_policy(), + lp_optimal_solution.begin(), + lp_optimal_solution.end(), + [] __host__ __device__(f_t val) { return std::isfinite(val); }), + "LP optimal solution contains non-finite values"); + } + ls.lp_optimal_exists = true; + if (!use_staged_simplex_solution) { + if (!root_lp_usable) { + // The concurrent root LP produced no usable solution. Do not hand an empty + // solution to B&B; instead release its root-relaxation wait loop so it falls + // back to its own dual-simplex root rather than deadlocking. + CUOPT_LOG_DEBUG( + "Root LP produced no usable solution (status %d); releasing B&B root solve", + (int)lp_result.get_termination_status()); + ls.lp_optimal_exists = false; + if (context.branch_and_bound_ptr != nullptr) { + context.branch_and_bound_ptr->set_root_concurrent_halt(1); + } + } else if (lp_result.get_termination_status() == pdlp_termination_status_t::Optimal) { + solution_t lp_sol(*problem_ptr); + lp_sol.copy_new_assignment(lp_optimal_solution); + const bool consider_integrality = false; + lp_sol.compute_feasibility(consider_integrality); + if (lp_sol.get_feasible()) { set_new_user_bound(lp_result.get_objective_value()); } + } else if (lp_result.get_termination_status() == + pdlp_termination_status_t::PrimalInfeasible) { + CUOPT_LOG_ERROR("Problem is primal infeasible, continuing anyway!"); + ls.lp_optimal_exists = false; + } else if (lp_result.get_termination_status() == + pdlp_termination_status_t::DualInfeasible) { + CUOPT_LOG_ERROR("PDLP detected dual infeasibility, continuing anyway!"); + ls.lp_optimal_exists = false; + } else if (lp_result.get_termination_status() == pdlp_termination_status_t::TimeLimit) { + CUOPT_LOG_DEBUG( + "Initial LP run exceeded time limit, continuing solver with partial LP result!"); + // note to developer, in debug mode the LP run might be too slow and it might cause PDLP + // not to bring variables within the bounds } - } else if (lp_result.get_termination_status() == pdlp_termination_status_t::Optimal) { - solution_t lp_sol(*problem_ptr); - lp_sol.copy_new_assignment(lp_optimal_solution); - const bool consider_integrality = false; - lp_sol.compute_feasibility(consider_integrality); - if (lp_sol.get_feasible()) { set_new_user_bound(lp_result.get_objective_value()); } - } else if (lp_result.get_termination_status() == - pdlp_termination_status_t::PrimalInfeasible) { - CUOPT_LOG_ERROR("Problem is primal infeasible, continuing anyway!"); - ls.lp_optimal_exists = false; - } else if (lp_result.get_termination_status() == pdlp_termination_status_t::DualInfeasible) { - CUOPT_LOG_ERROR("PDLP detected dual infeasibility, continuing anyway!"); - ls.lp_optimal_exists = false; - } else if (lp_result.get_termination_status() == pdlp_termination_status_t::TimeLimit) { - CUOPT_LOG_DEBUG( - "Initial LP run exceeded time limit, continuing solver with partial LP result!"); - // note to developer, in debug mode the LP run might be too slow and it might cause PDLP - // not to bring variables within the bounds } - } - // Hand the root relaxation off to branch and bound when we have a usable solution - // (sets root_crossover_solution_set_, releasing B&B's wait). When the root LP failed - // the wait is instead released above via set_root_concurrent_halt, and a staged - // dual-simplex solution is owned by B&B already, so neither needs this hand-off. - if (!use_staged_simplex_solution && root_lp_usable && - problem_ptr->set_root_relaxation_solution_callback != nullptr) { - auto& d_primal_solution = lp_result.get_primal_solution(); - auto& d_dual_solution = lp_result.get_dual_solution(); - auto& d_reduced_costs = lp_result.get_reduced_cost(); - - std::vector host_primal(d_primal_solution.size()); - std::vector host_dual(d_dual_solution.size()); - std::vector host_reduced_costs(d_reduced_costs.size()); - raft::copy(host_primal.data(), - d_primal_solution.data(), - d_primal_solution.size(), - problem_ptr->handle_ptr->get_stream()); - raft::copy(host_dual.data(), - d_dual_solution.data(), - d_dual_solution.size(), - problem_ptr->handle_ptr->get_stream()); - raft::copy(host_reduced_costs.data(), - d_reduced_costs.data(), - d_reduced_costs.size(), - problem_ptr->handle_ptr->get_stream()); - problem_ptr->handle_ptr->sync_stream(); - - // PDLP returns user-space objective (it applies objective_scaling_factor internally) - auto user_obj = lp_result.get_objective_value(); - auto solver_obj = problem_ptr->get_solver_obj_from_user_obj(user_obj); - auto iterations = lp_result.get_additional_termination_information().number_of_steps_taken; - auto method = lp_result.get_additional_termination_information().solved_by; - // Set for the B&B (param4 expects solver space, param5 expects user space) - problem_ptr->set_root_relaxation_solution_callback( - host_primal, host_dual, host_reduced_costs, solver_obj, user_obj, iterations, method); - } + // Hand the root relaxation off to branch and bound when we have a usable solution + // (sets root_crossover_solution_set_, releasing B&B's wait). When the root LP failed + // the wait is instead released above via set_root_concurrent_halt, and a staged + // dual-simplex solution is owned by B&B already, so neither needs this hand-off. + if (!use_staged_simplex_solution && root_lp_usable && + problem_ptr->set_root_relaxation_solution_callback != nullptr) { + auto& d_primal_solution = lp_result.get_primal_solution(); + auto& d_dual_solution = lp_result.get_dual_solution(); + auto& d_reduced_costs = lp_result.get_reduced_cost(); + + std::vector host_primal(d_primal_solution.size()); + std::vector host_dual(d_dual_solution.size()); + std::vector host_reduced_costs(d_reduced_costs.size()); + raft::copy(host_primal.data(), + d_primal_solution.data(), + d_primal_solution.size(), + problem_ptr->handle_ptr->get_stream()); + raft::copy(host_dual.data(), + d_dual_solution.data(), + d_dual_solution.size(), + problem_ptr->handle_ptr->get_stream()); + raft::copy(host_reduced_costs.data(), + d_reduced_costs.data(), + d_reduced_costs.size(), + problem_ptr->handle_ptr->get_stream()); + problem_ptr->handle_ptr->sync_stream(); + + // PDLP returns user-space objective (it applies objective_scaling_factor internally) + auto user_obj = lp_result.get_objective_value(); + auto solver_obj = problem_ptr->get_solver_obj_from_user_obj(user_obj); + auto iterations = lp_result.get_additional_termination_information().number_of_steps_taken; + auto method = lp_result.get_additional_termination_information().solved_by; + // Set for the B&B (param4 expects solver space, param5 expects user space) + problem_ptr->set_root_relaxation_solution_callback( + host_primal, host_dual, host_reduced_costs, solver_obj, user_obj, iterations, method); + } - if (!use_staged_simplex_solution && root_lp_usable) { - // in case the pdlp returned var boudns that are out of bounds - clamp_within_var_bounds(lp_optimal_solution, problem_ptr, problem_ptr->handle_ptr); + if (!use_staged_simplex_solution && root_lp_usable) { + // in case the pdlp returned var boudns that are out of bounds + clamp_within_var_bounds(lp_optimal_solution, problem_ptr, problem_ptr->handle_ptr); + } } } @@ -1041,6 +1067,77 @@ void diversity_manager_t::set_simplex_solution(const std::vector& CUOPT_LOG_DEBUG("Staged simplex solution and requested concurrent halt"); } +template +void diversity_manager_t::set_root_lp_solution(const std::vector& solution, + const std::vector& dual_solution, + f_t user_objective, + bool optimal) +{ + // Called from the B&B thread, so only stage the host values here. The device copies + // happen on the heuristics thread in consume_staged_root_lp_solution. + std::lock_guard lock(relaxed_solution_mutex); + cuopt_assert(solution.empty() || lp_optimal_solution.size() == solution.size(), + "Assignment size mismatch"); + cuopt_assert(solution.empty() || problem_ptr->n_constraints == dual_solution.size(), + "Dual assignment size mismatch"); + staged_root_lp_solution = solution; + staged_root_lp_dual_solution = dual_solution; + staged_root_lp_objective = user_objective; + staged_root_lp_optimal = optimal; + root_lp_solution_exists.store(true, std::memory_order_release); + CUOPT_LOG_DEBUG( + "Staged GPU root LP solution with objective %f (optimal %d)", user_objective, (int)optimal); +} + +template +bool diversity_manager_t::consume_staged_root_lp_solution() +{ + std::vector primal_local; + std::vector dual_local; + f_t objective_local = std::numeric_limits::infinity(); + bool optimal_local = false; + { + std::lock_guard guard(relaxed_solution_mutex); + cuopt_assert(root_lp_solution_exists.load(), + "Root LP solution flag set without a staged root LP solution"); + primal_local = staged_root_lp_solution; + dual_local = staged_root_lp_dual_solution; + objective_local = staged_root_lp_objective; + optimal_local = staged_root_lp_optimal; + } + // An empty hand-off means the root LP produced nothing usable; the heuristics run + // without LP guidance instead of waiting for a root that may never arrive. + if (primal_local.empty()) { + CUOPT_LOG_DEBUG("Root LP produced no usable solution; continuing without LP guidance"); + return false; + } + cuopt_assert(lp_optimal_solution.size() == primal_local.size(), "Assignment size mismatch"); + cuopt_assert(lp_dual_optimal_solution.size() == dual_local.size(), + "Dual assignment size mismatch"); + auto stream = problem_ptr->handle_ptr->get_stream(); + raft::copy(lp_optimal_solution.data(), primal_local.data(), lp_optimal_solution.size(), stream); + raft::copy( + lp_dual_optimal_solution.data(), dual_local.data(), lp_dual_optimal_solution.size(), stream); + problem_ptr->handle_ptr->sync_stream(); + cuopt_assert(thrust::all_of(problem_ptr->handle_ptr->get_thrust_policy(), + lp_optimal_solution.begin(), + lp_optimal_solution.end(), + [] __host__ __device__(f_t val) { return std::isfinite(val); }), + "LP optimal solution contains non-finite values"); + // Only an LP that proved optimality gives a valid dual bound for the MIP. A partial + // relaxation is still useful to guide the heuristics, but must not move the bound. + if (optimal_local) { + solution_t lp_sol(*problem_ptr); + lp_sol.copy_new_assignment(lp_optimal_solution); + const bool consider_integrality = false; + lp_sol.compute_feasibility(consider_integrality); + if (lp_sol.get_feasible()) { set_new_user_bound(objective_local); } + } + // PDLP can return values slightly outside the variable bounds. + clamp_within_var_bounds(lp_optimal_solution, problem_ptr, problem_ptr->handle_ptr); + return true; +} + #if MIP_INSTANTIATE_FLOAT template class diversity_manager_t; #endif diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cuh b/cpp/src/mip_heuristics/diversity/diversity_manager.cuh index 7655a58eb8..8ee9603968 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cuh +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cuh @@ -71,6 +71,12 @@ class diversity_manager_t { void set_simplex_solution(const std::vector& solution, const std::vector& dual_solution, f_t objective); + // Returns whether a usable LP relaxation was consumed. + bool consume_staged_root_lp_solution(); + void set_root_lp_solution(const std::vector& solution, + const std::vector& dual_solution, + f_t user_objective, + bool optimal); mip_solver_context_t& context; mip::branch_and_bound_t* branch_and_bound_ptr; problem_t* problem_ptr; @@ -82,6 +88,12 @@ class diversity_manager_t { std::vector staged_simplex_solution; std::vector staged_simplex_dual_solution; f_t staged_simplex_objective{std::numeric_limits::infinity()}; + // Concurrent GPU root LP handed over by branch and bound before crossover runs. + std::atomic root_lp_solution_exists{false}; + std::vector staged_root_lp_solution; + std::vector staged_root_lp_dual_solution; + f_t staged_root_lp_objective{std::numeric_limits::infinity()}; + bool staged_root_lp_optimal{false}; local_search_t ls; cuopt::timer_t timer; bound_prop_recombiner_t bound_prop_recombiner; diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index f8eac0c4d8..f716eb4e75 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -442,7 +442,36 @@ solution_t mip_solver_t::run_solver() // Set the primal heuristics -> branch and bound callback if (context.settings.determinism_mode == CUOPT_MODE_OPPORTUNISTIC) { - branch_and_bound->set_concurrent_lp_root_solve(true); + pdlp_solver_settings_t concurrent_root_settings{}; + concurrent_root_settings.tolerances.absolute_dual_tolerance = + context.settings.tolerances.absolute_tolerance; + concurrent_root_settings.tolerances.relative_dual_tolerance = + context.settings.tolerances.relative_tolerance; + concurrent_root_settings.tolerances.absolute_primal_tolerance = + context.settings.tolerances.absolute_tolerance; + concurrent_root_settings.tolerances.relative_primal_tolerance = + context.settings.tolerances.relative_tolerance; + concurrent_root_settings.first_primal_feasible = false; + concurrent_root_settings.method = method_t::Concurrent; + concurrent_root_settings.inside_mip = true; + concurrent_root_settings.pdlp_solver_mode = pdlp_solver_mode_t::Stable2; + concurrent_root_settings.num_gpus = context.settings.num_gpus; + concurrent_root_settings.presolver = presolver_t::None; + concurrent_root_settings.per_constraint_residual = true; + set_pdlp_solver_mode(concurrent_root_settings); + branch_and_bound->configure_concurrent_lp_root_solve( + context.problem_ptr, + concurrent_root_settings, + context.settings.heuristic_params.root_lp_max_time, + context.settings.heuristic_params.root_lp_time_ratio); + + branch_and_bound->set_root_lp_solution_callback( + std::bind(&diversity_manager_t::set_root_lp_solution, + &dm, + std::placeholders::_1, + std::placeholders::_2, + std::placeholders::_3, + std::placeholders::_4)); context.problem_ptr->branch_and_bound_callback = std::bind(&mip::branch_and_bound_t::set_solution_from_heuristics, @@ -461,17 +490,6 @@ solution_t mip_solver_t::run_solver() context.work_unit_scheduler_.register_context(branch_and_bound->get_work_unit_context()); // context.work_unit_scheduler_.verbose = true; - context.problem_ptr->set_root_relaxation_solution_callback = - std::bind(&mip::branch_and_bound_t::set_root_relaxation_solution, - branch_and_bound.get(), - std::placeholders::_1, - std::placeholders::_2, - std::placeholders::_3, - std::placeholders::_4, - std::placeholders::_5, - std::placeholders::_6, - std::placeholders::_7); - if (timer_.check_time_limit()) { CUOPT_LOG_INFO("Time limit reached during B&B setup"); context.stats.total_solve_time = timer_.elapsed_time();