From 9698d00fe1bb6b19b11506b56a3acbe352948992 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Mon, 27 Apr 2026 05:28:31 -0700 Subject: [PATCH 1/6] run a single burst --- cpp/src/branch_and_bound/branch_and_bound.cpp | 24 ++ .../feasibility_jump/cpu_fj_thread.cuh | 59 +++ .../mip_heuristics/feasibility_jump/fj_cpu.cu | 357 ++++++++++++++++-- .../feasibility_jump/fj_cpu.cuh | 20 +- 4 files changed, 401 insertions(+), 59 deletions(-) create mode 100644 cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index e69ff7b9a5..2ecd01558e 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include @@ -29,6 +30,7 @@ #include #include +#include #include #include #include @@ -2224,6 +2226,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t last_objective = root_objective_; f_t root_relax_objective = root_objective_; + constexpr bool enable_root_cut_cpufj = true; + constexpr double root_cut_cpufj_work_units = 0.1; + f_t cut_generation_start_time = tic(); i_t cut_pool_size = 0; for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { @@ -2504,6 +2509,25 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut original_lp_.A.col_start[original_lp_.A.n]); } + if (enable_root_cut_cpufj && original_lp_.num_rows > original_rows) { + std::atomic root_cut_cpufj_preemption{false}; + detail::run_fj_cpu_from_host_lp( + original_lp_, + var_types_, + root_relax_soln_.x, + settings_, + root_cut_cpufj_preemption, + f_t{1}, + 1, + [this](f_t, const std::vector& assignment, double obj) { + std::vector user_assignment(assignment.begin(), + assignment.begin() + original_problem_.num_cols); + CUOPT_LOG_INFO("Root cut CPUFJ found solution with objective %.16e\n", obj); + set_new_solution(user_assignment); + }, + "[RootCut CPUFJ] "); + } + set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); pc_.resize(original_lp_.num_cols); diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh b/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh new file mode 100644 index 0000000000..8de191cc71 --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh @@ -0,0 +1,59 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace cuopt::linear_programming::detail { + +template +struct fj_cpu_climber_t; + +template +class fj_t; + +template +struct cpu_fj_thread_t : public cpu_worker_thread_base_t> { + ~cpu_fj_thread_t(); + + void run_worker(); + void on_terminate(); + void on_start(); + bool get_result() { return cpu_fj_solution_found; } + + void stop_cpu_solver(); + + std::atomic cpu_fj_solution_found{false}; + f_t time_limit{+std::numeric_limits::infinity()}; + double work_unit_limit{std::numeric_limits::infinity()}; + std::unique_ptr> fj_cpu; + fj_t* fj_ptr{nullptr}; +}; + +template +bool run_fj_cpu_from_host_lp( + const dual_simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const dual_simplex::simplex_solver_settings_t& settings, + std::atomic& preemption_flag, + f_t time_limit, + double work_unit_limit, + std::function&, double)> improvement_callback, + std::string log_prefix); + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 4eaa5b6a21..9175e4291c 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -7,6 +7,10 @@ #include +#include +#include + +#include "cpu_fj_thread.cuh" #include "feasibility_jump.cuh" #include "feasibility_jump_impl_common.cuh" #include "fj_cpu.cuh" @@ -18,9 +22,12 @@ #include #include +#include #include +#include #include #include +#include #include #include #include @@ -41,6 +48,53 @@ namespace cuopt::linear_programming::detail { +template +static f_t clamp_value(f_t value, f_t lower, f_t upper) +{ + return std::min(std::max(value, lower), upper); +} + +template +static void rebuild_reverse_matrix(i_t n_variables, + i_t n_constraints, + i_t nnz, + const std::vector& coefficients, + const std::vector& variables, + const std::vector& offsets, + std::vector& reverse_coefficients, + std::vector& reverse_constraints, + std::vector& reverse_offsets) +{ + reverse_offsets.assign(n_variables + 1, 0); + for (i_t row = 0; row < n_constraints; ++row) { + for (i_t p = offsets[row]; p < offsets[row + 1]; ++p) { + ++reverse_offsets[variables[p] + 1]; + } + } + std::partial_sum(reverse_offsets.begin(), reverse_offsets.end(), reverse_offsets.begin()); + + reverse_constraints.resize(nnz); + reverse_coefficients.resize(nnz); + std::vector next = reverse_offsets; + for (i_t row = 0; row < n_constraints; ++row) { + for (i_t p = offsets[row]; p < offsets[row + 1]; ++p) { + const i_t col = variables[p]; + const i_t dst = next[col]++; + reverse_constraints[dst] = row; + reverse_coefficients[dst] = coefficients[p]; + } + } +} + +template +void finalize_fj_cpu_host_initialization( + fj_cpu_climber_t& fj_cpu, + i_t n_variables, + i_t n_constraints, + i_t n_integer_vars, + i_t nnz, + const typename mip_solver_settings_t::tolerances_t& tolerances); + template thrust::tuple get_mtm_for_bound(const typename fj_t::climber_data_t::view_t& fj, i_t var_idx, @@ -792,9 +846,8 @@ static void apply_move(fj_cpu_climber_t& fj_cpu, fj_cpu.h_incumbent_objective - fj_cpu.settings.parameters.breakthrough_move_epsilon; fj_cpu.h_best_assignment = fj_cpu.h_assignment; fj_cpu.iterations_since_best = 0; - CUOPT_LOG_TRACE("%sCPUFJ: new best objective: %g", - fj_cpu.log_prefix.c_str(), - fj_cpu.pb_ptr->get_user_obj_from_solver_obj(fj_cpu.h_incumbent_objective)); + CUOPT_LOG_TRACE( + "%sCPUFJ: new best objective: %g", fj_cpu.log_prefix.c_str(), fj_cpu.h_incumbent_objective); if (fj_cpu.improvement_callback) { double current_work_units = fj_cpu.work_units_elapsed.load(std::memory_order_acquire); fj_cpu.improvement_callback( @@ -829,7 +882,6 @@ static thrust::tuple find_mtm_move( fj_cpu_climber_t& fj_cpu, const std::vector& target_cstrs, bool localmin = false) { CPUFJ_NVTX_RANGE("CPUFJ::find_mtm_move"); - auto& problem = *fj_cpu.pb_ptr; raft::random::PCGenerator rng(fj_cpu.settings.seed + fj_cpu.iterations, 0, 0); @@ -1258,33 +1310,29 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, fj_cpu.h_tabu_lastinc.resize(fj_cpu.pb_ptr->n_variables, 0); fj_cpu.iterations = 0; - // set pointers to host copies - // technically not 'device_span's but raft doesn't have a universal span. - // cuda::std::span? - fj_cpu.view.cstr_left_weights = - raft::device_span(fj_cpu.h_cstr_left_weights.data(), fj_cpu.h_cstr_left_weights.size()); - fj_cpu.view.cstr_right_weights = - raft::device_span(fj_cpu.h_cstr_right_weights.data(), fj_cpu.h_cstr_right_weights.size()); - fj_cpu.view.objective_weight = &fj_cpu.h_objective_weight; - fj_cpu.view.incumbent_assignment = - raft::device_span(fj_cpu.h_assignment.data(), fj_cpu.h_assignment.size()); - fj_cpu.view.incumbent_lhs = raft::device_span(fj_cpu.h_lhs.data(), fj_cpu.h_lhs.size()); - fj_cpu.view.incumbent_lhs_sumcomp = - raft::device_span(fj_cpu.h_lhs_sumcomp.data(), fj_cpu.h_lhs_sumcomp.size()); - fj_cpu.view.tabu_nodec_until = - raft::device_span(fj_cpu.h_tabu_nodec_until.data(), fj_cpu.h_tabu_nodec_until.size()); - fj_cpu.view.tabu_noinc_until = - raft::device_span(fj_cpu.h_tabu_noinc_until.data(), fj_cpu.h_tabu_noinc_until.size()); - fj_cpu.view.tabu_lastdec = - raft::device_span(fj_cpu.h_tabu_lastdec.data(), fj_cpu.h_tabu_lastdec.size()); - fj_cpu.view.tabu_lastinc = - raft::device_span(fj_cpu.h_tabu_lastinc.data(), fj_cpu.h_tabu_lastinc.size()); - fj_cpu.view.objective_vars = - raft::device_span(fj_cpu.h_objective_vars.data(), fj_cpu.h_objective_vars.size()); - fj_cpu.view.incumbent_objective = &fj_cpu.h_incumbent_objective; - fj_cpu.view.best_objective = &fj_cpu.h_best_objective; + finalize_fj_cpu_host_initialization(fj_cpu, + problem.n_variables, + problem.n_constraints, + problem.n_integer_vars, + problem.nnz, + problem.tolerances); +} + +template +static void set_host_data_view( + fj_cpu_climber_t& fj_cpu, + i_t n_variables, + i_t n_constraints, + i_t n_integer_vars, + i_t nnz, + const typename mip_solver_settings_t::tolerances_t& tolerances) +{ + fj_cpu.view.pb.tolerances = tolerances; + fj_cpu.view.pb.n_variables = n_variables; + fj_cpu.view.pb.n_integer_vars = n_integer_vars; + fj_cpu.view.pb.n_constraints = n_constraints; + fj_cpu.view.pb.nnz = nnz; - fj_cpu.view.settings = &fj_cpu.settings; fj_cpu.view.pb.constraint_lower_bounds = raft::device_span(fj_cpu.h_cstr_lb.data(), fj_cpu.h_cstr_lb.size()); fj_cpu.view.pb.constraint_upper_bounds = @@ -1295,6 +1343,8 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, raft::device_span(fj_cpu.h_var_types.data(), fj_cpu.h_var_types.size()); fj_cpu.view.pb.is_binary_variable = raft::device_span(fj_cpu.h_is_binary_variable.data(), fj_cpu.h_is_binary_variable.size()); + fj_cpu.view.pb.binary_indices = + raft::device_span(fj_cpu.h_binary_indices.data(), fj_cpu.h_binary_indices.size()); fj_cpu.view.pb.coefficients = raft::device_span(fj_cpu.h_coefficients.data(), fj_cpu.h_coefficients.size()); fj_cpu.view.pb.offsets = raft::device_span(fj_cpu.h_offsets.data(), fj_cpu.h_offsets.size()); @@ -1308,22 +1358,69 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, raft::device_span(fj_cpu.h_reverse_offsets.data(), fj_cpu.h_reverse_offsets.size()); fj_cpu.view.pb.objective_coefficients = raft::device_span(fj_cpu.h_obj_coeffs.data(), fj_cpu.h_obj_coeffs.size()); - fj_cpu.h_objective_vars.resize(problem.n_variables); +} + +template +void finalize_fj_cpu_host_initialization( + fj_cpu_climber_t& fj_cpu, + i_t n_variables, + i_t n_constraints, + i_t n_integer_vars, + i_t nnz, + const typename mip_solver_settings_t::tolerances_t& tolerances) +{ + raft::common::nvtx::range scope("finalize_fj_cpu_host_initialization"); + + cuopt_assert(n_variables >= 0, "invalid variable count"); + cuopt_assert(n_constraints >= 0, "invalid constraint count"); + cuopt_assert(fj_cpu.h_offsets.size() == static_cast(n_constraints + 1), + "invalid CSR offsets"); + cuopt_assert(fj_cpu.h_reverse_offsets.size() == static_cast(n_variables + 1), + "invalid reverse offsets"); + cuopt_assert(fj_cpu.h_assignment.size() == static_cast(n_variables), + "seed assignment size mismatch"); + + set_host_data_view(fj_cpu, n_variables, n_constraints, n_integer_vars, nnz, tolerances); + + fj_cpu.view.cstr_left_weights = + raft::device_span(fj_cpu.h_cstr_left_weights.data(), fj_cpu.h_cstr_left_weights.size()); + fj_cpu.view.cstr_right_weights = + raft::device_span(fj_cpu.h_cstr_right_weights.data(), fj_cpu.h_cstr_right_weights.size()); + fj_cpu.view.objective_weight = &fj_cpu.h_objective_weight; + fj_cpu.view.incumbent_assignment = + raft::device_span(fj_cpu.h_assignment.data(), fj_cpu.h_assignment.size()); + fj_cpu.view.incumbent_lhs = raft::device_span(fj_cpu.h_lhs.data(), fj_cpu.h_lhs.size()); + fj_cpu.view.incumbent_lhs_sumcomp = + raft::device_span(fj_cpu.h_lhs_sumcomp.data(), fj_cpu.h_lhs_sumcomp.size()); + fj_cpu.view.tabu_nodec_until = + raft::device_span(fj_cpu.h_tabu_nodec_until.data(), fj_cpu.h_tabu_nodec_until.size()); + fj_cpu.view.tabu_noinc_until = + raft::device_span(fj_cpu.h_tabu_noinc_until.data(), fj_cpu.h_tabu_noinc_until.size()); + fj_cpu.view.tabu_lastdec = + raft::device_span(fj_cpu.h_tabu_lastdec.data(), fj_cpu.h_tabu_lastdec.size()); + fj_cpu.view.tabu_lastinc = + raft::device_span(fj_cpu.h_tabu_lastinc.data(), fj_cpu.h_tabu_lastinc.size()); + fj_cpu.view.incumbent_objective = &fj_cpu.h_incumbent_objective; + fj_cpu.view.best_objective = &fj_cpu.h_best_objective; + fj_cpu.view.settings = &fj_cpu.settings; + + fj_cpu.h_objective_vars.resize(n_variables); auto end = std::copy_if( thrust::counting_iterator(0), - thrust::counting_iterator(problem.n_variables), + thrust::counting_iterator(n_variables), fj_cpu.h_objective_vars.begin(), [&fj_cpu](i_t idx) { return !fj_cpu.view.pb.integer_equal(fj_cpu.h_obj_coeffs[idx], (f_t)0); }); fj_cpu.h_objective_vars.resize(end - fj_cpu.h_objective_vars.begin()); + fj_cpu.view.objective_vars = + raft::device_span(fj_cpu.h_objective_vars.data(), fj_cpu.h_objective_vars.size()); fj_cpu.h_best_objective = +std::numeric_limits::infinity(); // nnz count fj_cpu.cached_mtm_moves.resize(fj_cpu.h_coefficients.size(), std::make_pair(0, fj_staged_score_t::zero())); - fj_cpu.cached_cstr_bounds.resize(fj_cpu.h_reverse_coefficients.size()); - for (i_t var_idx = 0; var_idx < (i_t)fj_cpu.view.pb.n_variables; ++var_idx) { + for (i_t var_idx = 0; var_idx < n_variables; ++var_idx) { auto [offset_begin, offset_end] = reverse_range_for_var(fj_cpu, var_idx); for (i_t i = offset_begin; i < offset_end; ++i) { fj_cpu.cached_cstr_bounds[i] = @@ -1332,9 +1429,9 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, } } - fj_cpu.flip_move_computed.resize(fj_cpu.view.pb.n_variables, false); - fj_cpu.var_bitmap.resize(fj_cpu.view.pb.n_variables, false); - fj_cpu.iter_mtm_vars.reserve(fj_cpu.view.pb.n_variables); + fj_cpu.flip_move_computed.resize(n_variables, false); + fj_cpu.var_bitmap.resize(n_variables, false); + fj_cpu.iter_mtm_vars.reserve(n_variables); recompute_lhs(fj_cpu); @@ -1342,6 +1439,124 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, precompute_problem_features(fj_cpu); } +template +static std::unique_ptr> init_fj_cpu_from_host_lp( + const dual_simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const dual_simplex::simplex_solver_settings_t& settings, + std::atomic& preemption_flag) +{ + using f_t2 = typename type_2::type; + + cuopt_assert(variable_types.size() >= static_cast(problem.num_cols), + "variable type size mismatch"); + + typename mip_solver_settings_t::tolerances_t tolerances{}; + tolerances.absolute_tolerance = settings.primal_tol; + tolerances.integrality_tolerance = settings.integer_tol; + tolerances.absolute_mip_gap = settings.absolute_mip_gap_tol; + tolerances.relative_mip_gap = settings.relative_mip_gap_tol; + + const i_t n_variables = problem.num_cols; + const i_t n_constraints = problem.num_rows; + + dual_simplex::csr_matrix_t csr_A(problem.num_rows, problem.num_cols, problem.A.nnz()); + problem.A.to_compressed_row(csr_A); + std::vector coefficients = csr_A.x; + std::vector variables = csr_A.j; + std::vector offsets = csr_A.row_start; + std::vector constraint_lower_bounds = problem.rhs; + std::vector constraint_upper_bounds = problem.rhs; + std::vector variable_bounds(n_variables); + std::vector cpufj_variable_types(n_variables); + std::vector is_binary_variable(n_variables, 0); + i_t n_integer_vars = 0; + + for (i_t j = 0; j < n_variables; ++j) { + variable_bounds[j] = f_t2{problem.lower[j], problem.upper[j]}; + const auto var_type = variable_types[j]; + cpufj_variable_types[j] = + var_type == dual_simplex::variable_type_t::CONTINUOUS ? var_t::CONTINUOUS : var_t::INTEGER; + + const bool is_integer = cpufj_variable_types[j] == var_t::INTEGER; + const bool is_binary = is_integer && + integer_equal(problem.lower[j], f_t{0}, settings.integer_tol) && + integer_equal(problem.upper[j], f_t{1}, settings.integer_tol); + if (is_integer) { ++n_integer_vars; } + if (is_binary) { is_binary_variable[j] = 1; } + } + + const i_t nnz = static_cast(variables.size()); + std::vector reverse_coefficients; + std::vector reverse_constraints; + std::vector reverse_offsets; + rebuild_reverse_matrix(n_variables, + n_constraints, + nnz, + coefficients, + variables, + offsets, + reverse_coefficients, + reverse_constraints, + reverse_offsets); + + std::vector projected_seed(n_variables, f_t{0}); + for (i_t j = 0; j < n_variables; ++j) { + f_t value = j < static_cast(seed_assignment.size()) ? seed_assignment[j] : f_t{0}; + value = clamp_value(value, problem.lower[j], problem.upper[j]); + if (variable_types[j] != dual_simplex::variable_type_t::CONTINUOUS) { + value = clamp_value(std::round(value), problem.lower[j], problem.upper[j]); + } + projected_seed[j] = value; + } + + fj_settings_t fj_settings; + fj_settings.mode = fj_mode_t::EXIT_NON_IMPROVING; + fj_settings.n_of_minimums_for_exit = std::numeric_limits::max(); + fj_settings.time_limit = std::numeric_limits::infinity(); + fj_settings.iteration_limit = std::numeric_limits::max(); + fj_settings.update_weights = true; + fj_settings.feasibility_run = false; + fj_settings.seed = cuopt::seed_generator::get_seed(); + + auto fj_cpu = std::make_unique>(preemption_flag); + fj_cpu->view = typename fj_t::climber_data_t::view_t{}; + fj_cpu->pb_ptr = nullptr; + fj_cpu->settings = fj_settings; + + fj_cpu->h_reverse_coefficients = std::move(reverse_coefficients); + fj_cpu->h_reverse_constraints = std::move(reverse_constraints); + fj_cpu->h_reverse_offsets = std::move(reverse_offsets); + fj_cpu->h_coefficients = std::move(coefficients); + fj_cpu->h_offsets = std::move(offsets); + fj_cpu->h_variables = std::move(variables); + fj_cpu->h_obj_coeffs = problem.objective; + fj_cpu->h_var_bounds = std::move(variable_bounds); + fj_cpu->h_cstr_lb = std::move(constraint_lower_bounds); + fj_cpu->h_cstr_ub = std::move(constraint_upper_bounds); + fj_cpu->h_var_types = std::move(cpufj_variable_types); + fj_cpu->h_is_binary_variable = std::move(is_binary_variable); + + fj_cpu->h_cstr_left_weights.resize(n_constraints, 1.0); + fj_cpu->h_cstr_right_weights.resize(n_constraints, 1.0); + fj_cpu->max_weight = 1.0; + fj_cpu->h_objective_weight = 0.0; + fj_cpu->h_assignment = projected_seed; + fj_cpu->h_best_assignment = std::move(projected_seed); + fj_cpu->h_lhs.resize(n_constraints); + fj_cpu->h_lhs_sumcomp.resize(n_constraints, 0); + fj_cpu->h_tabu_nodec_until.resize(n_variables, 0); + fj_cpu->h_tabu_noinc_until.resize(n_variables, 0); + fj_cpu->h_tabu_lastdec.resize(n_variables, 0); + fj_cpu->h_tabu_lastinc.resize(n_variables, 0); + fj_cpu->iterations = 0; + + finalize_fj_cpu_host_initialization( + *fj_cpu, n_variables, n_constraints, n_integer_vars, nnz, tolerances); + return fj_cpu; +} + template static void sanity_checks(fj_cpu_climber_t& fj_cpu) { @@ -1417,7 +1632,9 @@ std::unique_ptr> fj_t::create_cpu_climber( } template -static bool cpufj_solve_loop(fj_cpu_climber_t& fj_cpu, f_t in_time_limit) +static bool cpufj_solve_loop(fj_cpu_climber_t& fj_cpu, + f_t in_time_limit, + double work_unit_limit = std::numeric_limits::infinity()) { i_t local_mins = 0; auto loop_start = std::chrono::high_resolution_clock::now(); @@ -1518,7 +1735,7 @@ static bool cpufj_solve_loop(fj_cpu_climber_t& fj_cpu, f_t in_time_lim fj_cpu.total_violations += fj_cpu.view.excess_score(cstr_idx, fj_cpu.h_lhs[cstr_idx]); } if (fj_cpu.iterations % fj_cpu.log_interval == 0) { - CUOPT_LOG_TRACE( + CUOPT_LOG_DEBUG( "%sCPUFJ iteration: %d/%d, local mins: %d, best_objective: %g, viol: %zu, obj weight %g, " "maxw %g", fj_cpu.log_prefix.c_str(), @@ -1527,7 +1744,7 @@ static bool cpufj_solve_loop(fj_cpu_climber_t& fj_cpu, f_t in_time_lim ? fj_cpu.settings.iteration_limit : -1, local_mins, - fj_cpu.pb_ptr->get_user_obj_from_solver_obj(fj_cpu.h_best_objective), + fj_cpu.h_best_objective, fj_cpu.violated_constraints.size(), fj_cpu.h_objective_weight, fj_cpu.max_weight); @@ -1553,6 +1770,7 @@ static bool cpufj_solve_loop(fj_cpu_climber_t& fj_cpu, f_t in_time_lim fj_cpu.work_units_elapsed += biased_work; if (fj_cpu.producer_sync != nullptr) { fj_cpu.producer_sync->notify_progress(); } + if (fj_cpu.work_units_elapsed.load(std::memory_order_acquire) >= work_unit_limit) { break; } } cuopt_func_call(sanity_checks(fj_cpu)); @@ -1592,7 +1810,7 @@ cpu_fj_thread_t::~cpu_fj_thread_t() template void cpu_fj_thread_t::run_worker() { - cpu_fj_solution_found = cpufj_solve_loop(*fj_cpu, time_limit); + cpu_fj_solution_found = cpufj_solve_loop(*fj_cpu, time_limit, work_unit_limit); } template @@ -1633,6 +1851,29 @@ std::unique_ptr> init_fj_cpu_standalone( return fj_cpu; } +template +bool run_fj_cpu_from_host_lp( + const dual_simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const dual_simplex::simplex_solver_settings_t& settings, + std::atomic& preemption_flag, + f_t time_limit, + double work_unit_limit, + std::function&, double)> improvement_callback, + std::string log_prefix) +{ + cpu_fj_thread_t cpu_fj; + cpu_fj.fj_cpu = + init_fj_cpu_from_host_lp(problem, variable_types, seed_assignment, settings, preemption_flag); + cpu_fj.fj_cpu->log_prefix = std::move(log_prefix); + cpu_fj.fj_cpu->improvement_callback = std::move(improvement_callback); + cpu_fj.time_limit = time_limit; + cpu_fj.work_unit_limit = work_unit_limit; + cpu_fj.start_cpu_solver(); + return cpu_fj.wait_for_cpu_solver(); +} + #if MIP_INSTANTIATE_FLOAT template class fj_t; template class cpu_fj_thread_t; @@ -1641,6 +1882,23 @@ template std::unique_ptr> init_fj_cpu_standalone( solution_t& solution, std::atomic& preemption_flag, fj_settings_t settings); +template bool run_fj_cpu_from_host_lp( + const dual_simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const dual_simplex::simplex_solver_settings_t& settings, + std::atomic& preemption_flag, + float time_limit, + double work_unit_limit, + std::function&, double)> improvement_callback, + std::string log_prefix); +template void finalize_fj_cpu_host_initialization( + fj_cpu_climber_t& fj_cpu, + int n_variables, + int n_constraints, + int n_integer_vars, + int nnz, + const typename mip_solver_settings_t::tolerances_t& tolerances); #endif #if MIP_INSTANTIATE_DOUBLE @@ -1651,6 +1909,23 @@ template std::unique_ptr> init_fj_cpu_standalone( solution_t& solution, std::atomic& preemption_flag, fj_settings_t settings); +template bool run_fj_cpu_from_host_lp( + const dual_simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const dual_simplex::simplex_solver_settings_t& settings, + std::atomic& preemption_flag, + double time_limit, + double work_unit_limit, + std::function&, double)> improvement_callback, + std::string log_prefix); +template void finalize_fj_cpu_host_initialization( + fj_cpu_climber_t& fj_cpu, + int n_variables, + int n_constraints, + int n_integer_vars, + int nnz, + const typename mip_solver_settings_t::tolerances_t& tolerances); #endif } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh index 3263609a2b..8982dfa5f7 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh @@ -11,13 +11,14 @@ #include #include #include +#include #include #include #include #include +#include #include -#include #include #include @@ -193,23 +194,6 @@ struct fj_cpu_climber_t { std::atomic& preemption_flag; }; -template -struct cpu_fj_thread_t : public cpu_worker_thread_base_t> { - ~cpu_fj_thread_t(); - - void run_worker(); - void on_terminate(); - void on_start(); - bool get_result() { return cpu_fj_solution_found; } - - void stop_cpu_solver(); - - std::atomic cpu_fj_solution_found{false}; - f_t time_limit{+std::numeric_limits::infinity()}; - std::unique_ptr> fj_cpu; - fj_t* fj_ptr{nullptr}; -}; - // Standalone CPUFJ init for running without full fj_t infrastructure (avoids GPU allocations). // Used for early CPUFJ during presolve. template From 80ee80a5349a23f4c14c54db39334405511619fe Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Mon, 27 Apr 2026 06:56:53 -0700 Subject: [PATCH 2/6] run cpufj in parallel during cuts --- cpp/src/branch_and_bound/branch_and_bound.cpp | 667 +++++++++++------- cpp/src/branch_and_bound/branch_and_bound.hpp | 27 + .../feasibility_jump/cpu_fj_thread.cuh | 20 +- .../mip_heuristics/feasibility_jump/fj_cpu.cu | 61 +- 4 files changed, 482 insertions(+), 293 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 2ecd01558e..6057f0c299 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -26,6 +26,7 @@ #include #include +#include #include @@ -2016,6 +2017,288 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( return root_status; } +template +typename branch_and_bound_t::root_cut_pass_result_t +branch_and_bound_t::run_root_cut_pass(i_t cut_pass, + mip_solution_t& solution, + cut_generation_t& cut_generation, + cut_pool_t& cut_pool, + variable_bounds_t& variable_bounds, + basis_update_mpf_t& basis_update, + simplex_solver_settings_t& lp_settings, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& edge_norms, + std::vector& fractional, + i_t& num_fractional, + cut_info_t& cut_info, + const std::vector& saved_solution, + f_t& last_upper_bound, + f_t& last_objective, + f_t root_relax_objective, + i_t original_rows, + i_t& cut_pool_size) +{ +#ifdef PRINT_FRACTIONAL_INFO + settings_.log.printf("Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); + for (i_t j : fractional) { + settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", + j, + original_lp_.lower[j], + root_relax_soln_.x[j], + original_lp_.upper[j]); + } +#endif + + f_t cut_start_time = tic(); + bool problem_feasible = cut_generation.generate_cuts(original_lp_, + settings_, + Arow_, + new_slacks_, + var_types_, + basis_update, + root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, + basic_list, + nonbasic_list, + variable_bounds, + exploration_stats_.start_time); + if (!problem_feasible) { + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); + } + if (clique_table_future_.valid()) { + signal_extend_cliques_.store(true, std::memory_order_release); + clique_table_ = clique_table_future_.get(); + } + return {root_cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; + } + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); + } + + f_t score_start_time = tic(); + cut_pool.score_cuts(root_relax_soln_.x); + f_t score_time = toc(score_start_time); + if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } + + csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); + std::vector cut_rhs; + std::vector cut_types; + i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); + if (num_cuts == 0) { return {root_cut_pass_action_t::BREAK, mip_status_t::UNSET}; } + cut_info.record_cut_types(cut_types); +#ifdef PRINT_CUT_POOL_TYPES + cut_pool.print_cutpool_types(); + print_cut_types("In LP ", cut_types, settings_); + printf("Cut pool size: %d\n", cut_pool.pool_size()); +#endif + +#ifdef CHECK_CUT_MATRIX + if (cuts_to_add.check_matrix() != 0) { + settings_.log.printf("Bad cuts matrix\n"); + for (i_t i = 0; i < static_cast(cut_types.size()); ++i) { + settings_.log.printf("row %d cut type %d\n", i, cut_types[i]); + } + return {root_cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; + } +#endif + +#ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION + verify_cuts_against_saved_solution(cuts_to_add, cut_rhs, saved_solution); +#else + (void)saved_solution; +#endif + cut_pool_size = cut_pool.pool_size(); + + settings_.log.debug( + "Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", + num_cuts, + cuts_to_add.row_start[cuts_to_add.m], + cut_pool.pool_size(), + cuts_to_add.m + original_lp_.num_rows); + lp_settings.log.log = false; + + f_t add_cuts_start_time = tic(); + mutex_original_lp_.lock(); + i_t add_cuts_status = add_cuts(settings_, + cuts_to_add, + cut_rhs, + original_lp_, + new_slacks_, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms); + var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); + variable_bounds.resize(original_lp_.num_cols); + mutex_original_lp_.unlock(); + f_t add_cuts_time = toc(add_cuts_start_time); + if (add_cuts_time > 1.0) { settings_.log.debug("Add cuts time %.2f seconds\n", add_cuts_time); } + if (add_cuts_status != 0) { + settings_.log.printf("Failed to add cuts\n"); + return {root_cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; + } + + if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) { + mutex_upper_.lock(); + last_upper_bound = upper_bound_.load(); + std::vector lower_bounds; + std::vector upper_bounds; + find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); + mutex_upper_.unlock(); + mutex_original_lp_.lock(); + original_lp_.lower = lower_bounds; + original_lp_.upper = upper_bounds; + mutex_original_lp_.unlock(); + } + + std::vector bounds_changed(original_lp_.num_cols, true); + std::vector row_sense; +#ifdef CHECK_MATRICES + settings_.log.printf("Before A check\n"); + original_lp_.A.check_matrix(); +#endif + original_lp_.A.to_compressed_row(Arow_); + + f_t node_presolve_start_time = tic(); + bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); + std::vector new_lower = original_lp_.lower; + std::vector new_upper = original_lp_.upper; + bool feasible = + node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); + mutex_original_lp_.lock(); + original_lp_.lower = new_lower; + original_lp_.upper = new_upper; + mutex_original_lp_.unlock(); + f_t node_presolve_time = toc(node_presolve_start_time); + if (node_presolve_time > 1.0) { + settings_.log.debug("Node presolve time %.2f seconds\n", node_presolve_time); + } + if (!feasible) { + settings_.log.printf("Bound strengthening detected infeasibility\n"); +#ifdef WRITE_BOUND_STRENGTHENING_INFEASIBLE_MPS + original_lp_.write_mps("bound_strengthening_infeasible.mps"); +#endif + return {root_cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; + } + + i_t iter = 0; + bool initialize_basis = false; + lp_settings.concurrent_halt = NULL; + f_t dual_phase2_start_time = tic(); + dual::status_t cut_status = dual_phase2_with_advanced_basis(2, + 0, + initialize_basis, + exploration_stats_.start_time, + original_lp_, + lp_settings, + root_vstatus_, + basis_update, + basic_list, + nonbasic_list, + root_relax_soln_, + iter, + edge_norms); + exploration_stats_.total_lp_iters += iter; + f_t dual_phase2_time = toc(dual_phase2_start_time); + if (dual_phase2_time > 1.0) { + settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); + } + if (cut_status == dual::status_t::TIME_LIMIT) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, root_objective_); + return {root_cut_pass_action_t::RETURN, solver_status_}; + } + + if (cut_status != dual::status_t::OPTIMAL) { + settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); + lp_status_t scratch_status = + solve_linear_program_with_advanced_basis(original_lp_, + exploration_stats_.start_time, + lp_settings, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms); + if (scratch_status == lp_status_t::OPTIMAL) { + cut_status = convert_lp_status_to_dual_status(scratch_status); + exploration_stats_.total_lp_iters += root_relax_soln_.iterations; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + } else { + settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); +#ifdef WRITE_CUT_INFEASIBLE_MPS + original_lp_.write_mps("cut_infeasible.mps"); +#endif + return {root_cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; + } + } + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + + f_t remove_cuts_start_time = tic(); + mutex_original_lp_.lock(); + remove_cuts(original_lp_, + settings_, + exploration_stats_.start_time, + Arow_, + new_slacks_, + original_rows, + var_types_, + root_vstatus_, + edge_norms, + root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, + basic_list, + nonbasic_list, + basis_update); + variable_bounds.resize(original_lp_.num_cols); + mutex_original_lp_.unlock(); + f_t remove_cuts_time = toc(remove_cuts_start_time); + if (remove_cuts_time > 1.0) { + settings_.log.debug("Remove cuts time %.2f seconds\n", remove_cuts_time); + } + fractional.clear(); + num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + + if (num_fractional == 0) { + upper_bound_ = root_objective_; + mutex_upper_.lock(); + incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); + mutex_upper_.unlock(); + } + f_t obj = upper_bound_.load(); + report(' ', obj, root_objective_, 0, num_fractional); + + f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); + f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound_.load(), root_objective_); + if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { + if (num_fractional == 0) { set_solution_at_root(solution, cut_info); } + set_final_solution(solution, root_objective_); + return {root_cut_pass_action_t::RETURN, mip_status_t::OPTIMAL}; + } + + f_t change_in_objective = root_objective_ - last_objective; + const f_t factor = settings_.cut_change_threshold; + const f_t min_objective = 1e-3; + if (factor > 0.0 && + change_in_objective <= factor * std::max(min_objective, std::abs(root_relax_objective))) { + settings_.log.printf( + "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", + change_in_objective, + root_relax_objective); + return {root_cut_pass_action_t::BREAK, mip_status_t::UNSET}; + } + last_objective = root_objective_; + return {root_cut_pass_action_t::CONTINUE, mip_status_t::UNSET}; +} + template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { @@ -2228,6 +2511,22 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut constexpr bool enable_root_cut_cpufj = true; constexpr double root_cut_cpufj_work_units = 0.1; + std::unique_ptr> root_cut_cpufj_task; + std::atomic root_cut_cpufj_preemption{false}; + auto root_cut_cpufj_improvement_callback = + [this](f_t obj, const std::vector& assignment, double) { + std::vector user_assignment(assignment.begin(), + assignment.begin() + original_problem_.num_cols); + CUOPT_LOG_INFO("Root cut CPUFJ found solution with objective %.16e\n", obj); + set_new_solution(user_assignment); + }; + auto stop_root_cut_cpufj = [&]() { + if (!root_cut_cpufj_task) { return; } + CUOPT_LOG_DEBUG("Stopping CPUFJ for this cut pass"); + detail::stop_fj_cpu_task(*root_cut_cpufj_task); + root_cut_cpufj_task.reset(); + }; + cuopt::scope_guard root_cut_cpufj_guard([&]() { stop_root_cut_cpufj(); }); f_t cut_generation_start_time = tic(); i_t cut_pool_size = 0; @@ -2235,266 +2534,90 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut if (num_fractional == 0) { set_solution_at_root(solution, cut_info); return mip_status_t::OPTIMAL; - } else { -#ifdef PRINT_FRACTIONAL_INFO - settings_.log.printf( - "Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); - for (i_t j : fractional) { - settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", - j, - original_lp_.lower[j], - root_relax_soln_.x[j], - original_lp_.upper[j]); - } -#endif - - // Generate cuts and add them to the cut pool - f_t cut_start_time = tic(); - bool problem_feasible = cut_generation.generate_cuts(original_lp_, - settings_, - Arow_, - new_slacks_, - var_types_, - basis_update, - root_relax_soln_.x, - root_relax_soln_.y, - root_relax_soln_.z, - basic_list, - nonbasic_list, - variable_bounds, - exploration_stats_.start_time); - if (!problem_feasible) { - if (settings_.heuristic_preemption_callback != nullptr) { - settings_.heuristic_preemption_callback(); - } - finish_clique_thread(); - return mip_status_t::INFEASIBLE; - } - f_t cut_generation_time = toc(cut_start_time); - if (cut_generation_time > 1.0) { - settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); - } - // Score the cuts - f_t score_start_time = tic(); - cut_pool.score_cuts(root_relax_soln_.x); - f_t score_time = toc(score_start_time); - if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } - // Get the best cuts from the cut pool - csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); - std::vector cut_rhs; - std::vector cut_types; - i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); - if (num_cuts == 0) { break; } - cut_info.record_cut_types(cut_types); -#ifdef PRINT_CUT_POOL_TYPES - cut_pool.print_cutpool_types(); - print_cut_types("In LP ", cut_types, settings_); - printf("Cut pool size: %d\n", cut_pool.pool_size()); -#endif - -#ifdef CHECK_CUT_MATRIX - if (cuts_to_add.check_matrix() != 0) { - settings_.log.printf("Bad cuts matrix\n"); - for (i_t i = 0; i < static_cast(cut_types.size()); ++i) { - settings_.log.printf("row %d cut type %d\n", i, cut_types[i]); - } - return mip_status_t::NUMERICAL; - } -#endif - // Check against saved solution -#ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION - verify_cuts_against_saved_solution(cuts_to_add, cut_rhs, saved_solution); -#endif - cut_pool_size = cut_pool.pool_size(); - - // Resolve the LP with the new cuts - settings_.log.debug( - "Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", - num_cuts, - cuts_to_add.row_start[cuts_to_add.m], - cut_pool.pool_size(), - cuts_to_add.m + original_lp_.num_rows); - lp_settings.log.log = false; - - f_t add_cuts_start_time = tic(); - mutex_original_lp_.lock(); - i_t add_cuts_status = add_cuts(settings_, - cuts_to_add, - cut_rhs, - original_lp_, - new_slacks_, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms_); - var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); - variable_bounds.resize(original_lp_.num_cols); - mutex_original_lp_.unlock(); - f_t add_cuts_time = toc(add_cuts_start_time); - if (add_cuts_time > 1.0) { - settings_.log.debug("Add cuts time %.2f seconds\n", add_cuts_time); - } - if (add_cuts_status != 0) { - settings_.log.printf("Failed to add cuts\n"); - return mip_status_t::NUMERICAL; - } - - if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) { - mutex_upper_.lock(); - last_upper_bound = upper_bound_.load(); - std::vector lower_bounds; - std::vector upper_bounds; - find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); - mutex_upper_.unlock(); - mutex_original_lp_.lock(); - original_lp_.lower = lower_bounds; - original_lp_.upper = upper_bounds; - mutex_original_lp_.unlock(); - } - - // Try to do bound strengthening - std::vector bounds_changed(original_lp_.num_cols, true); - std::vector row_sense; -#ifdef CHECK_MATRICES - settings_.log.printf("Before A check\n"); - original_lp_.A.check_matrix(); -#endif - original_lp_.A.to_compressed_row(Arow_); - - f_t node_presolve_start_time = tic(); - bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); - std::vector new_lower = original_lp_.lower; - std::vector new_upper = original_lp_.upper; - bool feasible = - node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); - mutex_original_lp_.lock(); - original_lp_.lower = new_lower; - original_lp_.upper = new_upper; - mutex_original_lp_.unlock(); - f_t node_presolve_time = toc(node_presolve_start_time); - if (node_presolve_time > 1.0) { - settings_.log.debug("Node presolve time %.2f seconds\n", node_presolve_time); - } - if (!feasible) { - settings_.log.printf("Bound strengthening detected infeasibility\n"); -#ifdef WRITE_BOUND_STRENGTHENING_INFEASIBLE_MPS - original_lp_.write_mps("bound_strengthening_infeasible.mps"); -#endif - return mip_status_t::INFEASIBLE; - } - - i_t iter = 0; - bool initialize_basis = false; - lp_settings.concurrent_halt = NULL; - f_t dual_phase2_start_time = tic(); - dual::status_t cut_status = dual_phase2_with_advanced_basis(2, - 0, - initialize_basis, - exploration_stats_.start_time, - original_lp_, - lp_settings, - root_vstatus_, - basis_update, - basic_list, - nonbasic_list, - root_relax_soln_, - iter, - edge_norms_); - exploration_stats_.total_lp_iters += iter; - f_t dual_phase2_time = toc(dual_phase2_start_time); - if (dual_phase2_time > 1.0) { - settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); - } - if (cut_status == dual::status_t::TIME_LIMIT) { - solver_status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, root_objective_); - return solver_status_; - } + } - if (cut_status != dual::status_t::OPTIMAL) { - settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); - lp_status_t scratch_status = - solve_linear_program_with_advanced_basis(original_lp_, - exploration_stats_.start_time, - lp_settings, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms_); - if (scratch_status == lp_status_t::OPTIMAL) { - // We recovered - cut_status = convert_lp_status_to_dual_status(scratch_status); - exploration_stats_.total_lp_iters += root_relax_soln_.iterations; - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - } else { - settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); -#ifdef WRITE_CUT_INFEASIBLE_MPS - original_lp_.write_mps("cut_infeasible.mps"); -#endif - return mip_status_t::NUMERICAL; + root_cut_pass_result_t cut_pass_result; + if (root_cut_cpufj_task) { + f_t root_cut_cpufj_stop_start_time = 0.0; +#pragma omp parallel num_threads(settings_.num_threads) + { +#pragma omp single + { +#pragma omp taskgroup + { +#pragma omp task shared(root_cut_cpufj_task) firstprivate(root_cut_cpufj_work_units) default(none) + detail::run_fj_cpu_task(*root_cut_cpufj_task, + std::numeric_limits::infinity(), + std::numeric_limits::infinity()); + + cut_pass_result = run_root_cut_pass(cut_pass, + solution, + cut_generation, + cut_pool, + variable_bounds, + basis_update, + lp_settings, + basic_list, + nonbasic_list, + edge_norms_, + fractional, + num_fractional, + cut_info, + saved_solution, + last_upper_bound, + last_objective, + root_relax_objective, + original_rows, + cut_pool_size); + + root_cut_cpufj_stop_start_time = tic(); + detail::stop_fj_cpu_task(*root_cut_cpufj_task); + } + settings_.log.printf("Root cut CPUFJ stop wait time: %.6f seconds\n", + toc(root_cut_cpufj_stop_start_time)); } } - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - - f_t remove_cuts_start_time = tic(); - mutex_original_lp_.lock(); - remove_cuts(original_lp_, - settings_, - exploration_stats_.start_time, - Arow_, - new_slacks_, - original_rows, - var_types_, - root_vstatus_, - edge_norms_, - root_relax_soln_.x, - root_relax_soln_.y, - root_relax_soln_.z, - basic_list, - nonbasic_list, - basis_update); - variable_bounds.resize(original_lp_.num_cols); - mutex_original_lp_.unlock(); - f_t remove_cuts_time = toc(remove_cuts_start_time); - if (remove_cuts_time > 1.0) { - settings_.log.debug("Remove cuts time %.2f seconds\n", remove_cuts_time); - } - fractional.clear(); - num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); - - if (num_fractional == 0) { - upper_bound_ = root_objective_; - mutex_upper_.lock(); - incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); - mutex_upper_.unlock(); - } - f_t obj = upper_bound_.load(); - report(' ', obj, root_objective_, 0, num_fractional); - - f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); - f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound_.load(), root_objective_); - if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { - if (num_fractional == 0) { set_solution_at_root(solution, cut_info); } - set_final_solution(solution, root_objective_); - return mip_status_t::OPTIMAL; - } - - f_t change_in_objective = root_objective_ - last_objective; - const f_t factor = settings_.cut_change_threshold; - const f_t min_objective = 1e-3; - if (factor > 0.0 && - change_in_objective <= factor * std::max(min_objective, std::abs(root_relax_objective))) { - settings_.log.printf( - "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", - change_in_objective, - root_relax_objective); - break; - } - last_objective = root_objective_; + root_cut_cpufj_task.reset(); + } else { + cut_pass_result = run_root_cut_pass(cut_pass, + solution, + cut_generation, + cut_pool, + variable_bounds, + basis_update, + lp_settings, + basic_list, + nonbasic_list, + edge_norms_, + fractional, + num_fractional, + cut_info, + saved_solution, + last_upper_bound, + last_objective, + root_relax_objective, + original_rows, + cut_pool_size); + } + + if (cut_pass_result.action == root_cut_pass_action_t::RETURN) { return cut_pass_result.status; } + if (cut_pass_result.action == root_cut_pass_action_t::BREAK) { break; } + + if (enable_root_cut_cpufj && !settings_.deterministic && settings_.num_threads >= 2 && + cut_pass + 1 < settings_.max_cut_passes && original_lp_.num_rows > original_rows) { + root_cut_cpufj_preemption = false; + f_t root_cut_cpufj_build_start_time = tic(); + root_cut_cpufj_task = + detail::make_fj_cpu_task_from_host_lp(original_lp_, + var_types_, + root_relax_soln_.x, + settings_, + root_cut_cpufj_preemption, + root_cut_cpufj_improvement_callback, + "[RootCut CPUFJ] "); + settings_.log.printf("Root cut CPUFJ problem build time after pass %d: %.6f seconds\n", + cut_pass, + toc(root_cut_cpufj_build_start_time)); } } @@ -2510,22 +2633,20 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (enable_root_cut_cpufj && original_lp_.num_rows > original_rows) { - std::atomic root_cut_cpufj_preemption{false}; - detail::run_fj_cpu_from_host_lp( - original_lp_, - var_types_, - root_relax_soln_.x, - settings_, - root_cut_cpufj_preemption, - f_t{1}, - 1, - [this](f_t, const std::vector& assignment, double obj) { - std::vector user_assignment(assignment.begin(), - assignment.begin() + original_problem_.num_cols); - CUOPT_LOG_INFO("Root cut CPUFJ found solution with objective %.16e\n", obj); - set_new_solution(user_assignment); - }, - "[RootCut CPUFJ] "); + root_cut_cpufj_preemption = false; + f_t root_cut_cpufj_build_start_time = tic(); + root_cut_cpufj_task = + detail::make_fj_cpu_task_from_host_lp(original_lp_, + var_types_, + root_relax_soln_.x, + settings_, + root_cut_cpufj_preemption, + root_cut_cpufj_improvement_callback, + "[RootCut CPUFJ] "); + settings_.log.printf("Root cut CPUFJ final problem build time: %.6f seconds\n", + toc(root_cut_cpufj_build_start_time)); + detail::run_fj_cpu_task(*root_cut_cpufj_task, f_t{1}, f_t{1}); + root_cut_cpufj_task.reset(); } set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index f2917ba930..e69798bb58 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -158,6 +158,33 @@ class branch_and_bound_t { producer_sync_t& get_producer_sync() { return producer_sync_; } private: + enum class root_cut_pass_action_t { CONTINUE, BREAK, RETURN }; + + struct root_cut_pass_result_t { + root_cut_pass_action_t action{root_cut_pass_action_t::CONTINUE}; + mip_status_t status{mip_status_t::UNSET}; + }; + + root_cut_pass_result_t run_root_cut_pass(i_t cut_pass, + mip_solution_t& solution, + cut_generation_t& cut_generation, + cut_pool_t& cut_pool, + variable_bounds_t& variable_bounds, + basis_update_mpf_t& basis_update, + simplex_solver_settings_t& lp_settings, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& edge_norms, + std::vector& fractional, + i_t& num_fractional, + cut_info_t& cut_info, + const std::vector& saved_solution, + f_t& last_upper_bound, + f_t& last_objective, + f_t root_relax_objective, + i_t original_rows, + i_t& cut_pool_size); + const user_problem_t& original_problem_; const simplex_solver_settings_t settings_; const probing_implied_bound_t& probing_implied_bound_; diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh b/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh index 8de191cc71..a041957ece 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh @@ -45,15 +45,29 @@ struct cpu_fj_thread_t : public cpu_worker_thread_base_t -bool run_fj_cpu_from_host_lp( +struct fj_cpu_task_t { + struct fj_cpu_deleter_t { + void operator()(fj_cpu_climber_t* ptr) const; + }; + std::unique_ptr, fj_cpu_deleter_t> fj_cpu; +}; + +template +std::unique_ptr> make_fj_cpu_task_from_host_lp( const dual_simplex::lp_problem_t& problem, const std::vector& variable_types, const std::vector& seed_assignment, const dual_simplex::simplex_solver_settings_t& settings, std::atomic& preemption_flag, - f_t time_limit, - double work_unit_limit, std::function&, double)> improvement_callback, std::string log_prefix); +template +void run_fj_cpu_task(fj_cpu_task_t& task, + f_t time_limit = std::numeric_limits::infinity(), + double work_unit_limit = std::numeric_limits::infinity()); + +template +void stop_fj_cpu_task(fj_cpu_task_t& task); + } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 9175e4291c..6dde7e1686 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -1852,46 +1852,70 @@ std::unique_ptr> init_fj_cpu_standalone( } template -bool run_fj_cpu_from_host_lp( +void fj_cpu_task_t::fj_cpu_deleter_t::operator()(fj_cpu_climber_t* ptr) const +{ + delete ptr; +} + +template +std::unique_ptr> make_fj_cpu_task_from_host_lp( const dual_simplex::lp_problem_t& problem, const std::vector& variable_types, const std::vector& seed_assignment, const dual_simplex::simplex_solver_settings_t& settings, std::atomic& preemption_flag, - f_t time_limit, - double work_unit_limit, std::function&, double)> improvement_callback, std::string log_prefix) { - cpu_fj_thread_t cpu_fj; - cpu_fj.fj_cpu = + auto task = std::make_unique>(); + auto fj_cpu = init_fj_cpu_from_host_lp(problem, variable_types, seed_assignment, settings, preemption_flag); - cpu_fj.fj_cpu->log_prefix = std::move(log_prefix); - cpu_fj.fj_cpu->improvement_callback = std::move(improvement_callback); - cpu_fj.time_limit = time_limit; - cpu_fj.work_unit_limit = work_unit_limit; - cpu_fj.start_cpu_solver(); - return cpu_fj.wait_for_cpu_solver(); + fj_cpu->log_prefix = std::move(log_prefix); + fj_cpu->improvement_callback = std::move(improvement_callback); + task->fj_cpu.reset(fj_cpu.release()); + return task; +} + +template +void run_fj_cpu_task(fj_cpu_task_t& task, f_t time_limit, double work_unit_limit) +{ + cuopt_assert(task.fj_cpu != nullptr, "CPUFJ task has no climber"); + auto& fj_cpu = *task.fj_cpu; + fj_cpu.halted = false; + cpufj_solve_loop(fj_cpu, time_limit, work_unit_limit); +} + +template +void stop_fj_cpu_task(fj_cpu_task_t& task) +{ + if (task.fj_cpu) { + auto& fj_cpu = *task.fj_cpu; + fj_cpu.preemption_flag.store(true, std::memory_order_release); + fj_cpu.halted = true; + } } #if MIP_INSTANTIATE_FLOAT template class fj_t; template class cpu_fj_thread_t; +template struct fj_cpu_task_t; template std::unique_ptr> init_fj_cpu_standalone( problem_t& problem, solution_t& solution, std::atomic& preemption_flag, fj_settings_t settings); -template bool run_fj_cpu_from_host_lp( +template std::unique_ptr> make_fj_cpu_task_from_host_lp( const dual_simplex::lp_problem_t& problem, const std::vector& variable_types, const std::vector& seed_assignment, const dual_simplex::simplex_solver_settings_t& settings, std::atomic& preemption_flag, - float time_limit, - double work_unit_limit, std::function&, double)> improvement_callback, std::string log_prefix); +template void run_fj_cpu_task(fj_cpu_task_t& task, + float time_limit, + double work_unit_limit); +template void stop_fj_cpu_task(fj_cpu_task_t& task); template void finalize_fj_cpu_host_initialization( fj_cpu_climber_t& fj_cpu, int n_variables, @@ -1904,21 +1928,24 @@ template void finalize_fj_cpu_host_initialization( #if MIP_INSTANTIATE_DOUBLE template class fj_t; template class cpu_fj_thread_t; +template struct fj_cpu_task_t; template std::unique_ptr> init_fj_cpu_standalone( problem_t& problem, solution_t& solution, std::atomic& preemption_flag, fj_settings_t settings); -template bool run_fj_cpu_from_host_lp( +template std::unique_ptr> make_fj_cpu_task_from_host_lp( const dual_simplex::lp_problem_t& problem, const std::vector& variable_types, const std::vector& seed_assignment, const dual_simplex::simplex_solver_settings_t& settings, std::atomic& preemption_flag, - double time_limit, - double work_unit_limit, std::function&, double)> improvement_callback, std::string log_prefix); +template void run_fj_cpu_task(fj_cpu_task_t& task, + double time_limit, + double work_unit_limit); +template void stop_fj_cpu_task(fj_cpu_task_t& task); template void finalize_fj_cpu_host_initialization( fj_cpu_climber_t& fj_cpu, int n_variables, From 9e464bb335c1d871b53eb54df3bc0d6b574e86f8 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Thu, 30 Apr 2026 02:02:41 -0700 Subject: [PATCH 3/6] adjust logs, fix time limit --- cpp/src/branch_and_bound/branch_and_bound.cpp | 21 +++++++++++-------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 6057f0c299..a8347ea8b4 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2517,7 +2517,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut [this](f_t obj, const std::vector& assignment, double) { std::vector user_assignment(assignment.begin(), assignment.begin() + original_problem_.num_cols); - CUOPT_LOG_INFO("Root cut CPUFJ found solution with objective %.16e\n", obj); + CUOPT_LOG_DEBUG("Root cut CPUFJ found solution with objective %.16e\n", obj); set_new_solution(user_assignment); }; auto stop_root_cut_cpufj = [&]() { @@ -2573,8 +2573,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_cut_cpufj_stop_start_time = tic(); detail::stop_fj_cpu_task(*root_cut_cpufj_task); } - settings_.log.printf("Root cut CPUFJ stop wait time: %.6f seconds\n", - toc(root_cut_cpufj_stop_start_time)); + settings_.log.debug("Root cut CPUFJ stop wait time: %.6f seconds\n", + toc(root_cut_cpufj_stop_start_time)); } } root_cut_cpufj_task.reset(); @@ -2615,9 +2615,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_cut_cpufj_preemption, root_cut_cpufj_improvement_callback, "[RootCut CPUFJ] "); - settings_.log.printf("Root cut CPUFJ problem build time after pass %d: %.6f seconds\n", - cut_pass, - toc(root_cut_cpufj_build_start_time)); + settings_.log.debug("Root cut CPUFJ problem build time after pass %d: %.6f seconds\n", + cut_pass, + toc(root_cut_cpufj_build_start_time)); } } @@ -2643,9 +2643,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_cut_cpufj_preemption, root_cut_cpufj_improvement_callback, "[RootCut CPUFJ] "); - settings_.log.printf("Root cut CPUFJ final problem build time: %.6f seconds\n", - toc(root_cut_cpufj_build_start_time)); - detail::run_fj_cpu_task(*root_cut_cpufj_task, f_t{1}, f_t{1}); + settings_.log.debug("Root cut CPUFJ final problem build time: %.6f seconds\n", + toc(root_cut_cpufj_build_start_time)); + f_t fj_time_limit = settings_.deterministic + ? f_t(settings_.time_limit - toc(exploration_stats_.start_time)) + : f_t{1}; + detail::run_fj_cpu_task(*root_cut_cpufj_task, fj_time_limit, 0.5); root_cut_cpufj_task.reset(); } From e2300a4510c9a4d6abcb9d14b6fe743fbb244b9f Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Thu, 30 Apr 2026 03:43:08 -0700 Subject: [PATCH 4/6] cleanup --- cpp/src/branch_and_bound/branch_and_bound.cpp | 645 ++++++++---------- cpp/src/branch_and_bound/branch_and_bound.hpp | 27 - .../feasibility_jump/cpu_fj_thread.cuh | 2 +- .../mip_heuristics/feasibility_jump/fj_cpu.cu | 11 +- .../feasibility_jump/fj_cpu.cuh | 1 - .../local_search/local_search.cu | 2 +- 6 files changed, 295 insertions(+), 393 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index a8347ea8b4..f56d0d3e79 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -31,7 +31,6 @@ #include #include -#include #include #include #include @@ -2017,288 +2016,6 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( return root_status; } -template -typename branch_and_bound_t::root_cut_pass_result_t -branch_and_bound_t::run_root_cut_pass(i_t cut_pass, - mip_solution_t& solution, - cut_generation_t& cut_generation, - cut_pool_t& cut_pool, - variable_bounds_t& variable_bounds, - basis_update_mpf_t& basis_update, - simplex_solver_settings_t& lp_settings, - std::vector& basic_list, - std::vector& nonbasic_list, - std::vector& edge_norms, - std::vector& fractional, - i_t& num_fractional, - cut_info_t& cut_info, - const std::vector& saved_solution, - f_t& last_upper_bound, - f_t& last_objective, - f_t root_relax_objective, - i_t original_rows, - i_t& cut_pool_size) -{ -#ifdef PRINT_FRACTIONAL_INFO - settings_.log.printf("Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); - for (i_t j : fractional) { - settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", - j, - original_lp_.lower[j], - root_relax_soln_.x[j], - original_lp_.upper[j]); - } -#endif - - f_t cut_start_time = tic(); - bool problem_feasible = cut_generation.generate_cuts(original_lp_, - settings_, - Arow_, - new_slacks_, - var_types_, - basis_update, - root_relax_soln_.x, - root_relax_soln_.y, - root_relax_soln_.z, - basic_list, - nonbasic_list, - variable_bounds, - exploration_stats_.start_time); - if (!problem_feasible) { - if (settings_.heuristic_preemption_callback != nullptr) { - settings_.heuristic_preemption_callback(); - } - if (clique_table_future_.valid()) { - signal_extend_cliques_.store(true, std::memory_order_release); - clique_table_ = clique_table_future_.get(); - } - return {root_cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; - } - f_t cut_generation_time = toc(cut_start_time); - if (cut_generation_time > 1.0) { - settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); - } - - f_t score_start_time = tic(); - cut_pool.score_cuts(root_relax_soln_.x); - f_t score_time = toc(score_start_time); - if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } - - csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); - std::vector cut_rhs; - std::vector cut_types; - i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); - if (num_cuts == 0) { return {root_cut_pass_action_t::BREAK, mip_status_t::UNSET}; } - cut_info.record_cut_types(cut_types); -#ifdef PRINT_CUT_POOL_TYPES - cut_pool.print_cutpool_types(); - print_cut_types("In LP ", cut_types, settings_); - printf("Cut pool size: %d\n", cut_pool.pool_size()); -#endif - -#ifdef CHECK_CUT_MATRIX - if (cuts_to_add.check_matrix() != 0) { - settings_.log.printf("Bad cuts matrix\n"); - for (i_t i = 0; i < static_cast(cut_types.size()); ++i) { - settings_.log.printf("row %d cut type %d\n", i, cut_types[i]); - } - return {root_cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; - } -#endif - -#ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION - verify_cuts_against_saved_solution(cuts_to_add, cut_rhs, saved_solution); -#else - (void)saved_solution; -#endif - cut_pool_size = cut_pool.pool_size(); - - settings_.log.debug( - "Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", - num_cuts, - cuts_to_add.row_start[cuts_to_add.m], - cut_pool.pool_size(), - cuts_to_add.m + original_lp_.num_rows); - lp_settings.log.log = false; - - f_t add_cuts_start_time = tic(); - mutex_original_lp_.lock(); - i_t add_cuts_status = add_cuts(settings_, - cuts_to_add, - cut_rhs, - original_lp_, - new_slacks_, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms); - var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); - variable_bounds.resize(original_lp_.num_cols); - mutex_original_lp_.unlock(); - f_t add_cuts_time = toc(add_cuts_start_time); - if (add_cuts_time > 1.0) { settings_.log.debug("Add cuts time %.2f seconds\n", add_cuts_time); } - if (add_cuts_status != 0) { - settings_.log.printf("Failed to add cuts\n"); - return {root_cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; - } - - if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) { - mutex_upper_.lock(); - last_upper_bound = upper_bound_.load(); - std::vector lower_bounds; - std::vector upper_bounds; - find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); - mutex_upper_.unlock(); - mutex_original_lp_.lock(); - original_lp_.lower = lower_bounds; - original_lp_.upper = upper_bounds; - mutex_original_lp_.unlock(); - } - - std::vector bounds_changed(original_lp_.num_cols, true); - std::vector row_sense; -#ifdef CHECK_MATRICES - settings_.log.printf("Before A check\n"); - original_lp_.A.check_matrix(); -#endif - original_lp_.A.to_compressed_row(Arow_); - - f_t node_presolve_start_time = tic(); - bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); - std::vector new_lower = original_lp_.lower; - std::vector new_upper = original_lp_.upper; - bool feasible = - node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); - mutex_original_lp_.lock(); - original_lp_.lower = new_lower; - original_lp_.upper = new_upper; - mutex_original_lp_.unlock(); - f_t node_presolve_time = toc(node_presolve_start_time); - if (node_presolve_time > 1.0) { - settings_.log.debug("Node presolve time %.2f seconds\n", node_presolve_time); - } - if (!feasible) { - settings_.log.printf("Bound strengthening detected infeasibility\n"); -#ifdef WRITE_BOUND_STRENGTHENING_INFEASIBLE_MPS - original_lp_.write_mps("bound_strengthening_infeasible.mps"); -#endif - return {root_cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; - } - - i_t iter = 0; - bool initialize_basis = false; - lp_settings.concurrent_halt = NULL; - f_t dual_phase2_start_time = tic(); - dual::status_t cut_status = dual_phase2_with_advanced_basis(2, - 0, - initialize_basis, - exploration_stats_.start_time, - original_lp_, - lp_settings, - root_vstatus_, - basis_update, - basic_list, - nonbasic_list, - root_relax_soln_, - iter, - edge_norms); - exploration_stats_.total_lp_iters += iter; - f_t dual_phase2_time = toc(dual_phase2_start_time); - if (dual_phase2_time > 1.0) { - settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); - } - if (cut_status == dual::status_t::TIME_LIMIT) { - solver_status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, root_objective_); - return {root_cut_pass_action_t::RETURN, solver_status_}; - } - - if (cut_status != dual::status_t::OPTIMAL) { - settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); - lp_status_t scratch_status = - solve_linear_program_with_advanced_basis(original_lp_, - exploration_stats_.start_time, - lp_settings, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms); - if (scratch_status == lp_status_t::OPTIMAL) { - cut_status = convert_lp_status_to_dual_status(scratch_status); - exploration_stats_.total_lp_iters += root_relax_soln_.iterations; - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - } else { - settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); -#ifdef WRITE_CUT_INFEASIBLE_MPS - original_lp_.write_mps("cut_infeasible.mps"); -#endif - return {root_cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; - } - } - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - - f_t remove_cuts_start_time = tic(); - mutex_original_lp_.lock(); - remove_cuts(original_lp_, - settings_, - exploration_stats_.start_time, - Arow_, - new_slacks_, - original_rows, - var_types_, - root_vstatus_, - edge_norms, - root_relax_soln_.x, - root_relax_soln_.y, - root_relax_soln_.z, - basic_list, - nonbasic_list, - basis_update); - variable_bounds.resize(original_lp_.num_cols); - mutex_original_lp_.unlock(); - f_t remove_cuts_time = toc(remove_cuts_start_time); - if (remove_cuts_time > 1.0) { - settings_.log.debug("Remove cuts time %.2f seconds\n", remove_cuts_time); - } - fractional.clear(); - num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); - - if (num_fractional == 0) { - upper_bound_ = root_objective_; - mutex_upper_.lock(); - incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); - mutex_upper_.unlock(); - } - f_t obj = upper_bound_.load(); - report(' ', obj, root_objective_, 0, num_fractional); - - f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); - f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound_.load(), root_objective_); - if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { - if (num_fractional == 0) { set_solution_at_root(solution, cut_info); } - set_final_solution(solution, root_objective_); - return {root_cut_pass_action_t::RETURN, mip_status_t::OPTIMAL}; - } - - f_t change_in_objective = root_objective_ - last_objective; - const f_t factor = settings_.cut_change_threshold; - const f_t min_objective = 1e-3; - if (factor > 0.0 && - change_in_objective <= factor * std::max(min_objective, std::abs(root_relax_objective))) { - settings_.log.printf( - "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", - change_in_objective, - root_relax_objective); - return {root_cut_pass_action_t::BREAK, mip_status_t::UNSET}; - } - last_objective = root_objective_; - return {root_cut_pass_action_t::CONTINUE, mip_status_t::UNSET}; -} - template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { @@ -2509,10 +2226,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t last_objective = root_objective_; f_t root_relax_objective = root_objective_; - constexpr bool enable_root_cut_cpufj = true; - constexpr double root_cut_cpufj_work_units = 0.1; + constexpr bool enable_root_cut_cpufj = true; std::unique_ptr> root_cut_cpufj_task; - std::atomic root_cut_cpufj_preemption{false}; auto root_cut_cpufj_improvement_callback = [this](f_t obj, const std::vector& assignment, double) { std::vector user_assignment(assignment.begin(), @@ -2522,12 +2237,17 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut }; auto stop_root_cut_cpufj = [&]() { if (!root_cut_cpufj_task) { return; } - CUOPT_LOG_DEBUG("Stopping CPUFJ for this cut pass"); detail::stop_fj_cpu_task(*root_cut_cpufj_task); root_cut_cpufj_task.reset(); }; cuopt::scope_guard root_cut_cpufj_guard([&]() { stop_root_cut_cpufj(); }); + enum class cut_pass_action_t { CONTINUE, BREAK, RETURN }; + struct cut_pass_result_t { + cut_pass_action_t action{cut_pass_action_t::CONTINUE}; + mip_status_t status{mip_status_t::UNSET}; + }; + f_t cut_generation_start_time = tic(); i_t cut_pool_size = 0; for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { @@ -2536,83 +2256,298 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::OPTIMAL; } - root_cut_pass_result_t cut_pass_result; - if (root_cut_cpufj_task) { - f_t root_cut_cpufj_stop_start_time = 0.0; -#pragma omp parallel num_threads(settings_.num_threads) - { + auto do_cut_pass = [&]() -> cut_pass_result_t { +#ifdef PRINT_FRACTIONAL_INFO + settings_.log.printf( + "Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); + for (i_t j : fractional) { + settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", + j, + original_lp_.lower[j], + root_relax_soln_.x[j], + original_lp_.upper[j]); + } +#endif + + f_t cut_start_time = tic(); + bool problem_feasible = cut_generation.generate_cuts(original_lp_, + settings_, + Arow_, + new_slacks_, + var_types_, + basis_update, + root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, + basic_list, + nonbasic_list, + variable_bounds, + exploration_stats_.start_time); + if (!problem_feasible) { + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); + } + finish_clique_thread(); + return {cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; + } + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); + } + + f_t score_start_time = tic(); + cut_pool.score_cuts(root_relax_soln_.x); + f_t score_time = toc(score_start_time); + if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } + + csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); + std::vector cut_rhs; + std::vector cut_types; + i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); + if (num_cuts == 0) { return {cut_pass_action_t::BREAK, mip_status_t::UNSET}; } + cut_info.record_cut_types(cut_types); +#ifdef PRINT_CUT_POOL_TYPES + cut_pool.print_cutpool_types(); + print_cut_types("In LP ", cut_types, settings_); + printf("Cut pool size: %d\n", cut_pool.pool_size()); +#endif + +#ifdef CHECK_CUT_MATRIX + if (cuts_to_add.check_matrix() != 0) { + settings_.log.printf("Bad cuts matrix\n"); + for (i_t i = 0; i < static_cast(cut_types.size()); ++i) { + settings_.log.printf("row %d cut type %d\n", i, cut_types[i]); + } + return {cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; + } +#endif + // Check against saved solution +#ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION + verify_cuts_against_saved_solution(cuts_to_add, cut_rhs, saved_solution); +#endif + cut_pool_size = cut_pool.pool_size(); + + // Resolve the LP with the new cuts + settings_.log.debug( + "Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", + num_cuts, + cuts_to_add.row_start[cuts_to_add.m], + cut_pool.pool_size(), + cuts_to_add.m + original_lp_.num_rows); + lp_settings.log.log = false; + + f_t add_cuts_start_time = tic(); + mutex_original_lp_.lock(); + i_t add_cuts_status = add_cuts(settings_, + cuts_to_add, + cut_rhs, + original_lp_, + new_slacks_, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); + variable_bounds.resize(original_lp_.num_cols); + mutex_original_lp_.unlock(); + f_t add_cuts_time = toc(add_cuts_start_time); + if (add_cuts_time > 1.0) { + settings_.log.debug("Add cuts time %.2f seconds\n", add_cuts_time); + } + if (add_cuts_status != 0) { + settings_.log.printf("Failed to add cuts\n"); + return {cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; + } + + if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) { + mutex_upper_.lock(); + last_upper_bound = upper_bound_.load(); + std::vector lower_bounds; + std::vector upper_bounds; + find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); + mutex_upper_.unlock(); + mutex_original_lp_.lock(); + original_lp_.lower = lower_bounds; + original_lp_.upper = upper_bounds; + mutex_original_lp_.unlock(); + } + + // Try to do bound strengthening + std::vector bounds_changed(original_lp_.num_cols, true); + std::vector row_sense; +#ifdef CHECK_MATRICES + settings_.log.printf("Before A check\n"); + original_lp_.A.check_matrix(); +#endif + original_lp_.A.to_compressed_row(Arow_); + + f_t node_presolve_start_time = tic(); + bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); + std::vector new_lower = original_lp_.lower; + std::vector new_upper = original_lp_.upper; + bool feasible = + node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); + mutex_original_lp_.lock(); + original_lp_.lower = new_lower; + original_lp_.upper = new_upper; + mutex_original_lp_.unlock(); + f_t node_presolve_time = toc(node_presolve_start_time); + if (node_presolve_time > 1.0) { + settings_.log.debug("Node presolve time %.2f seconds\n", node_presolve_time); + } + if (!feasible) { + settings_.log.printf("Bound strengthening detected infeasibility\n"); +#ifdef WRITE_BOUND_STRENGTHENING_INFEASIBLE_MPS + original_lp_.write_mps("bound_strengthening_infeasible.mps"); +#endif + return {cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; + } + + i_t iter = 0; + bool initialize_basis = false; + lp_settings.concurrent_halt = NULL; + f_t dual_phase2_start_time = tic(); + dual::status_t cut_status = dual_phase2_with_advanced_basis(2, + 0, + initialize_basis, + exploration_stats_.start_time, + original_lp_, + lp_settings, + root_vstatus_, + basis_update, + basic_list, + nonbasic_list, + root_relax_soln_, + iter, + edge_norms_); + exploration_stats_.total_lp_iters += iter; + f_t dual_phase2_time = toc(dual_phase2_start_time); + if (dual_phase2_time > 1.0) { + settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); + } + if (cut_status == dual::status_t::TIME_LIMIT) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, root_objective_); + return {cut_pass_action_t::RETURN, solver_status_}; + } + + if (cut_status != dual::status_t::OPTIMAL) { + settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); + lp_status_t scratch_status = + solve_linear_program_with_advanced_basis(original_lp_, + exploration_stats_.start_time, + lp_settings, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + if (scratch_status == lp_status_t::OPTIMAL) { + // We recovered + cut_status = convert_lp_status_to_dual_status(scratch_status); + exploration_stats_.total_lp_iters += root_relax_soln_.iterations; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + } else { + settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); +#ifdef WRITE_CUT_INFEASIBLE_MPS + original_lp_.write_mps("cut_infeasible.mps"); +#endif + return {cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; + } + } + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + + f_t remove_cuts_start_time = tic(); + mutex_original_lp_.lock(); + remove_cuts(original_lp_, + settings_, + exploration_stats_.start_time, + Arow_, + new_slacks_, + original_rows, + var_types_, + root_vstatus_, + edge_norms_, + root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, + basic_list, + nonbasic_list, + basis_update); + variable_bounds.resize(original_lp_.num_cols); + mutex_original_lp_.unlock(); + f_t remove_cuts_time = toc(remove_cuts_start_time); + if (remove_cuts_time > 1.0) { + settings_.log.debug("Remove cuts time %.2f seconds\n", remove_cuts_time); + } + fractional.clear(); + num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + + if (num_fractional == 0) { + upper_bound_ = root_objective_; + mutex_upper_.lock(); + incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); + mutex_upper_.unlock(); + } + f_t obj = upper_bound_.load(); + report(' ', obj, root_objective_, 0, num_fractional); + + f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); + f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound_.load(), root_objective_); + if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { + if (num_fractional == 0) { set_solution_at_root(solution, cut_info); } + set_final_solution(solution, root_objective_); + return {cut_pass_action_t::RETURN, mip_status_t::OPTIMAL}; + } + + f_t change_in_objective = root_objective_ - last_objective; + const f_t factor = settings_.cut_change_threshold; + const f_t min_objective = 1e-3; + if (factor > 0.0 && + change_in_objective <= factor * std::max(min_objective, std::abs(root_relax_objective))) { + settings_.log.printf( + "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", + change_in_objective, + root_relax_objective); + return {cut_pass_action_t::BREAK, mip_status_t::UNSET}; + } + last_objective = root_objective_; + return {cut_pass_action_t::CONTINUE, mip_status_t::UNSET}; + }; + + cut_pass_result_t cut_pass_result; +#pragma omp parallel num_threads(root_cut_cpufj_task ? 2 : 1) + { #pragma omp single - { -#pragma omp taskgroup - { -#pragma omp task shared(root_cut_cpufj_task) firstprivate(root_cut_cpufj_work_units) default(none) - detail::run_fj_cpu_task(*root_cut_cpufj_task, - std::numeric_limits::infinity(), - std::numeric_limits::infinity()); - - cut_pass_result = run_root_cut_pass(cut_pass, - solution, - cut_generation, - cut_pool, - variable_bounds, - basis_update, - lp_settings, - basic_list, - nonbasic_list, - edge_norms_, - fractional, - num_fractional, - cut_info, - saved_solution, - last_upper_bound, - last_objective, - root_relax_objective, - original_rows, - cut_pool_size); - - root_cut_cpufj_stop_start_time = tic(); - detail::stop_fj_cpu_task(*root_cut_cpufj_task); - } - settings_.log.debug("Root cut CPUFJ stop wait time: %.6f seconds\n", - toc(root_cut_cpufj_stop_start_time)); + { + if (root_cut_cpufj_task) { +#pragma omp task shared(root_cut_cpufj_task) default(none) + detail::run_fj_cpu_task(*root_cut_cpufj_task, + std::numeric_limits::infinity(), + std::numeric_limits::infinity()); } + + cut_pass_result = do_cut_pass(); + + if (root_cut_cpufj_task) { detail::stop_fj_cpu_task(*root_cut_cpufj_task); } +#pragma omp taskwait } - root_cut_cpufj_task.reset(); - } else { - cut_pass_result = run_root_cut_pass(cut_pass, - solution, - cut_generation, - cut_pool, - variable_bounds, - basis_update, - lp_settings, - basic_list, - nonbasic_list, - edge_norms_, - fractional, - num_fractional, - cut_info, - saved_solution, - last_upper_bound, - last_objective, - root_relax_objective, - original_rows, - cut_pool_size); - } - - if (cut_pass_result.action == root_cut_pass_action_t::RETURN) { return cut_pass_result.status; } - if (cut_pass_result.action == root_cut_pass_action_t::BREAK) { break; } + } + + if (cut_pass_result.action == cut_pass_action_t::RETURN) { return cut_pass_result.status; } + if (cut_pass_result.action == cut_pass_action_t::BREAK) { break; } if (enable_root_cut_cpufj && !settings_.deterministic && settings_.num_threads >= 2 && - cut_pass + 1 < settings_.max_cut_passes && original_lp_.num_rows > original_rows) { - root_cut_cpufj_preemption = false; + cut_pass + 1 < settings_.max_cut_passes) { f_t root_cut_cpufj_build_start_time = tic(); root_cut_cpufj_task = detail::make_fj_cpu_task_from_host_lp(original_lp_, var_types_, root_relax_soln_.x, settings_, - root_cut_cpufj_preemption, root_cut_cpufj_improvement_callback, "[RootCut CPUFJ] "); settings_.log.debug("Root cut CPUFJ problem build time after pass %d: %.6f seconds\n", @@ -2632,15 +2567,13 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut original_lp_.A.col_start[original_lp_.A.n]); } - if (enable_root_cut_cpufj && original_lp_.num_rows > original_rows) { - root_cut_cpufj_preemption = false; + if (enable_root_cut_cpufj && cut_info.has_cuts()) { f_t root_cut_cpufj_build_start_time = tic(); root_cut_cpufj_task = detail::make_fj_cpu_task_from_host_lp(original_lp_, var_types_, root_relax_soln_.x, settings_, - root_cut_cpufj_preemption, root_cut_cpufj_improvement_callback, "[RootCut CPUFJ] "); settings_.log.debug("Root cut CPUFJ final problem build time: %.6f seconds\n", diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index e69798bb58..f2917ba930 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -158,33 +158,6 @@ class branch_and_bound_t { producer_sync_t& get_producer_sync() { return producer_sync_; } private: - enum class root_cut_pass_action_t { CONTINUE, BREAK, RETURN }; - - struct root_cut_pass_result_t { - root_cut_pass_action_t action{root_cut_pass_action_t::CONTINUE}; - mip_status_t status{mip_status_t::UNSET}; - }; - - root_cut_pass_result_t run_root_cut_pass(i_t cut_pass, - mip_solution_t& solution, - cut_generation_t& cut_generation, - cut_pool_t& cut_pool, - variable_bounds_t& variable_bounds, - basis_update_mpf_t& basis_update, - simplex_solver_settings_t& lp_settings, - std::vector& basic_list, - std::vector& nonbasic_list, - std::vector& edge_norms, - std::vector& fractional, - i_t& num_fractional, - cut_info_t& cut_info, - const std::vector& saved_solution, - f_t& last_upper_bound, - f_t& last_objective, - f_t root_relax_objective, - i_t original_rows, - i_t& cut_pool_size); - const user_problem_t& original_problem_; const simplex_solver_settings_t settings_; const probing_implied_bound_t& probing_implied_bound_; diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh b/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh index a041957ece..9a828ea43a 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh @@ -49,6 +49,7 @@ struct fj_cpu_task_t { struct fj_cpu_deleter_t { void operator()(fj_cpu_climber_t* ptr) const; }; + std::atomic preemption_flag{false}; std::unique_ptr, fj_cpu_deleter_t> fj_cpu; }; @@ -58,7 +59,6 @@ std::unique_ptr> make_fj_cpu_task_from_host_lp( const std::vector& variable_types, const std::vector& seed_assignment, const dual_simplex::simplex_solver_settings_t& settings, - std::atomic& preemption_flag, std::function&, double)> improvement_callback, std::string log_prefix); diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 09ded977f9..797e0354e3 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -1770,7 +1770,7 @@ static bool cpufj_solve_loop(fj_cpu_climber_t& fj_cpu, fj_cpu.work_units_elapsed += biased_work; if (fj_cpu.producer_sync != nullptr) { fj_cpu.producer_sync->notify_progress(); } - if (fj_cpu.work_units_elapsed.load(std::memory_order_acquire) >= work_unit_limit) { break; } + if (fj_cpu.work_units_elapsed >= work_unit_limit) { break; } } cuopt_func_call(sanity_checks(fj_cpu)); @@ -1863,13 +1863,12 @@ std::unique_ptr> make_fj_cpu_task_from_host_lp( const std::vector& variable_types, const std::vector& seed_assignment, const dual_simplex::simplex_solver_settings_t& settings, - std::atomic& preemption_flag, std::function&, double)> improvement_callback, std::string log_prefix) { - auto task = std::make_unique>(); - auto fj_cpu = - init_fj_cpu_from_host_lp(problem, variable_types, seed_assignment, settings, preemption_flag); + auto task = std::make_unique>(); + auto fj_cpu = init_fj_cpu_from_host_lp( + problem, variable_types, seed_assignment, settings, task->preemption_flag); fj_cpu->log_prefix = std::move(log_prefix); fj_cpu->improvement_callback = std::move(improvement_callback); task->fj_cpu.reset(fj_cpu.release()); @@ -1909,7 +1908,6 @@ template std::unique_ptr> make_fj_cpu_task_from_host_l const std::vector& variable_types, const std::vector& seed_assignment, const dual_simplex::simplex_solver_settings_t& settings, - std::atomic& preemption_flag, std::function&, double)> improvement_callback, std::string log_prefix); template void run_fj_cpu_task(fj_cpu_task_t& task, @@ -1939,7 +1937,6 @@ template std::unique_ptr> make_fj_cpu_task_from_host_ const std::vector& variable_types, const std::vector& seed_assignment, const dual_simplex::simplex_solver_settings_t& settings, - std::atomic& preemption_flag, std::function&, double)> improvement_callback, std::string log_prefix); template void run_fj_cpu_task(fj_cpu_task_t& task, diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh index 8982dfa5f7..cfac5121dd 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh @@ -11,7 +11,6 @@ #include #include #include -#include #include #include #include diff --git a/cpp/src/mip_heuristics/local_search/local_search.cu b/cpp/src/mip_heuristics/local_search/local_search.cu index b96b48a413..135f889921 100644 --- a/cpp/src/mip_heuristics/local_search/local_search.cu +++ b/cpp/src/mip_heuristics/local_search/local_search.cu @@ -125,7 +125,7 @@ void local_search_t::start_cpufj_lptopt_scratch_threads( solution_lp, default_weights, default_weights, 0., context.preempt_heuristic_solver_); scratch_cpu_fj_on_lp_opt.fj_cpu->log_prefix = "******* scratch on LP optimal: "; scratch_cpu_fj_on_lp_opt.fj_cpu->improvement_callback = - [&population](f_t obj, const std::vector& h_vec, double /*work_units*/) { + [&population, this](f_t obj, const std::vector& h_vec, double /*work_units*/) { population.add_external_solution(h_vec, obj, solution_origin_t::CPUFJ); if (obj < local_search_best_obj) { CUOPT_LOG_DEBUG("******* New local search best obj %g, best overall %g", From 45c81542992dc1fa4b01c457199fcbf0280f6e9e Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 5 May 2026 01:33:29 -0700 Subject: [PATCH 5/6] some cleanup --- cpp/src/branch_and_bound/branch_and_bound.cpp | 6 +++--- cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu | 7 ++++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index f56d0d3e79..b4f18b0700 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2232,7 +2232,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut [this](f_t obj, const std::vector& assignment, double) { std::vector user_assignment(assignment.begin(), assignment.begin() + original_problem_.num_cols); - CUOPT_LOG_DEBUG("Root cut CPUFJ found solution with objective %.16e\n", obj); + settings_.log.debug("Root cut CPUFJ found solution with objective %.16e\n", obj); set_new_solution(user_assignment); }; auto stop_root_cut_cpufj = [&]() { @@ -2294,12 +2294,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut if (cut_generation_time > 1.0) { settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); } - + // Score the cuts f_t score_start_time = tic(); cut_pool.score_cuts(root_relax_soln_.x); f_t score_time = toc(score_start_time); if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } - + // Get the best cuts from the cut pool csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); std::vector cut_rhs; std::vector cut_types; diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 797e0354e3..4ef68abfb5 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -1419,6 +1419,7 @@ void finalize_fj_cpu_host_initialization( // nnz count fj_cpu.cached_mtm_moves.resize(fj_cpu.h_coefficients.size(), std::make_pair(0, fj_staged_score_t::zero())); + fj_cpu.cached_cstr_bounds.resize(fj_cpu.h_reverse_coefficients.size()); for (i_t var_idx = 0; var_idx < n_variables; ++var_idx) { auto [offset_begin, offset_end] = reverse_range_for_var(fj_cpu, var_idx); @@ -1888,9 +1889,9 @@ template void stop_fj_cpu_task(fj_cpu_task_t& task) { if (task.fj_cpu) { - auto& fj_cpu = *task.fj_cpu; - fj_cpu.preemption_flag.store(true, std::memory_order_release); - fj_cpu.halted = true; + auto& fj_cpu = *task.fj_cpu; + fj_cpu.preemption_flag = true; + fj_cpu.halted = true; } } From 64ba333e8142b955560f701f5c7050034ede7695 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 5 May 2026 01:51:08 -0700 Subject: [PATCH 6/6] AI review comments --- cpp/src/branch_and_bound/branch_and_bound.cpp | 4 ++-- cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index b4f18b0700..35225b042a 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2230,8 +2230,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut std::unique_ptr> root_cut_cpufj_task; auto root_cut_cpufj_improvement_callback = [this](f_t obj, const std::vector& assignment, double) { - std::vector user_assignment(assignment.begin(), - assignment.begin() + original_problem_.num_cols); + std::vector user_assignment; + uncrush_primal_solution(original_problem_, original_lp_, assignment, user_assignment); settings_.log.debug("Root cut CPUFJ found solution with objective %.16e\n", obj); set_new_solution(user_assignment); }; diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 4ef68abfb5..12758b509f 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -1455,6 +1455,7 @@ static std::unique_ptr> init_fj_cpu_from_host_lp( typename mip_solver_settings_t::tolerances_t tolerances{}; tolerances.absolute_tolerance = settings.primal_tol; + tolerances.relative_tolerance = settings.zero_tol; tolerances.integrality_tolerance = settings.integer_tol; tolerances.absolute_mip_gap = settings.absolute_mip_gap_tol; tolerances.relative_mip_gap = settings.relative_mip_gap_tol;