From a123da29b0af0ebb69e91b048fd94bf6ec7a05a0 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 19 Sep 2025 22:23:58 +0200 Subject: [PATCH 01/33] parallel mip solver with support for best first search with plunging and diving threads. --- cpp/src/dual_simplex/branch_and_bound.cpp | 603 +++++++++++------- cpp/src/dual_simplex/branch_and_bound.hpp | 112 ++-- cpp/src/dual_simplex/mip_node.hpp | 27 +- cpp/src/dual_simplex/pseudo_costs.cpp | 141 ++-- .../dual_simplex/simplex_solver_settings.hpp | 9 +- cpp/src/math_optimization/solver_settings.cu | 2 +- cpp/src/mip/diversity/recombiners/sub_mip.cuh | 1 + 7 files changed, 531 insertions(+), 364 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 3d41921054..5d18f253cb 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -17,7 +17,6 @@ #include #include - #include #include #include @@ -29,9 +28,10 @@ #include #include +#include #include #include -#include +#include #include #include #include @@ -229,7 +229,8 @@ branch_and_bound_t::branch_and_bound_t( original_lp_(1, 1, 1), incumbent_(1), root_relax_soln_(1, 1), - pc_(1) + pc_(1), + status_(mip_status_t::UNSET) { stats_.start_time = tic(); convert_user_problem(original_problem_, settings_, original_lp_, new_slacks_); @@ -238,14 +239,6 @@ branch_and_bound_t::branch_and_bound_t( mutex_upper_.lock(); upper_bound_ = inf; mutex_upper_.unlock(); - - mutex_lower_.lock(); - lower_bound_ = -inf; - mutex_lower_.unlock(); - - mutex_branching_.lock(); - currently_branching_ = false; - mutex_branching_.unlock(); } template @@ -257,15 +250,44 @@ f_t branch_and_bound_t::get_upper_bound() return upper_bound; } +// Calculates the global lower bound considering +// the lower bounds in each thread and the top of the heap. template f_t branch_and_bound_t::get_lower_bound() { - mutex_lower_.lock(); - const f_t lower_bound = lower_bound_; - mutex_lower_.unlock(); + f_t lower_bound = -inf; + + mutex_heap_.lock(); + if (heap_.size() > 0) { + lower_bound = heap_.top()->lower_bound; + } else if (search_tree_ != nullptr) { + lower_bound = search_tree_->lower_bound; + } + mutex_heap_.unlock(); + return lower_bound; } +template +mip_status_t branch_and_bound_t::get_status() +{ + mip_status_t status; + +#pragma omp atomic read + status = status_; + + return status; +} + +template +i_t branch_and_bound_t::get_heap_size() +{ + mutex_heap_.lock(); + i_t size = heap_.size(); + mutex_heap_.unlock(); + return size; +} + template void branch_and_bound_t::set_new_solution(const std::vector& solution) { @@ -304,16 +326,13 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu mutex_upper_.unlock(); if (is_feasible) { - mutex_branching_.lock(); - bool currently_branching = currently_branching_; - mutex_branching_.unlock(); - if (currently_branching) { + if (get_status() == mip_status_t::RUNNING) { f_t user_obj = compute_user_objective(original_lp_, obj); f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); std::string gap = user_mip_gap(user_obj, user_lower); settings_.log.printf( - "H %+13.6e %+10.6e %s %9.2f\n", + "H %+13.6e %+10.6e %s %9.2f\n", user_obj, user_lower, gap.c_str(), @@ -424,6 +443,7 @@ void branch_and_bound_t::repair_heuristic_solutions() f_t obj = compute_user_objective(original_lp_, repaired_obj); f_t lower = compute_user_objective(original_lp_, get_lower_bound()); std::string user_gap = user_mip_gap(obj, lower); + settings_.log.printf( "H %+13.6e %+10.6e %s %9.2f\n", obj, @@ -450,16 +470,24 @@ void branch_and_bound_t::branch(mip_node_t* parent_node, f_t branch_var_val, const std::vector& parent_vstatus) { + i_t id; + +#pragma omp atomic capture + { + id = stats_.num_nodes; + stats_.num_nodes += 2; + } + // down child auto down_child = std::make_unique>( - original_lp_, parent_node, ++stats_.num_nodes, branch_var, 0, branch_var_val, parent_vstatus); + original_lp_, parent_node, ++id, branch_var, 0, branch_var_val, parent_vstatus); graphviz_edge( settings_, parent_node, down_child.get(), branch_var, 0, std::floor(branch_var_val)); // up child auto up_child = std::make_unique>( - original_lp_, parent_node, ++stats_.num_nodes, branch_var, 1, branch_var_val, parent_vstatus); + original_lp_, parent_node, ++id, branch_var, 1, branch_var_val, parent_vstatus); graphviz_edge(settings_, parent_node, up_child.get(), branch_var, 1, std::ceil(branch_var_val)); @@ -472,8 +500,10 @@ template void branch_and_bound_t::update_tree(mip_node_t* node_ptr, node_status_t status) { std::vector*> stack; + mutex_search_tree_.lock(); node_ptr->set_status(status, stack); remove_fathomed_nodes(stack); + mutex_search_tree_.unlock(); } template @@ -482,20 +512,24 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, i_t leaf_depth, char symbol) { - bool send_solution = false; - i_t nodes_explored = stats_.nodes_explored; - i_t nodes_unexplored = stats_.nodes_unexplored; - f_t gap; + bool send_solution = false; + i_t nodes_explored; + i_t nodes_unexplored; + +#pragma omp atomic read + nodes_explored = stats_.nodes_explored; + +#pragma omp atomic read + nodes_unexplored = stats_.nodes_unexplored; mutex_upper_.lock(); if (leaf_objective < upper_bound_) { incumbent_.set_incumbent_solution(leaf_objective, leaf_solution); upper_bound_ = leaf_objective; f_t lower_bound = get_lower_bound(); - gap = upper_bound_ - lower_bound; f_t obj = compute_user_objective(original_lp_, upper_bound_); f_t lower = compute_user_objective(original_lp_, lower_bound); - settings_.log.printf("%c%8d %8lu %+13.6e %+10.6e %4d %7.1e %s %9.2f\n", + settings_.log.printf("%c%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", symbol, nodes_explored, nodes_unexplored, @@ -515,12 +549,6 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, settings_.solution_callback(original_x, upper_bound_); } mutex_upper_.unlock(); - - if (send_solution) { - mutex_gap_.lock(); - gap_ = gap; - mutex_gap_.unlock(); - } } template @@ -567,36 +595,36 @@ dual::status_t branch_and_bound_t::node_dual_simplex( leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms); lp_status = convert_lp_status_to_dual_status(second_status); } - } else { - log.printf("Infeasible after bounds strengthening. Fathoming node %d.\n", leaf_id); } - mutex_stats_.lock(); +#pragma omp atomic update stats_.total_lp_solve_time += toc(lp_start_time); + +#pragma omp atomic update stats_.total_lp_iters += node_iter; - mutex_stats_.unlock(); return lp_status; } template -mip_status_t branch_and_bound_t::solve_node_lp( - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - csc_matrix_t& Arow, - const std::vector& var_types, - f_t upper_bound) +mip_status_t branch_and_bound_t::solve_node_lp(mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow, + f_t upper_bound) { logger_t log; - log.log = false; + log.log = false; + + f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; + f_t rel_fathom_tol = settings_.relative_mip_gap_tol; + std::vector& leaf_vstatus = node_ptr->vstatus; lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); // Set the correct bounds for the leaf problem leaf_problem.lower = original_lp_.lower; leaf_problem.upper = original_lp_.upper; - - std::vector bounds_changed(original_lp_.num_cols, false); + std::vector bounds_changed(leaf_problem.num_cols, false); // Technically, we can get the already strengthened bounds from the node/parent instead of // getting it from the original problem and re-strengthening. But this requires storing // two vectors at each node and potentially cause memory issues @@ -615,6 +643,7 @@ mip_status_t branch_and_bound_t::solve_node_lp( node_ptr->lower_bound = inf; graphviz_node(settings_, node_ptr, "infeasible", 0.0); update_tree(node_ptr, node_status_t::INFEASIBLE); + // Node was infeasible. Do not branch } else if (lp_status == dual::status_t::CUTOFF) { @@ -622,12 +651,13 @@ mip_status_t branch_and_bound_t::solve_node_lp( f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); graphviz_node(settings_, node_ptr, "cut off", leaf_objective); update_tree(node_ptr, node_status_t::FATHOMED); + // Node was cut off. Do not branch } else if (lp_status == dual::status_t::OPTIMAL) { // LP was feasible - std::vector fractional; - const i_t leaf_fractional = - fractional_variables(settings_, leaf_solution.x, var_types_, fractional); + std::vector leaf_fractional; + i_t leaf_num_fractional = + fractional_variables(settings_, leaf_solution.x, var_types_, leaf_fractional); f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); graphviz_node(settings_, node_ptr, "lower bound", leaf_objective); @@ -637,17 +667,17 @@ mip_status_t branch_and_bound_t::solve_node_lp( node_ptr->lower_bound = leaf_objective; - constexpr f_t fathom_tol = 1e-5; - if (leaf_fractional == 0) { + if (leaf_num_fractional == 0) { add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, 'B'); graphviz_node(settings_, node_ptr, "integer feasible", leaf_objective); update_tree(node_ptr, node_status_t::INTEGER_FEASIBLE); - } else if (leaf_objective <= upper_bound + fathom_tol) { + } else if (leaf_objective <= upper_bound + abs_fathom_tol && + relative_gap(upper_bound, leaf_objective) > rel_fathom_tol) { // Choose fractional variable to branch on mutex_pc_.lock(); const i_t branch_var = pc_.variable_selection( - fractional, leaf_solution.x, leaf_problem.lower, leaf_problem.upper, log); + leaf_fractional, leaf_solution.x, leaf_problem.lower, leaf_problem.upper, log); mutex_pc_.unlock(); assert(leaf_vstatus.size() == leaf_problem.num_cols); @@ -658,210 +688,330 @@ mip_status_t branch_and_bound_t::solve_node_lp( graphviz_node(settings_, node_ptr, "fathomed", leaf_objective); update_tree(node_ptr, node_status_t::FATHOMED); } + } else if (lp_status == dual::status_t::TIME_LIMIT) { + graphviz_node(settings_, node_ptr, "timeout", 0.0); + return mip_status_t::TIME_LIMIT; + } else { graphviz_node(settings_, node_ptr, "numerical", 0.0); - settings_.log.printf("Encountered LP status %d. This indicates a numerical issue.\n", - lp_status); + settings_.log.printf("Encountered LP status %d on node %d. This indicates a numerical issue.\n", + lp_status, + node_ptr->node_id); return mip_status_t::NUMERICAL; } - return mip_status_t::UNSET; + return mip_status_t::RUNNING; } template -mip_status_t branch_and_bound_t::explore_tree(i_t branch_var, - mip_solution_t& solution) +void branch_and_bound_t::explore_subtree(mip_node_t* start_node) { - mip_status_t status = mip_status_t::UNSET; - logger_t log; - log.log = false; - auto compare = [](mip_node_t* a, mip_node_t* b) { - return a->lower_bound > - b->lower_bound; // True if a comes before b, elements that come before are output last - }; - - std::priority_queue*, std::vector*>, decltype(compare)> - heap(compare); - - mip_node_t root_node(root_objective_, root_vstatus_); - graphviz_node(settings_, &root_node, "lower bound", root_objective_); - - branch(&root_node, branch_var, root_relax_soln_.x[branch_var], root_vstatus_); - - // the stack does not own the unique_ptr the tree does - heap.push(root_node.get_down_child()); - heap.push(root_node.get_up_child()); - // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; csc_matrix_t Arow(1, 1, 1); - original_lp_.A.transpose(Arow); + leaf_problem.A.transpose(Arow); - f_t lower_bound = get_lower_bound(); - f_t gap = get_upper_bound() - lower_bound; - f_t last_log = 0; + std::deque*> stack; + stack.push_front(start_node); - while (gap > settings_.absolute_mip_gap_tol && - relative_gap(get_upper_bound(), lower_bound) > settings_.relative_mip_gap_tol && - heap.size() > 0) { + // It contains the lower bound of the parent for now + f_t lower_bound = start_node->lower_bound; + f_t upper_bound = get_upper_bound(); + f_t gap = upper_bound - lower_bound; + f_t last_log = 0; + i_t nodes_explored = 0; + i_t nodes_unexplored = 0; + + while (stack.size() > 0 && get_status() == mip_status_t::RUNNING) { repair_heuristic_solutions(); - // Get a node off the heap - mip_node_t* node_ptr = heap.top(); - heap.pop(); - stats_.nodes_explored++; + mip_node_t* node_ptr = stack.front(); + stack.pop_front(); - f_t upper_bound = get_upper_bound(); - if (upper_bound < node_ptr->lower_bound) { - // This node was put on the heap earlier but its lower bound is now greater than the current - // upper bound - update_tree(node_ptr, node_status_t::FATHOMED); + lower_bound = node_ptr->lower_bound; + upper_bound = get_upper_bound(); + gap = upper_bound - lower_bound; + +#pragma omp atomic capture + nodes_explored = stats_.nodes_explored++; + +#pragma omp atomic capture + nodes_unexplored = stats_.nodes_unexplored--; + + if (upper_bound < node_ptr->lower_bound || + relative_gap(upper_bound, lower_bound) < settings_.relative_mip_gap_tol) { graphviz_node(settings_, node_ptr, "cutoff", node_ptr->lower_bound); + update_tree(node_ptr, node_status_t::FATHOMED); continue; } - mutex_lower_.lock(); - lower_bound = lower_bound_ = node_ptr->lower_bound; - mutex_lower_.unlock(); - mutex_gap_.lock(); - gap_ = gap = upper_bound - lower_bound; - mutex_gap_.unlock(); + f_t now = toc(stats_.start_time); - i_t nodes_explored = stats_.nodes_explored; - f_t now = toc(stats_.start_time); f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); - if ((nodes_explored % 1000 == 0 || gap < 10 * settings_.absolute_mip_gap_tol || - nodes_explored < 1000) && - (time_since_log >= 1) || + if (((nodes_explored % 1000 == 0 || gap < 10 * settings_.absolute_mip_gap_tol || + nodes_explored < 1000) && + (time_since_log >= 1)) || (time_since_log > 60) || now > settings_.time_limit) { - f_t user_obj = compute_user_objective(original_lp_, upper_bound); - f_t user_lower = compute_user_objective(original_lp_, lower_bound); - std::string user_gap = user_mip_gap(user_obj, user_lower); + f_t obj = compute_user_objective(original_lp_, upper_bound); + f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); + std::string gap_user = user_mip_gap(obj, user_lower); - settings_.log.printf(" %8d %8lu %+13.6e %+10.6e %4d %7.1e %s %9.2f\n", + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", nodes_explored, - heap.size(), - user_obj, + nodes_unexplored, + obj, user_lower, node_ptr->depth, nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, - user_gap.c_str(), + gap_user.c_str(), now); last_log = tic(); } if (toc(stats_.start_time) > settings_.time_limit) { - settings_.log.printf("Hit time limit. Stopping\n"); - status = mip_status_t::TIME_LIMIT; +#pragma omp atomic write + status_ = mip_status_t::TIME_LIMIT; +#pragma omp single nowait + settings_.log.printf("Time limit reached. Stopping the solver...\n"); break; } - status = solve_node_lp(node_ptr, leaf_problem, Arow, var_types_, upper_bound); - - if (status == mip_status_t::NUMERICAL) { break; } + // We are just checking for numerical issues or timeout during the LP solve, otherwise + // lp_status is set to RUNNING. + mip_status_t lp_status = solve_node_lp(node_ptr, leaf_problem, Arow, upper_bound); + if (lp_status == mip_status_t::TIME_LIMIT) { +#pragma omp atomic write + status_ = lp_status; +#pragma omp single nowait + settings_.log.printf("Time limit reached. Stopping the solver...\n"); + break; + } if (node_ptr->status == node_status_t::HAS_CHILDREN) { - // the heap does not own the unique_ptr the tree does - heap.push(node_ptr->get_down_child()); - heap.push(node_ptr->get_up_child()); - } - } + // Martin's child selection + const i_t branch_var = node_ptr->get_down_child()->branch_var; + const f_t branch_var_val = node_ptr->get_down_child()->fractional_val; + const f_t down_val = std::floor(root_relax_soln_.x[branch_var]); + const f_t up_val = std::ceil(root_relax_soln_.x[branch_var]); + const f_t down_dist = branch_var_val - down_val; + const f_t up_dist = up_val - branch_var_val; - stats_.nodes_unexplored = heap.size(); + if (stack.size() > 0) { + mip_node_t* node = stack.back(); + stack.pop_back(); - if (stats_.nodes_unexplored == 0) { - mutex_lower_.lock(); - lower_bound = lower_bound_ = root_node.lower_bound; - mutex_lower_.unlock(); + mutex_heap_.lock(); + heap_.push(node); + mutex_heap_.unlock(); - mutex_gap_.lock(); - gap_ = gap = get_upper_bound() - lower_bound; - mutex_gap_.unlock(); - } + mutex_dive_queue_.lock(); + dive_queue_.push(node->create_detach_copy()); + mutex_dive_queue_.unlock(); + } - return status; + if (down_dist < up_dist) { + stack.push_front(node_ptr->get_up_child()); + stack.push_front(node_ptr->get_down_child()); + } else { + stack.push_front(node_ptr->get_down_child()); + stack.push_front(node_ptr->get_up_child()); + } + +#pragma omp atomic update + stats_.nodes_unexplored += 2; + } + } } template -mip_status_t branch_and_bound_t::dive(i_t branch_var, mip_solution_t& solution) +void branch_and_bound_t::dive(std::shared_ptr> start_node) { - mip_status_t status = mip_status_t::UNSET; - logger_t log; log.log = false; - std::vector*> node_stack; - - mip_node_t root_node(root_objective_, root_vstatus_); - graphviz_node(settings_, &root_node, "lower bound", root_objective_); - - branch(&root_node, branch_var, root_relax_soln_.x[branch_var], root_vstatus_); - - // the stack does not own the unique_ptr the tree does - node_stack.push_back(root_node.get_down_child()); - node_stack.push_back(root_node.get_up_child()); + f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; + f_t rel_fathom_tol = settings_.relative_mip_gap_tol; // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; - csc_matrix_t Arow(1, 1, 1); - original_lp_.A.transpose(Arow); + leaf_problem.A.transpose(Arow); - f_t lower_bound = get_lower_bound(); - f_t gap = get_upper_bound() - lower_bound; - i_t nodes_explored = 0; + std::deque*> stack; + stack.push_front(start_node.get()); - while (node_stack.size() > 0) { - // Get a node off the stack - mip_node_t* node_ptr = node_stack.back(); - node_stack.pop_back(); - nodes_explored++; + // It contains the lower bound of the parent for now + f_t upper_bound = get_upper_bound(); + i_t num_nodes = 1; - f_t upper_bound = get_upper_bound(); - lower_bound = get_lower_bound(); - gap = upper_bound - lower_bound; + while (stack.size() > 0 && get_status() == mip_status_t::RUNNING) { + if (stack.size() > 1 && dive_queue_.size() < 4 * settings_.num_diving_threads) { + std::lock_guard lock(mutex_dive_queue_); + mip_node_t* node = stack.back(); + stack.pop_back(); + dive_queue_.push(node->create_detach_copy()); + } - if (gap < settings_.absolute_mip_gap_tol && - relative_gap(get_upper_bound(), lower_bound) < settings_.relative_mip_gap_tol) { - update_tree(node_ptr, node_status_t::FATHOMED); + mip_node_t* node_ptr = stack.front(); + stack.pop_front(); + upper_bound = get_upper_bound(); + + if (upper_bound < node_ptr->lower_bound || + relative_gap(upper_bound, node_ptr->lower_bound) < settings_.relative_mip_gap_tol) { continue; } - if (toc(stats_.start_time) > settings_.time_limit) { - status = mip_status_t::TIME_LIMIT; - break; - } + std::vector& leaf_vstatus = node_ptr->vstatus; + lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); - status = solve_node_lp(node_ptr, leaf_problem, Arow, var_types_, upper_bound); + leaf_problem.lower = original_lp_.lower; + leaf_problem.upper = original_lp_.upper; + std::vector bounds_changed(leaf_problem.num_cols, false); + node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); - if (status == mip_status_t::NUMERICAL) { continue; } + dual::status_t lp_status = node_dual_simplex(node_ptr->node_id, + leaf_problem, + leaf_vstatus, + leaf_solution, + bounds_changed, + Arow, + upper_bound, + log); + + if (lp_status != dual::status_t::OPTIMAL) { continue; } + + std::vector fractional; + const i_t leaf_fractional = + fractional_variables(settings_, leaf_solution.x, var_types_, fractional); + f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); + + mutex_pc_.lock(); + pc_.update_pseudo_costs(node_ptr, leaf_objective); + mutex_pc_.unlock(); + + node_ptr->lower_bound = leaf_objective; + + if (leaf_fractional == 0) { + add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, 'D'); + } else if (leaf_objective <= upper_bound + abs_fathom_tol && + relative_gap(upper_bound, leaf_objective) > rel_fathom_tol) { + // Choose fractional variable to branch on + mutex_pc_.lock(); + const i_t branch_var = pc_.variable_selection( + fractional, leaf_solution.x, leaf_problem.lower, leaf_problem.upper, log); + mutex_pc_.unlock(); + const f_t branch_var_val = leaf_solution.x[branch_var]; + + assert(leaf_vstatus.size() == leaf_problem.num_cols); + assert(leaf_vstatus.size() == original_lp_.num_cols); + + // down child + auto down_child = std::make_unique>( + original_lp_, node_ptr, ++num_nodes, branch_var, 0, branch_var_val, leaf_vstatus); + + // up child + auto up_child = std::make_unique>( + original_lp_, node_ptr, ++num_nodes, branch_var, 1, branch_var_val, leaf_vstatus); - if (node_ptr->status == node_status_t::HAS_CHILDREN) { // Martin's child selection - const i_t branch_var = node_ptr->get_down_child()->branch_var; - const f_t branch_var_val = node_ptr->get_down_child()->fractional_val; - const f_t down_val = std::floor(root_relax_soln_.x[branch_var]); - const f_t up_val = std::ceil(root_relax_soln_.x[branch_var]); - const f_t down_dist = branch_var_val - down_val; - const f_t up_dist = up_val - branch_var_val; + const f_t down_val = std::floor(root_relax_soln_.x[branch_var]); + const f_t up_val = std::ceil(root_relax_soln_.x[branch_var]); + const f_t down_dist = -down_val; + const f_t up_dist = up_val - leaf_solution.x[branch_var]; if (down_dist < up_dist) { - node_stack.push_back(node_ptr->get_up_child()); - node_stack.push_back(node_ptr->get_down_child()); + stack.push_front(up_child.get()); + stack.push_front(down_child.get()); } else { - node_stack.push_back(node_ptr->get_down_child()); - node_stack.push_back(node_ptr->get_up_child()); + stack.push_front(down_child.get()); + stack.push_front(up_child.get()); } + + node_ptr->add_children(std::move(down_child), std::move(up_child)); } } +} - return status; +template +void branch_and_bound_t::best_first_thread() +{ + i_t active = 0; + f_t lower_bound = -inf; + f_t upper_bound = inf; + f_t gap = inf; + + do { + mip_node_t* node_ptr = nullptr; + + // If there any node left in the heap, we pop the top node and explore it. + mutex_heap_.lock(); + if (heap_.size() > 0) { + node_ptr = heap_.top(); + heap_.pop(); + } + mutex_heap_.unlock(); + + if (node_ptr != nullptr) { + if (get_upper_bound() < node_ptr->lower_bound) { + // This node was put on the heap earlier but its lower bound is now greater than the + // current upper bound + graphviz_node(settings_, node_ptr, "cutoff", node_ptr->lower_bound); + update_tree(node_ptr, node_status_t::FATHOMED); + continue; + } + +#pragma omp atomic update + active_best_first_tasks_++; + + // Best-first search with plunging + explore_subtree(node_ptr); + +#pragma omp atomic update + active_best_first_tasks_--; + } + +#pragma omp atomic read + active = active_best_first_tasks_; + + lower_bound = get_lower_bound(); + upper_bound = get_upper_bound(); + gap = upper_bound - lower_bound; + + } while (get_status() == mip_status_t::RUNNING && gap > settings_.absolute_mip_gap_tol && + relative_gap(upper_bound, lower_bound) > settings_.relative_mip_gap_tol && + (active > 0 || get_heap_size() > 0)); +} + +template +void branch_and_bound_t::diving_thread() +{ + while (get_status() == mip_status_t::RUNNING) { + std::shared_ptr> node_ptr; + + mutex_dive_queue_.lock(); + if (dive_queue_.size() > 0) { node_ptr = dive_queue_.pop(); } + mutex_dive_queue_.unlock(); + + if (node_ptr) { + if (get_upper_bound() < node_ptr->lower_bound) { continue; } + dive(node_ptr); + } + } } template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { - mip_status_t status = mip_status_t::UNSET; + logger_t log; + + log.log = false; + status_ = mip_status_t::UNSET; + + stats_.total_lp_iters = 0; + stats_.nodes_explored = 0; + stats_.nodes_unexplored = 2; + stats_.num_nodes = 0; + active_best_first_tasks_ = 0; if (guess_.size() != 0) { std::vector crushed_guess; @@ -910,6 +1060,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.log.printf("Hit time limit\n"); return mip_status_t::TIME_LIMIT; } + set_uninitialized_steepest_edge_norms(edge_norms_); root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); @@ -927,9 +1078,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.set_simplex_solution_callback( original_x, original_dual, compute_user_objective(original_lp_, root_objective_)); } - mutex_lower_.lock(); - lower_bound_ = root_objective_; - mutex_lower_.unlock(); std::vector fractional; const i_t num_fractional = @@ -943,7 +1091,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut // We should be done here uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); solution.objective = incumbent_.objective; - solution.lower_bound = lower_bound_; + solution.lower_bound = root_objective_; solution.nodes_explored = 0; solution.simplex_iterations = root_relax_soln_.iterations; settings_.log.printf("Optimal solution found at root node. Objective %.16e. Time %.2f.\n", @@ -958,7 +1106,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::OPTIMAL; } - pc_.resize(original_lp_.num_cols); strong_branching(original_lp_, settings_, stats_.start_time, @@ -970,55 +1117,67 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut edge_norms_, pc_); + if (toc(stats_.start_time) > settings_.time_limit) { + settings_.log.printf("Hit time limit\n"); + return mip_status_t::TIME_LIMIT; + } + // Choose variable to branch on - logger_t log; - log.log = false; i_t branch_var = pc_.variable_selection( fractional, root_relax_soln_.x, original_lp_.lower, original_lp_.upper, log); - stats_.total_lp_iters = 0; - stats_.nodes_explored = 0; - stats_.nodes_unexplored = 0; - stats_.num_nodes = 1; + search_tree_ = std::make_unique>(root_objective_, root_vstatus_); + graphviz_node(settings_, search_tree_.get(), "lower bound", root_objective_); + branch(search_tree_.get(), branch_var, root_relax_soln_.x[branch_var], root_vstatus_); + heap_.push(search_tree_->get_down_child()); + heap_.push(search_tree_->get_up_child()); + + settings_.log.printf("Exploring the B&B tree using %d threads and %d diving threads\n", + settings_.num_bfs_threads, + settings_.num_threads - settings_.num_bfs_threads); settings_.log.printf( - "| Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap | " + "| Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap | " " Time \n"); - mutex_branching_.lock(); - currently_branching_ = true; - mutex_branching_.unlock(); +#pragma omp atomic write + status_ = mip_status_t::RUNNING; - std::future diving_thread; +#pragma omp parallel num_threads(settings_.num_threads) + { + if (omp_get_thread_num() < settings_.num_bfs_threads) { + best_first_thread(); + } else { + diving_thread(); + } - if (settings_.num_threads > 0) { - diving_thread = std::async(std::launch::async, [&]() { return dive(branch_var, solution); }); +// Check if the solver exited normally. +#pragma omp single nowait + { + if (status_ == mip_status_t::RUNNING) { +#pragma omp atomic write + status_ = mip_status_t::FINISHED; + } + } } - status = explore_tree(branch_var, solution); - - if (settings_.num_threads > 0) { mip_status_t diving_status = diving_thread.get(); } - - mutex_branching_.lock(); - currently_branching_ = false; - mutex_branching_.unlock(); + f_t lower_bound = get_lower_bound(); + f_t upper_bound = get_upper_bound(); + f_t gap = upper_bound - lower_bound; + f_t obj = compute_user_objective(original_lp_, upper_bound); + f_t user_lower = compute_user_objective(original_lp_, lower_bound); + f_t gap_rel = relative_gap(get_upper_bound(), lower_bound); settings_.log.printf( - "Explored %d nodes in %.2fs.\nAbsolute Gap %e Objective %.16e Lower Bound %.16e\n", - stats_.nodes_explored.load(), - toc(stats_.start_time), - gap_, - compute_user_objective(original_lp_, get_upper_bound()), - compute_user_objective(original_lp_, lower_bound_)); - - if (gap_ <= settings_.absolute_mip_gap_tol || - relative_gap(get_upper_bound(), lower_bound_) <= settings_.relative_mip_gap_tol) { - status = mip_status_t::OPTIMAL; - if (gap_ > 0 && gap_ <= settings_.absolute_mip_gap_tol) { + "Explored %d nodes in %.2fs.\n", stats_.nodes_explored, toc(stats_.start_time)); + settings_.log.printf("Absolute Gap %e Objective %.16e Lower Bound %.16e\n", gap, obj, user_lower); + + if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) { + status_ = mip_status_t::OPTIMAL; + if (gap > 0 && gap <= settings_.absolute_mip_gap_tol) { settings_.log.printf("Optimal solution found within absolute MIP gap tolerance (%.1e)\n", settings_.absolute_mip_gap_tol); - } else if (gap_ > 0 && - relative_gap(get_upper_bound(), lower_bound_) <= settings_.relative_mip_gap_tol) { + } else if (gap > 0 && gap_rel <= settings_.relative_mip_gap_tol) { settings_.log.printf("Optimal solution found within relative MIP gap tolerance (%.1e)\n", settings_.relative_mip_gap_tol); } else { @@ -1029,9 +1188,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } } - if (stats_.nodes_unexplored == 0 && get_upper_bound() == inf) { + if (get_heap_size() == 0 && upper_bound == inf) { settings_.log.printf("Integer infeasible.\n"); - status = mip_status_t::INFEASIBLE; + status_ = mip_status_t::INFEASIBLE; if (settings_.heuristic_preemption_callback != nullptr) { settings_.heuristic_preemption_callback(); } @@ -1039,10 +1198,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); solution.objective = incumbent_.objective; - solution.lower_bound = get_lower_bound(); + solution.lower_bound = lower_bound; solution.nodes_explored = stats_.nodes_explored; solution.simplex_iterations = stats_.total_lp_iters; - return status; + return get_status(); } #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index efe28000bd..1bb67f024c 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -17,7 +17,9 @@ #pragma once +#include #include +#include #include #include #include @@ -25,14 +27,10 @@ #include #include -#include -#include +#include #include #include -#include #include -#include "cuopt/linear_programming/mip/solver_settings.hpp" -#include "dual_simplex/mip_node.hpp" namespace cuopt::linear_programming::dual_simplex { @@ -43,15 +41,53 @@ enum class mip_status_t { TIME_LIMIT = 3, NODE_LIMIT = 4, NUMERICAL = 5, - UNSET = 6 + UNSET = 6, + RUNNING = 7, + FINISHED = 8 }; template void upper_bound_callback(f_t upper_bound); +template +class dive_queue_t { + private: + std::vector>> buffer; + std::mutex mutex; + static constexpr i_t max_size_ = 8192; + + public: + void push(std::shared_ptr>&& node) + { + buffer.push_back(std::forward>>(node)); + std::push_heap(buffer.begin(), buffer.end(), node_compare_t()); + if (buffer.size() > max_size()) { buffer.pop_back(); } + } + + std::shared_ptr> pop() + { + if (buffer.size() == 0) { return std::shared_ptr>(); } + std::shared_ptr> node = buffer.front(); + std::pop_heap(buffer.begin(), buffer.end(), node_compare_t()); + buffer.pop_back(); + return node; + } + + i_t size() const { return buffer.size(); } + constexpr i_t max_size() const { return max_size_; } + std::shared_ptr>& top() { return buffer.front(); } + const std::shared_ptr>& top() const { return buffer.front(); } + void clear() { buffer.clear(); } +}; + +// Note that floating point atomics are only supported in C++20. So, we +// are using omp atomic operations instead. template class branch_and_bound_t { public: + template + using mip_node_heap_t = std::priority_queue, node_compare_t>; + branch_and_bound_t(const user_problem_t& user_problem, const simplex_solver_settings_t& solver_settings); @@ -68,8 +104,9 @@ class branch_and_bound_t { std::vector& repaired_solution) const; f_t get_upper_bound(); - f_t get_lower_bound(); + mip_status_t get_status(); + i_t get_heap_size(); // The main entry routine. Returns the solver status and populates solution with the incumbent. mip_status_t solve(mip_solution_t& solution); @@ -81,16 +118,11 @@ class branch_and_bound_t { // Initial guess. std::vector guess_; + // LP relaxation lp_problem_t original_lp_; std::vector new_slacks_; std::vector var_types_; - // Mutex for lower bound - std::mutex mutex_lower_; - - // Global variable for lower bound - f_t lower_bound_; - // Mutex for upper bound std::mutex mutex_upper_; @@ -100,27 +132,14 @@ class branch_and_bound_t { // Global variable for incumbent. The incumbent should be updated with the upper bound mip_solution_t incumbent_; - // Mutex for gap - std::mutex mutex_gap_; - - // Global variable for gap - f_t gap_; - - // Mutex for branching - std::mutex mutex_branching_; - bool currently_branching_; - - // Global variable for stats - std::mutex mutex_stats_; - - // Note that floating point atomics are only supported in C++20. + // Structure with the general info of the solver. struct stats_t { - f_t start_time = 0.0; - f_t total_lp_solve_time = 0.0; - std::atomic nodes_explored = 0; - std::atomic nodes_unexplored = 0; - f_t total_lp_iters = 0; - std::atomic num_nodes = 0; + f_t start_time = 0.0; + f_t total_lp_solve_time = 0.0; + i_t nodes_explored = 0; + i_t nodes_unexplored = 0; + f_t total_lp_iters = 0; + i_t num_nodes = 0; } stats_; // Mutex for repair @@ -137,6 +156,22 @@ class branch_and_bound_t { pseudo_costs_t pc_; std::mutex mutex_pc_; + // Search tree + std::mutex mutex_search_tree_; + std::unique_ptr> search_tree_; + + // Heap storing the nodes to be explored. + std::mutex mutex_heap_; + mip_node_heap_t*> heap_; + + // Global status + mip_status_t status_; + + i_t active_best_first_tasks_; + + std::mutex mutex_dive_queue_; + dive_queue_t dive_queue_; + // Update the status of the nodes in the search tree. void update_tree(mip_node_t* node_ptr, node_status_t status); @@ -150,11 +185,13 @@ class branch_and_bound_t { // Repairs low-quality solutions from the heuristics, if it is applicable. void repair_heuristic_solutions(); - // Explore the search tree using the best-first search strategy. - mip_status_t explore_tree(i_t branch_var, mip_solution_t& solution); + // Explore the search tree using the best-first search with plunging strategy. + void explore_subtree(mip_node_t* start_node); + void best_first_thread(); - // Explore the search tree using the depth-first search strategy. - mip_status_t dive(i_t branch_var, mip_solution_t& solution); + // Perform a deep dive starting with a given node. + void dive(std::shared_ptr> start_node); + void diving_thread(); // Branch the current node, creating two children. void branch(mip_node_t* parent_node, @@ -166,7 +203,6 @@ class branch_and_bound_t { mip_status_t solve_node_lp(mip_node_t* node_ptr, lp_problem_t& leaf_problem, csc_matrix_t& Arow, - const std::vector& var_types, f_t upper_bound); // Solve the LP relaxation of a leaf node using the dual simplex method. diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 2d315fe0e3..d31efeb86d 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -201,6 +201,18 @@ class mip_node_t { } } + std::shared_ptr> create_detach_copy() const + { + std::shared_ptr> copy = + std::make_shared>(lower_bound, vstatus); + copy->branch_var = branch_var; + copy->branch_dir = branch_dir; + copy->branch_var_lower = branch_var_lower; + copy->branch_var_upper = branch_var_upper; + copy->fractional_val = fractional_val; + return copy; + } + node_status_t status; f_t lower_bound; i_t depth; @@ -230,11 +242,24 @@ void remove_fathomed_nodes(std::vector*>& stack) template class node_compare_t { public: - bool operator()(mip_node_t& a, mip_node_t& b) + bool operator()(mip_node_t& a, mip_node_t& b) const { return a.lower_bound > b.lower_bound; // True if a comes before b, elements that come before are output last } + + bool operator()(mip_node_t* a, mip_node_t* b) const + { + return a->lower_bound > + b->lower_bound; // True if a comes before b, elements that come before are output last + } + + bool operator()(std::shared_ptr> a, + std::shared_ptr> b) const + { + return a->lower_bound > + b->lower_bound; // True if a comes before b, elements that come before are output last + } }; } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/pseudo_costs.cpp b/cpp/src/dual_simplex/pseudo_costs.cpp index 95c96f859f..4d3db6f06b 100644 --- a/cpp/src/dual_simplex/pseudo_costs.cpp +++ b/cpp/src/dual_simplex/pseudo_costs.cpp @@ -21,15 +21,14 @@ #include #include -#include +#include namespace cuopt::linear_programming::dual_simplex { namespace { template -void strong_branch_helper(i_t thread_id, - i_t start, +void strong_branch_helper(i_t start, i_t end, f_t start_time, const lp_problem_t& original_lp, @@ -46,6 +45,7 @@ void strong_branch_helper(i_t thread_id, constexpr bool verbose = false; f_t last_log = tic(); + i_t thread_id = omp_get_thread_num(); for (i_t k = start; k < end; ++k) { const i_t j = fractional[k]; @@ -79,7 +79,7 @@ void strong_branch_helper(i_t thread_id, solution, iter, child_edge_norms); - const f_t lp_solve_time = toc(lp_start_time); + // const f_t lp_solve_time = toc(lp_start_time); f_t obj = std::numeric_limits::infinity(); if (status == dual::status_t::DUAL_UNBOUNDED) { @@ -127,14 +127,11 @@ void strong_branch_helper(i_t thread_id, if (toc(start_time) > settings.time_limit) { break; } pc.strong_branches_completed.lock(); - pc.num_strong_branches_completed++; + const i_t completed = pc.num_strong_branches_completed++; pc.strong_branches_completed.unlock(); if (thread_id == 0 && toc(last_log) > 10) { last_log = tic(); - pc.strong_branches_completed.lock(); - const i_t completed = pc.num_strong_branches_completed; - pc.strong_branches_completed.unlock(); settings.log.printf("%d of %ld strong branches completed in %.1fs\n", completed, fractional.size(), @@ -148,58 +145,6 @@ void strong_branch_helper(i_t thread_id, } } -template -void get_partitioning(i_t N, i_t num_threads, i_t tid, i_t& begin, i_t& end) -{ - int base_tasks = N / num_threads; - int remaining_tasks = N % num_threads; - - if (tid < remaining_tasks) { - begin = tid * (base_tasks + 1); - end = begin + base_tasks + 1; - } else { - begin = tid * base_tasks + remaining_tasks; - end = begin + base_tasks; - } -} - -template -void thread_strong_branch(i_t thread_id, - i_t num_threads, - f_t start_time, - const lp_problem_t& original_lp, - const simplex_solver_settings_t& settings, - const std::vector& var_types, - const std::vector& fractional, - f_t root_obj, - const std::vector& root_soln, - const std::vector& root_vstatus, - const std::vector& edge_norms, - pseudo_costs_t& pc) -{ - i_t start, end; - get_partitioning(static_cast(fractional.size()), num_threads, thread_id, start, end); - constexpr bool verbose = false; - if (verbose) { - settings.log.printf( - "Thread id %d start %d end %d. size %d\n", thread_id, start, end, end - start); - } - - strong_branch_helper(thread_id, - start, - end, - start_time, - original_lp, - settings, - var_types, - fractional, - root_obj, - root_soln, - root_vstatus, - edge_norms, - pc); -} - } // namespace template @@ -214,51 +159,49 @@ void strong_branching(const lp_problem_t& original_lp, const std::vector& edge_norms, pseudo_costs_t& pc) { + pc.resize(original_lp.num_cols); pc.strong_branch_down.resize(fractional.size()); pc.strong_branch_up.resize(fractional.size()); - pc.num_strong_branches_completed = 0; - if (fractional.size() > settings.num_threads) { - int num_threads = settings.num_threads; - std::thread threads[num_threads]; - settings.log.printf("Strong branching using %d threads and %ld fractional variables\n", - num_threads, - fractional.size()); - for (int thread_id = 0; thread_id < num_threads; thread_id++) { - threads[thread_id] = std::thread(thread_strong_branch, - thread_id, - num_threads, - start_time, - original_lp, - settings, - var_types, - fractional, - root_obj, - root_soln, - root_vstatus, - edge_norms, - std::ref(pc)); - } + settings.log.printf("Strong branching using %d threads and %ld fractional variables\n", + settings.num_threads, + fractional.size()); + +#pragma omp parallel num_threads(settings.num_threads) + { + i_t n = std::min(4 * settings.num_threads, fractional.size()); + + // Here we are creating more tasks than the number of threads + // such that they can be scheduled dynamically to the threads. +#pragma omp for schedule(dynamic, 1) + for (i_t k = 0; k < n; k++) { + i_t start = std::floor(k * fractional.size() / n); + i_t end = std::floor((k + 1) * fractional.size() / n); + + constexpr bool verbose = false; + if (verbose) { + settings.log.printf("Thread id %d task id %d start %d end %d. size %d\n", + omp_get_thread_num(), + k, + start, + end, + end - start); + } - for (int thread_id = 0; thread_id < num_threads; thread_id++) { - threads[thread_id].join(); + strong_branch_helper(start, + end, + start_time, + original_lp, + settings, + var_types, + fractional, + root_obj, + root_soln, + root_vstatus, + edge_norms, + pc); } - } else { - settings.log.printf("Strong branching on %ld fractional variables\n", fractional.size()); - strong_branch_helper(0, - 0, - static_cast(fractional.size()), - start_time, - original_lp, - settings, - var_types, - fractional, - root_obj, - root_soln, - root_vstatus, - edge_norms, - pc); } pc.update_pseudo_costs_from_strong_branching(fractional, root_soln); diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 5b7e8bf0fa..3fad1d8ad1 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -59,9 +59,9 @@ struct simplex_solver_settings_t { refactor_frequency(100), iteration_log_frequency(1000), first_iteration_log(2), - num_threads(std::thread::hardware_concurrency() > 8 - ? (std::thread::hardware_concurrency() / 8) - : std::thread::hardware_concurrency()), + num_threads(std::thread::hardware_concurrency() - 1), + num_bfs_threads(std::min(num_threads, 4)), + num_diving_threads(num_threads - num_bfs_threads), random_seed(0), inside_mip(0), solution_callback(nullptr), @@ -113,6 +113,9 @@ struct simplex_solver_settings_t { mutable logger_t log; std::atomic* concurrent_halt; // if nullptr ignored, if !nullptr, 0 if solver should // continue, 1 if solver should halt + + i_t num_bfs_threads; + i_t num_diving_threads; }; } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index a049d0d09b..6400f50cf1 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -90,7 +90,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_ITERATION_LIMIT, &pdlp_settings.iteration_limit, 0, std::numeric_limits::max(), std::numeric_limits::max()}, {CUOPT_PDLP_SOLVER_MODE, reinterpret_cast(&pdlp_settings.pdlp_solver_mode), CUOPT_PDLP_SOLVER_MODE_STABLE1, CUOPT_PDLP_SOLVER_MODE_FAST1, CUOPT_PDLP_SOLVER_MODE_STABLE2}, {CUOPT_METHOD, reinterpret_cast(&pdlp_settings.method), CUOPT_METHOD_CONCURRENT, CUOPT_METHOD_DUAL_SIMPLEX, CUOPT_METHOD_CONCURRENT}, - {CUOPT_NUM_CPU_THREADS, &mip_settings.num_cpu_threads, -1, std::numeric_limits::max(), -1} + {CUOPT_NUM_CPU_THREADS, &mip_settings.num_cpu_threads, -1, std::numeric_limits::max(), -1}, }; // Bool parameters diff --git a/cpp/src/mip/diversity/recombiners/sub_mip.cuh b/cpp/src/mip/diversity/recombiners/sub_mip.cuh index 1818092fc9..de6d917e7c 100644 --- a/cpp/src/mip/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip/diversity/recombiners/sub_mip.cuh @@ -112,6 +112,7 @@ class sub_mip_recombiner_t : public recombiner_t { branch_and_bound_settings.absolute_mip_gap_tol = context.settings.tolerances.absolute_mip_gap; branch_and_bound_settings.relative_mip_gap_tol = context.settings.tolerances.relative_mip_gap; branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance; + branch_and_bound_settings.num_threads = 1; branch_and_bound_settings.solution_callback = [this](std::vector& solution, f_t objective) { this->solution_callback(solution, objective); From 02d4b6042ff5ffbdeda2cc81eea7424bd5c340e5 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 19 Sep 2025 22:26:10 +0200 Subject: [PATCH 02/33] removed unused/duplicated includes --- cpp/src/dual_simplex/branch_and_bound.cpp | 1 - cpp/src/dual_simplex/branch_and_bound.hpp | 1 - 2 files changed, 2 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 5d18f253cb..7a9f912b78 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -15,7 +15,6 @@ * limitations under the License. */ -#include #include #include #include diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 1bb67f024c..a352169cf8 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -17,7 +17,6 @@ #pragma once -#include #include #include #include From bd53b945f1c54ce30a0b294d2dccf029aaa98e17 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 22 Sep 2025 13:49:27 +0200 Subject: [PATCH 03/33] added a separated class for the search tree. code cleanup. --- cpp/src/dual_simplex/branch_and_bound.cpp | 427 ++++++++-------------- cpp/src/dual_simplex/branch_and_bound.hpp | 51 ++- cpp/src/dual_simplex/mip_node.hpp | 103 +++++- 3 files changed, 267 insertions(+), 314 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 7a9f912b78..924e9d5455 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -32,6 +32,7 @@ #include #include #include +#include #include #include @@ -132,37 +133,6 @@ void set_uninitialized_steepest_edge_norms(std::vector& edge_norms) } } -constexpr int write_graphviz = false; - -template -void graphviz_node(const simplex_solver_settings_t& settings, - mip_node_t* node_ptr, - std::string label, - f_t val) -{ - if (write_graphviz) { - settings.log.printf("Node%d [label=\"%s %.16e\"]\n", node_ptr->node_id, label.c_str(), val); - } -} - -template -void graphviz_edge(const simplex_solver_settings_t& settings, - mip_node_t* origin_ptr, - mip_node_t* dest_ptr, - i_t branch_var, - i_t branch_dir, - f_t bound) -{ - if (write_graphviz) { - settings.log.printf("Node%d -> Node%d [label=\"x%d %s %e\"]\n", - origin_ptr->node_id, - dest_ptr->node_id, - branch_var, - branch_dir == 0 ? "<=" : ">=", - bound); - } -} - dual::status_t convert_lp_status_to_dual_status(lp_status_t status) { if (status == lp_status_t::OPTIMAL) { @@ -254,14 +224,10 @@ f_t branch_and_bound_t::get_upper_bound() template f_t branch_and_bound_t::get_lower_bound() { - f_t lower_bound = -inf; + f_t lower_bound = root_objective_; mutex_heap_.lock(); - if (heap_.size() > 0) { - lower_bound = heap_.top()->lower_bound; - } else if (search_tree_ != nullptr) { - lower_bound = search_tree_->lower_bound; - } + if (heap_.size() > 0) { lower_bound = heap_.top()->lower_bound; } mutex_heap_.unlock(); return lower_bound; @@ -287,6 +253,15 @@ i_t branch_and_bound_t::get_heap_size() return size; } +template +i_t branch_and_bound_t::get_active_subtrees() +{ + i_t active_subtrees; +#pragma omp atomic read + active_subtrees = active_subtrees_; + return active_subtrees; +} + template void branch_and_bound_t::set_new_solution(const std::vector& solution) { @@ -463,48 +438,6 @@ void branch_and_bound_t::repair_heuristic_solutions() } } -template -void branch_and_bound_t::branch(mip_node_t* parent_node, - i_t branch_var, - f_t branch_var_val, - const std::vector& parent_vstatus) -{ - i_t id; - -#pragma omp atomic capture - { - id = stats_.num_nodes; - stats_.num_nodes += 2; - } - - // down child - auto down_child = std::make_unique>( - original_lp_, parent_node, ++id, branch_var, 0, branch_var_val, parent_vstatus); - - graphviz_edge( - settings_, parent_node, down_child.get(), branch_var, 0, std::floor(branch_var_val)); - - // up child - auto up_child = std::make_unique>( - original_lp_, parent_node, ++id, branch_var, 1, branch_var_val, parent_vstatus); - - graphviz_edge(settings_, parent_node, up_child.get(), branch_var, 1, std::ceil(branch_var_val)); - - assert(parent_vstatus.size() == original_lp_.num_cols); - parent_node->add_children(std::move(down_child), - std::move(up_child)); // child pointers moved into the tree -} - -template -void branch_and_bound_t::update_tree(mip_node_t* node_ptr, node_status_t status) -{ - std::vector*> stack; - mutex_search_tree_.lock(); - node_ptr->set_status(status, stack); - remove_fathomed_nodes(stack); - mutex_search_tree_.unlock(); -} - template void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, const std::vector& leaf_solution, @@ -550,6 +483,26 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, mutex_upper_.unlock(); } +template +std::pair*, mip_node_t*> +branch_and_bound_t::child_selection(mip_node_t* node_ptr) +{ + // Martin's child selection + const i_t branch_var = node_ptr->get_down_child()->branch_var; + const f_t branch_var_val = node_ptr->get_down_child()->fractional_val; + const f_t down_val = std::floor(root_relax_soln_.x[branch_var]); + const f_t up_val = std::ceil(root_relax_soln_.x[branch_var]); + const f_t down_dist = branch_var_val - down_val; + const f_t up_dist = up_val - branch_var_val; + + if (down_dist < up_dist) { + return std::make_pair(node_ptr->get_down_child(), node_ptr->get_up_child()); + + } else { + return std::make_pair(node_ptr->get_up_child(), node_ptr->get_down_child()); + } +} + template dual::status_t branch_and_bound_t::node_dual_simplex( i_t leaf_id, @@ -606,13 +559,15 @@ dual::status_t branch_and_bound_t::node_dual_simplex( } template -mip_status_t branch_and_bound_t::solve_node_lp(mip_node_t* node_ptr, +mip_status_t branch_and_bound_t::solve_node_lp(search_tree_t& search_tree, + mip_node_t* node_ptr, lp_problem_t& leaf_problem, csc_matrix_t& Arow, - f_t upper_bound) + f_t upper_bound, + logger_t& log) { - logger_t log; - log.log = false; + logger_t pc_log; + pc_log.log = false; f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; f_t rel_fathom_tol = settings_.relative_mip_gap_tol; @@ -636,20 +591,20 @@ mip_status_t branch_and_bound_t::solve_node_lp(mip_node_t* n bounds_changed, Arow, upper_bound, - settings_.log); + log); if (lp_status == dual::status_t::DUAL_UNBOUNDED) { node_ptr->lower_bound = inf; - graphviz_node(settings_, node_ptr, "infeasible", 0.0); - update_tree(node_ptr, node_status_t::INFEASIBLE); + search_tree.graphviz_node(node_ptr, "infeasible", 0.0); + search_tree.update_tree(node_ptr, node_status_t::INFEASIBLE); // Node was infeasible. Do not branch } else if (lp_status == dual::status_t::CUTOFF) { node_ptr->lower_bound = upper_bound; f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); - graphviz_node(settings_, node_ptr, "cut off", leaf_objective); - update_tree(node_ptr, node_status_t::FATHOMED); + search_tree.graphviz_node(node_ptr, "cut off", leaf_objective); + search_tree.update_tree(node_ptr, node_status_t::FATHOMED); // Node was cut off. Do not branch } else if (lp_status == dual::status_t::OPTIMAL) { @@ -658,7 +613,7 @@ mip_status_t branch_and_bound_t::solve_node_lp(mip_node_t* n i_t leaf_num_fractional = fractional_variables(settings_, leaf_solution.x, var_types_, leaf_fractional); f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); - graphviz_node(settings_, node_ptr, "lower bound", leaf_objective); + search_tree.graphviz_node(node_ptr, "lower bound", leaf_objective); mutex_pc_.lock(); pc_.update_pseudo_costs(node_ptr, leaf_objective); @@ -668,34 +623,35 @@ mip_status_t branch_and_bound_t::solve_node_lp(mip_node_t* n if (leaf_num_fractional == 0) { add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, 'B'); - graphviz_node(settings_, node_ptr, "integer feasible", leaf_objective); - update_tree(node_ptr, node_status_t::INTEGER_FEASIBLE); + search_tree.graphviz_node(node_ptr, "integer feasible", leaf_objective); + search_tree.update_tree(node_ptr, node_status_t::INTEGER_FEASIBLE); } else if (leaf_objective <= upper_bound + abs_fathom_tol && relative_gap(upper_bound, leaf_objective) > rel_fathom_tol) { // Choose fractional variable to branch on mutex_pc_.lock(); const i_t branch_var = pc_.variable_selection( - leaf_fractional, leaf_solution.x, leaf_problem.lower, leaf_problem.upper, log); + leaf_fractional, leaf_solution.x, leaf_problem.lower, leaf_problem.upper, pc_log); mutex_pc_.unlock(); assert(leaf_vstatus.size() == leaf_problem.num_cols); - branch(node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus); + search_tree.branch( + node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus, original_lp_); node_ptr->status = node_status_t::HAS_CHILDREN; } else { - graphviz_node(settings_, node_ptr, "fathomed", leaf_objective); - update_tree(node_ptr, node_status_t::FATHOMED); + search_tree.graphviz_node(node_ptr, "fathomed", leaf_objective); + search_tree.update_tree(node_ptr, node_status_t::FATHOMED); } } else if (lp_status == dual::status_t::TIME_LIMIT) { - graphviz_node(settings_, node_ptr, "timeout", 0.0); + search_tree.graphviz_node(node_ptr, "timeout", 0.0); return mip_status_t::TIME_LIMIT; } else { - graphviz_node(settings_, node_ptr, "numerical", 0.0); - settings_.log.printf("Encountered LP status %d on node %d. This indicates a numerical issue.\n", - lp_status, - node_ptr->node_id); + search_tree.graphviz_node(node_ptr, "numerical", 0.0); + log.printf("Encountered LP status %d on node %d. This indicates a numerical issue.\n", + lp_status, + node_ptr->node_id); return mip_status_t::NUMERICAL; } @@ -703,13 +659,11 @@ mip_status_t branch_and_bound_t::solve_node_lp(mip_node_t* n } template -void branch_and_bound_t::explore_subtree(mip_node_t* start_node) +void branch_and_bound_t::explore_subtree(search_tree_t& search_tree, + mip_node_t* start_node, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow) { - // Make a copy of the original LP. We will modify its bounds at each leaf - lp_problem_t leaf_problem = original_lp_; - csc_matrix_t Arow(1, 1, 1); - leaf_problem.A.transpose(Arow); - std::deque*> stack; stack.push_front(start_node); @@ -739,8 +693,8 @@ void branch_and_bound_t::explore_subtree(mip_node_t* start_n if (upper_bound < node_ptr->lower_bound || relative_gap(upper_bound, lower_bound) < settings_.relative_mip_gap_tol) { - graphviz_node(settings_, node_ptr, "cutoff", node_ptr->lower_bound); - update_tree(node_ptr, node_status_t::FATHOMED); + search_tree.graphviz_node(node_ptr, "cutoff", node_ptr->lower_bound); + search_tree.update_tree(node_ptr, node_status_t::FATHOMED); continue; } @@ -777,7 +731,8 @@ void branch_and_bound_t::explore_subtree(mip_node_t* start_n // We are just checking for numerical issues or timeout during the LP solve, otherwise // lp_status is set to RUNNING. - mip_status_t lp_status = solve_node_lp(node_ptr, leaf_problem, Arow, upper_bound); + mip_status_t lp_status = + solve_node_lp(search_tree, node_ptr, leaf_problem, Arow, upper_bound, settings_.log); if (lp_status == mip_status_t::TIME_LIMIT) { #pragma omp atomic write status_ = lp_status; @@ -787,14 +742,6 @@ void branch_and_bound_t::explore_subtree(mip_node_t* start_n } if (node_ptr->status == node_status_t::HAS_CHILDREN) { - // Martin's child selection - const i_t branch_var = node_ptr->get_down_child()->branch_var; - const f_t branch_var_val = node_ptr->get_down_child()->fractional_val; - const f_t down_val = std::floor(root_relax_soln_.x[branch_var]); - const f_t up_val = std::ceil(root_relax_soln_.x[branch_var]); - const f_t down_dist = branch_var_val - down_val; - const f_t up_dist = up_val - branch_var_val; - if (stack.size() > 0) { mip_node_t* node = stack.back(); stack.pop_back(); @@ -804,141 +751,32 @@ void branch_and_bound_t::explore_subtree(mip_node_t* start_n mutex_heap_.unlock(); mutex_dive_queue_.lock(); - dive_queue_.push(node->create_detach_copy()); + dive_queue_.push(node->detach_copy()); mutex_dive_queue_.unlock(); } - if (down_dist < up_dist) { - stack.push_front(node_ptr->get_up_child()); - stack.push_front(node_ptr->get_down_child()); - } else { - stack.push_front(node_ptr->get_down_child()); - stack.push_front(node_ptr->get_up_child()); - } - #pragma omp atomic update stats_.nodes_unexplored += 2; + + auto [first, second] = child_selection(node_ptr); + stack.push_front(second); + stack.push_front(first); } } } template -void branch_and_bound_t::dive(std::shared_ptr> start_node) +void branch_and_bound_t::best_first_thread(search_tree_t& search_tree) { - logger_t log; - log.log = false; - - f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; - f_t rel_fathom_tol = settings_.relative_mip_gap_tol; + f_t lower_bound = -inf; + f_t upper_bound = inf; + f_t gap = inf; // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; csc_matrix_t Arow(1, 1, 1); leaf_problem.A.transpose(Arow); - std::deque*> stack; - stack.push_front(start_node.get()); - - // It contains the lower bound of the parent for now - f_t upper_bound = get_upper_bound(); - i_t num_nodes = 1; - - while (stack.size() > 0 && get_status() == mip_status_t::RUNNING) { - if (stack.size() > 1 && dive_queue_.size() < 4 * settings_.num_diving_threads) { - std::lock_guard lock(mutex_dive_queue_); - mip_node_t* node = stack.back(); - stack.pop_back(); - dive_queue_.push(node->create_detach_copy()); - } - - mip_node_t* node_ptr = stack.front(); - stack.pop_front(); - upper_bound = get_upper_bound(); - - if (upper_bound < node_ptr->lower_bound || - relative_gap(upper_bound, node_ptr->lower_bound) < settings_.relative_mip_gap_tol) { - continue; - } - - std::vector& leaf_vstatus = node_ptr->vstatus; - lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); - - leaf_problem.lower = original_lp_.lower; - leaf_problem.upper = original_lp_.upper; - std::vector bounds_changed(leaf_problem.num_cols, false); - node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); - - dual::status_t lp_status = node_dual_simplex(node_ptr->node_id, - leaf_problem, - leaf_vstatus, - leaf_solution, - bounds_changed, - Arow, - upper_bound, - log); - - if (lp_status != dual::status_t::OPTIMAL) { continue; } - - std::vector fractional; - const i_t leaf_fractional = - fractional_variables(settings_, leaf_solution.x, var_types_, fractional); - f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); - - mutex_pc_.lock(); - pc_.update_pseudo_costs(node_ptr, leaf_objective); - mutex_pc_.unlock(); - - node_ptr->lower_bound = leaf_objective; - - if (leaf_fractional == 0) { - add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, 'D'); - } else if (leaf_objective <= upper_bound + abs_fathom_tol && - relative_gap(upper_bound, leaf_objective) > rel_fathom_tol) { - // Choose fractional variable to branch on - mutex_pc_.lock(); - const i_t branch_var = pc_.variable_selection( - fractional, leaf_solution.x, leaf_problem.lower, leaf_problem.upper, log); - mutex_pc_.unlock(); - const f_t branch_var_val = leaf_solution.x[branch_var]; - - assert(leaf_vstatus.size() == leaf_problem.num_cols); - assert(leaf_vstatus.size() == original_lp_.num_cols); - - // down child - auto down_child = std::make_unique>( - original_lp_, node_ptr, ++num_nodes, branch_var, 0, branch_var_val, leaf_vstatus); - - // up child - auto up_child = std::make_unique>( - original_lp_, node_ptr, ++num_nodes, branch_var, 1, branch_var_val, leaf_vstatus); - - // Martin's child selection - const f_t down_val = std::floor(root_relax_soln_.x[branch_var]); - const f_t up_val = std::ceil(root_relax_soln_.x[branch_var]); - const f_t down_dist = -down_val; - const f_t up_dist = up_val - leaf_solution.x[branch_var]; - - if (down_dist < up_dist) { - stack.push_front(up_child.get()); - stack.push_front(down_child.get()); - } else { - stack.push_front(down_child.get()); - stack.push_front(up_child.get()); - } - - node_ptr->add_children(std::move(down_child), std::move(up_child)); - } - } -} - -template -void branch_and_bound_t::best_first_thread() -{ - i_t active = 0; - f_t lower_bound = -inf; - f_t upper_bound = inf; - f_t gap = inf; - do { mip_node_t* node_ptr = nullptr; @@ -947,6 +785,9 @@ void branch_and_bound_t::best_first_thread() if (heap_.size() > 0) { node_ptr = heap_.top(); heap_.pop(); + +#pragma omp atomic update + active_subtrees_++; } mutex_heap_.unlock(); @@ -954,48 +795,92 @@ void branch_and_bound_t::best_first_thread() if (get_upper_bound() < node_ptr->lower_bound) { // This node was put on the heap earlier but its lower bound is now greater than the // current upper bound - graphviz_node(settings_, node_ptr, "cutoff", node_ptr->lower_bound); - update_tree(node_ptr, node_status_t::FATHOMED); - continue; - } + search_tree.graphviz_node(node_ptr, "cutoff", node_ptr->lower_bound); + search_tree.update_tree(node_ptr, node_status_t::FATHOMED); #pragma omp atomic update - active_best_first_tasks_++; + active_subtrees_--; + continue; + } // Best-first search with plunging - explore_subtree(node_ptr); + explore_subtree(search_tree, node_ptr, leaf_problem, Arow); #pragma omp atomic update - active_best_first_tasks_--; + active_subtrees_--; } -#pragma omp atomic read - active = active_best_first_tasks_; - lower_bound = get_lower_bound(); upper_bound = get_upper_bound(); gap = upper_bound - lower_bound; } while (get_status() == mip_status_t::RUNNING && gap > settings_.absolute_mip_gap_tol && relative_gap(upper_bound, lower_bound) > settings_.relative_mip_gap_tol && - (active > 0 || get_heap_size() > 0)); + (get_active_subtrees() > 0 || get_heap_size() > 0)); } template void branch_and_bound_t::diving_thread() { - while (get_status() == mip_status_t::RUNNING) { - std::shared_ptr> node_ptr; + logger_t log; + log.log = false; + + // Make a copy of the original LP. We will modify its bounds at each leaf + lp_problem_t leaf_problem = original_lp_; + csc_matrix_t Arow(1, 1, 1); + leaf_problem.A.transpose(Arow); + + do { + std::optional> start_node; mutex_dive_queue_.lock(); - if (dive_queue_.size() > 0) { node_ptr = dive_queue_.pop(); } + if (dive_queue_.size() > 0) { start_node = dive_queue_.pop(); } mutex_dive_queue_.unlock(); - if (node_ptr) { - if (get_upper_bound() < node_ptr->lower_bound) { continue; } - dive(node_ptr); + if (start_node) { + if (get_upper_bound() < start_node->lower_bound) { continue; } + + search_tree_t subtree(std::move(start_node.value()), settings_.log); + std::deque*> stack; + stack.push_front(&subtree.root); + + while (stack.size() > 0 && get_status() == mip_status_t::RUNNING) { + if (stack.size() > 1 && dive_queue_.size() < 4 * settings_.num_diving_threads) { + std::lock_guard lock(mutex_dive_queue_); + mip_node_t* node = stack.back(); + stack.pop_back(); + dive_queue_.push(node->detach_copy()); + } + + mip_node_t* node_ptr = stack.front(); + stack.pop_front(); + f_t upper_bound = get_upper_bound(); + + if (upper_bound < node_ptr->lower_bound || + relative_gap(upper_bound, node_ptr->lower_bound) < settings_.relative_mip_gap_tol) { + continue; + } + + mip_status_t lp_status = + solve_node_lp(subtree, node_ptr, leaf_problem, Arow, upper_bound, log); + if (lp_status == mip_status_t::TIME_LIMIT) { break; } + + if (node_ptr->status == node_status_t::HAS_CHILDREN) { + auto [first, second] = child_selection(node_ptr); + + if (dive_queue_.size() < 4 * settings_.num_diving_threads) { + dive_queue_.push(second->detach_copy()); + } else { + stack.push_front(second); + } + + stack.push_front(first); + } + } } - } + + } while (get_status() == mip_status_t::RUNNING && + (get_active_subtrees() > 0 || get_heap_size() > 0)); } template @@ -1006,11 +891,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut log.log = false; status_ = mip_status_t::UNSET; - stats_.total_lp_iters = 0; - stats_.nodes_explored = 0; - stats_.nodes_unexplored = 2; - stats_.num_nodes = 0; - active_best_first_tasks_ = 0; + stats_.total_lp_iters = 0; + stats_.nodes_explored = 0; + stats_.nodes_unexplored = 2; + active_subtrees_ = 0; if (guess_.size() != 0) { std::vector crushed_guess; @@ -1105,6 +989,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::OPTIMAL; } + pc_.resize(original_lp_.num_cols); strong_branching(original_lp_, settings_, stats_.start_time, @@ -1125,12 +1010,13 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut i_t branch_var = pc_.variable_selection( fractional, root_relax_soln_.x, original_lp_.lower, original_lp_.upper, log); - search_tree_ = std::make_unique>(root_objective_, root_vstatus_); - graphviz_node(settings_, search_tree_.get(), "lower bound", root_objective_); - branch(search_tree_.get(), branch_var, root_relax_soln_.x[branch_var], root_vstatus_); + search_tree_t search_tree(root_objective_, root_vstatus_, settings_.log); + search_tree.graphviz_node(&search_tree.root, "lower bound", root_objective_); + search_tree.branch( + &search_tree.root, branch_var, root_relax_soln_.x[branch_var], root_vstatus_, original_lp_); - heap_.push(search_tree_->get_down_child()); - heap_.push(search_tree_->get_up_child()); + heap_.push(search_tree.root.get_down_child()); + heap_.push(search_tree.root.get_up_child()); settings_.log.printf("Exploring the B&B tree using %d threads and %d diving threads\n", settings_.num_bfs_threads, @@ -1145,19 +1031,16 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut #pragma omp parallel num_threads(settings_.num_threads) { if (omp_get_thread_num() < settings_.num_bfs_threads) { - best_first_thread(); + best_first_thread(search_tree); } else { diving_thread(); } + } -// Check if the solver exited normally. -#pragma omp single nowait - { - if (status_ == mip_status_t::RUNNING) { + // Check if the solver exited naturally, or was due to a timeout or numerical error. + if (status_ == mip_status_t::RUNNING) { #pragma omp atomic write - status_ = mip_status_t::FINISHED; - } - } + status_ = mip_status_t::FINISHED; } f_t lower_bound = get_lower_bound(); @@ -1200,7 +1083,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut solution.lower_bound = lower_bound; solution.nodes_explored = stats_.nodes_explored; solution.simplex_iterations = stats_.total_lp_iters; - return get_status(); + return status_; } #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index a352169cf8..8299dc02c4 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -51,31 +51,31 @@ void upper_bound_callback(f_t upper_bound); template class dive_queue_t { private: - std::vector>> buffer; + std::vector> buffer; std::mutex mutex; static constexpr i_t max_size_ = 8192; public: - void push(std::shared_ptr>&& node) + dive_queue_t() { buffer.reserve(max_size_); } + + void push(mip_node_t&& node) { - buffer.push_back(std::forward>>(node)); + buffer.push_back(std::forward&&>(node)); std::push_heap(buffer.begin(), buffer.end(), node_compare_t()); if (buffer.size() > max_size()) { buffer.pop_back(); } } - std::shared_ptr> pop() + mip_node_t pop() { - if (buffer.size() == 0) { return std::shared_ptr>(); } - std::shared_ptr> node = buffer.front(); std::pop_heap(buffer.begin(), buffer.end(), node_compare_t()); + mip_node_t node = std::move(buffer.back()); buffer.pop_back(); return node; } i_t size() const { return buffer.size(); } constexpr i_t max_size() const { return max_size_; } - std::shared_ptr>& top() { return buffer.front(); } - const std::shared_ptr>& top() const { return buffer.front(); } + const mip_node_t& top() const { return buffer.front(); } void clear() { buffer.clear(); } }; @@ -106,6 +106,7 @@ class branch_and_bound_t { f_t get_lower_bound(); mip_status_t get_status(); i_t get_heap_size(); + i_t get_active_subtrees(); // The main entry routine. Returns the solver status and populates solution with the incumbent. mip_status_t solve(mip_solution_t& solution); @@ -138,7 +139,6 @@ class branch_and_bound_t { i_t nodes_explored = 0; i_t nodes_unexplored = 0; f_t total_lp_iters = 0; - i_t num_nodes = 0; } stats_; // Mutex for repair @@ -155,10 +155,6 @@ class branch_and_bound_t { pseudo_costs_t pc_; std::mutex mutex_pc_; - // Search tree - std::mutex mutex_search_tree_; - std::unique_ptr> search_tree_; - // Heap storing the nodes to be explored. std::mutex mutex_heap_; mip_node_heap_t*> heap_; @@ -166,14 +162,11 @@ class branch_and_bound_t { // Global status mip_status_t status_; - i_t active_best_first_tasks_; + i_t active_subtrees_; std::mutex mutex_dive_queue_; dive_queue_t dive_queue_; - // Update the status of the nodes in the search tree. - void update_tree(mip_node_t* node_ptr, node_status_t status); - // Update the incumbent solution with the new feasible solution. // found during branch and bound. void add_feasible_solution(f_t leaf_objective, @@ -185,24 +178,21 @@ class branch_and_bound_t { void repair_heuristic_solutions(); // Explore the search tree using the best-first search with plunging strategy. - void explore_subtree(mip_node_t* start_node); - void best_first_thread(); + void explore_subtree(search_tree_t& search_tree, + mip_node_t* start_node, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow); + void best_first_thread(search_tree_t& search_tree); - // Perform a deep dive starting with a given node. - void dive(std::shared_ptr> start_node); void diving_thread(); - // Branch the current node, creating two children. - void branch(mip_node_t* parent_node, - i_t branch_var, - f_t branch_var_val, - const std::vector& parent_vstatus); - // Solve the LP relaxation of a leaf node. - mip_status_t solve_node_lp(mip_node_t* node_ptr, + mip_status_t solve_node_lp(search_tree_t& search_tree, + mip_node_t* node_ptr, lp_problem_t& leaf_problem, csc_matrix_t& Arow, - f_t upper_bound); + f_t upper_bound, + logger_t& log); // Solve the LP relaxation of a leaf node using the dual simplex method. dual::status_t node_dual_simplex(i_t leaf_id, @@ -213,6 +203,9 @@ class branch_and_bound_t { csc_matrix_t& Arow, f_t upper_bound, logger_t& log); + + std::pair*, mip_node_t*> child_selection( + mip_node_t* node_ptr); }; } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index d31efeb86d..d59bf82f6e 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -20,10 +20,10 @@ #include #include -#include #include #include #include +#include #include namespace cuopt::linear_programming::dual_simplex { @@ -201,15 +201,14 @@ class mip_node_t { } } - std::shared_ptr> create_detach_copy() const + mip_node_t detach_copy() const { - std::shared_ptr> copy = - std::make_shared>(lower_bound, vstatus); - copy->branch_var = branch_var; - copy->branch_dir = branch_dir; - copy->branch_var_lower = branch_var_lower; - copy->branch_var_upper = branch_var_upper; - copy->fractional_val = fractional_val; + mip_node_t copy(lower_bound, vstatus); + copy.branch_var = branch_var; + copy.branch_dir = branch_dir; + copy.branch_var_lower = branch_var_lower; + copy.branch_var_upper = branch_var_upper; + copy.fractional_val = fractional_val; return copy; } @@ -253,13 +252,91 @@ class node_compare_t { return a->lower_bound > b->lower_bound; // True if a comes before b, elements that come before are output last } +}; - bool operator()(std::shared_ptr> a, - std::shared_ptr> b) const +template +class search_tree_t { + public: + search_tree_t(f_t root_lower_bound, const std::vector& basis, logger_t& log) + : root(root_lower_bound, basis), num_nodes(0), log(log) { - return a->lower_bound > - b->lower_bound; // True if a comes before b, elements that come before are output last } + + search_tree_t(mip_node_t&& node, logger_t& log) + : root(std::forward&&>(node)), num_nodes(0), log(log) + { + } + + f_t get_lower_bound() const { return root.lower_bound; } + + void update_tree(mip_node_t* node_ptr, node_status_t status) + { + std::lock_guard lock(mutex); + std::vector*> stack; + node_ptr->set_status(status, stack); + remove_fathomed_nodes(stack); + } + + void branch(mip_node_t* parent_node, + i_t branch_var, + f_t branch_var_val, + const std::vector& parent_vstatus, + const lp_problem_t& original_lp) + { + i_t id; + +#pragma omp atomic capture + { + id = num_nodes; + num_nodes += 2; + } + + // down child + auto down_child = std::make_unique>( + original_lp, parent_node, ++id, branch_var, 0, branch_var_val, parent_vstatus); + + graphviz_edge(parent_node, down_child.get(), branch_var, 0, std::floor(branch_var_val)); + + // up child + auto up_child = std::make_unique>( + original_lp, parent_node, ++id, branch_var, 1, branch_var_val, parent_vstatus); + + graphviz_edge(parent_node, up_child.get(), branch_var, 1, std::ceil(branch_var_val)); + + assert(parent_vstatus.size() == original_lp.num_cols); + parent_node->add_children(std::move(down_child), + std::move(up_child)); // child pointers moved into the tree + } + + void graphviz_node(mip_node_t* node_ptr, std::string label, f_t val) + { + if (write_graphviz) { + log.printf("Node%d [label=\"%s %.16e\"]\n", node_ptr->node_id, label.c_str(), val); + } + } + + void graphviz_edge(mip_node_t* origin_ptr, + mip_node_t* dest_ptr, + i_t branch_var, + i_t branch_dir, + f_t bound) + { + if (write_graphviz) { + log.printf("Node%d -> Node%d [label=\"x%d %s %e\"]\n", + origin_ptr->node_id, + dest_ptr->node_id, + branch_var, + branch_dir == 0 ? "<=" : ">=", + bound); + } + } + + mip_node_t root; + std::mutex mutex; + i_t num_nodes; + logger_t log; + + static constexpr int write_graphviz = false; }; } // namespace cuopt::linear_programming::dual_simplex From 3e6c21c8bbec79c78af21bc3707a0424552e9e81 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 22 Sep 2025 14:12:48 +0200 Subject: [PATCH 04/33] fixed race condition --- cpp/src/dual_simplex/branch_and_bound.cpp | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 924e9d5455..e0b832ae6f 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -750,9 +750,11 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear heap_.push(node); mutex_heap_.unlock(); - mutex_dive_queue_.lock(); - dive_queue_.push(node->detach_copy()); - mutex_dive_queue_.unlock(); + if (get_heap_size() > 4 * settings_.num_bfs_threads) { + mutex_dive_queue_.lock(); + dive_queue_.push(node->detach_copy()); + mutex_dive_queue_.unlock(); + } } #pragma omp atomic update @@ -845,13 +847,6 @@ void branch_and_bound_t::diving_thread() stack.push_front(&subtree.root); while (stack.size() > 0 && get_status() == mip_status_t::RUNNING) { - if (stack.size() > 1 && dive_queue_.size() < 4 * settings_.num_diving_threads) { - std::lock_guard lock(mutex_dive_queue_); - mip_node_t* node = stack.back(); - stack.pop_back(); - dive_queue_.push(node->detach_copy()); - } - mip_node_t* node_ptr = stack.front(); stack.pop_front(); f_t upper_bound = get_upper_bound(); @@ -869,7 +864,9 @@ void branch_and_bound_t::diving_thread() auto [first, second] = child_selection(node_ptr); if (dive_queue_.size() < 4 * settings_.num_diving_threads) { + mutex_dive_queue_.lock(); dive_queue_.push(second->detach_copy()); + mutex_dive_queue_.unlock(); } else { stack.push_front(second); } From 3c9b1921dba906eecb0c7e06bd78077e22a6baa5 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 22 Sep 2025 15:35:49 +0200 Subject: [PATCH 05/33] added ramp up phase --- cpp/src/dual_simplex/branch_and_bound.cpp | 205 ++++++++++++++++------ cpp/src/dual_simplex/branch_and_bound.hpp | 27 ++- 2 files changed, 174 insertions(+), 58 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index e0b832ae6f..6fa1f23c96 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -219,17 +219,13 @@ f_t branch_and_bound_t::get_upper_bound() return upper_bound; } -// Calculates the global lower bound considering -// the lower bounds in each thread and the top of the heap. template f_t branch_and_bound_t::get_lower_bound() { f_t lower_bound = root_objective_; - mutex_heap_.lock(); if (heap_.size() > 0) { lower_bound = heap_.top()->lower_bound; } mutex_heap_.unlock(); - return lower_bound; } @@ -487,7 +483,6 @@ template std::pair*, mip_node_t*> branch_and_bound_t::child_selection(mip_node_t* node_ptr) { - // Martin's child selection const i_t branch_var = node_ptr->get_down_child()->branch_var; const f_t branch_var_val = node_ptr->get_down_child()->fractional_val; const f_t down_val = std::floor(root_relax_soln_.x[branch_var]); @@ -658,6 +653,98 @@ mip_status_t branch_and_bound_t::solve_node_lp(search_tree_t return mip_status_t::RUNNING; } +template +void branch_and_bound_t::exploration_ramp_up(search_tree_t* search_tree, + mip_node_t* node, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow, + i_t max_depth) +{ + if (get_status() != mip_status_t::RUNNING) { return; } + repair_heuristic_solutions(); + + f_t lower_bound = node->lower_bound; + f_t upper_bound = get_upper_bound(); + f_t gap = upper_bound - lower_bound; + i_t nodes_explored = 0; + i_t nodes_unexplored = 0; + +#pragma omp atomic capture + nodes_explored = stats_.nodes_explored++; + +#pragma omp atomic capture + nodes_unexplored = stats_.nodes_unexplored--; + + if (upper_bound < node->lower_bound || + relative_gap(upper_bound, lower_bound) < settings_.relative_mip_gap_tol) { + search_tree->graphviz_node(node, "cutoff", node->lower_bound); + search_tree->update_tree(node, node_status_t::FATHOMED); + return; + } + + f_t now = toc(stats_.start_time); + +#pragma omp single nowait + { + if (nodes_explored % 1000 == 0 || gap < 10 * settings_.absolute_mip_gap_tol || + nodes_explored < 1000) { + f_t obj = compute_user_objective(original_lp_, upper_bound); + f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); + std::string gap_user = user_mip_gap(obj, user_lower); + + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + nodes_explored, + nodes_unexplored, + obj, + user_lower, + node->depth, + nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, + gap_user.c_str(), + now); + } + } + + if (toc(stats_.start_time) > settings_.time_limit) { +#pragma omp atomic write + status_ = mip_status_t::TIME_LIMIT; + +#pragma omp single nowait + settings_.log.printf("Time limit reached. Stopping the solver...\n"); + return; + } + + // We are just checking for numerical issues or timeout during the LP solve, otherwise + // lp_status is set to RUNNING. + mip_status_t lp_status = + solve_node_lp(*search_tree, node, leaf_problem, Arow, upper_bound, settings_.log); + if (lp_status == mip_status_t::TIME_LIMIT) { +#pragma omp atomic write + status_ = lp_status; +#pragma omp single nowait + settings_.log.printf("Time limit reached. Stopping the solver...\n"); + return; + } + + if (node->status == node_status_t::HAS_CHILDREN) { +#pragma omp atomic update + stats_.nodes_unexplored += 2; + + if (node->depth < max_depth) { +#pragma omp task + exploration_ramp_up(search_tree, node->get_down_child(), leaf_problem, Arow, max_depth); + +#pragma omp task + exploration_ramp_up(search_tree, node->get_up_child(), leaf_problem, Arow, max_depth); + + } else { + mutex_heap_.lock(); + heap_.push(node->get_down_child()); + heap_.push(node->get_up_child()); + mutex_heap_.unlock(); + } + } +} + template void branch_and_bound_t::explore_subtree(search_tree_t& search_tree, mip_node_t* start_node, @@ -700,25 +787,28 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear f_t now = toc(stats_.start_time); - f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); - if (((nodes_explored % 1000 == 0 || gap < 10 * settings_.absolute_mip_gap_tol || - nodes_explored < 1000) && - (time_since_log >= 1)) || - (time_since_log > 60) || now > settings_.time_limit) { - f_t obj = compute_user_objective(original_lp_, upper_bound); - f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); - std::string gap_user = user_mip_gap(obj, user_lower); - - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - nodes_explored, - nodes_unexplored, - obj, - user_lower, - node_ptr->depth, - nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, - gap_user.c_str(), - now); - last_log = tic(); +#pragma omp single nowait + { + f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); + if (((nodes_explored % 1000 == 0 || gap < 10 * settings_.absolute_mip_gap_tol || + nodes_explored < 1000) && + (time_since_log >= 1)) || + (time_since_log > 60) || now > settings_.time_limit) { + f_t obj = compute_user_objective(original_lp_, upper_bound); + f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); + std::string gap_user = user_mip_gap(obj, user_lower); + + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + nodes_explored, + nodes_unexplored, + obj, + user_lower, + node_ptr->depth, + nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, + gap_user.c_str(), + now); + last_log = tic(); + } } if (toc(stats_.start_time) > settings_.time_limit) { @@ -750,11 +840,9 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear heap_.push(node); mutex_heap_.unlock(); - if (get_heap_size() > 4 * settings_.num_bfs_threads) { - mutex_dive_queue_.lock(); - dive_queue_.push(node->detach_copy()); - mutex_dive_queue_.unlock(); - } + mutex_dive_queue_.lock(); + dive_queue_.push(node->detach_copy()); + mutex_dive_queue_.unlock(); } #pragma omp atomic update @@ -768,17 +856,14 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear } template -void branch_and_bound_t::best_first_thread(search_tree_t& search_tree) +void branch_and_bound_t::best_first_thread(search_tree_t& search_tree, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow) { f_t lower_bound = -inf; f_t upper_bound = inf; f_t gap = inf; - // Make a copy of the original LP. We will modify its bounds at each leaf - lp_problem_t leaf_problem = original_lp_; - csc_matrix_t Arow(1, 1, 1); - leaf_problem.A.transpose(Arow); - do { mip_node_t* node_ptr = nullptr; @@ -819,19 +904,21 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se } while (get_status() == mip_status_t::RUNNING && gap > settings_.absolute_mip_gap_tol && relative_gap(upper_bound, lower_bound) > settings_.relative_mip_gap_tol && (get_active_subtrees() > 0 || get_heap_size() > 0)); + + // Check if the solver exited naturally, or was due to a timeout or numerical error. + if (status_ == mip_status_t::RUNNING) { +#pragma omp atomic write + status_ = mip_status_t::FINISHED; + } } template -void branch_and_bound_t::diving_thread() +void branch_and_bound_t::diving_thread(lp_problem_t& leaf_problem, + csc_matrix_t& Arow) { logger_t log; log.log = false; - // Make a copy of the original LP. We will modify its bounds at each leaf - lp_problem_t leaf_problem = original_lp_; - csc_matrix_t Arow(1, 1, 1); - leaf_problem.A.transpose(Arow); - do { std::optional> start_node; @@ -884,7 +971,6 @@ template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { logger_t log; - log.log = false; status_ = mip_status_t::UNSET; @@ -1012,9 +1098,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut search_tree.branch( &search_tree.root, branch_var, root_relax_soln_.x[branch_var], root_vstatus_, original_lp_); - heap_.push(search_tree.root.get_down_child()); - heap_.push(search_tree.root.get_up_child()); - settings_.log.printf("Exploring the B&B tree using %d threads and %d diving threads\n", settings_.num_bfs_threads, settings_.num_threads - settings_.num_bfs_threads); @@ -1027,20 +1110,34 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut #pragma omp parallel num_threads(settings_.num_threads) { + // Make a copy of the original LP. We will modify its bounds at each leaf + lp_problem_t leaf_problem = original_lp_; + csc_matrix_t Arow(1, 1, 1); + leaf_problem.A.transpose(Arow); + +#pragma omp master + { + i_t max_depth = 4 * std::ceil(log2((double)settings_.num_bfs_threads)); + +#pragma omp task + exploration_ramp_up( + &search_tree, search_tree.root.get_down_child(), leaf_problem, Arow, max_depth); + +#pragma omp task + exploration_ramp_up( + &search_tree, search_tree.root.get_up_child(), leaf_problem, Arow, max_depth); + } + +#pragma omp barrier + if (omp_get_thread_num() < settings_.num_bfs_threads) { - best_first_thread(search_tree); + best_first_thread(search_tree, leaf_problem, Arow); } else { - diving_thread(); + diving_thread(leaf_problem, Arow); } } - // Check if the solver exited naturally, or was due to a timeout or numerical error. - if (status_ == mip_status_t::RUNNING) { -#pragma omp atomic write - status_ = mip_status_t::FINISHED; - } - - f_t lower_bound = get_lower_bound(); + f_t lower_bound = heap_.size() > 0 ? heap_.top()->lower_bound : search_tree.get_lower_bound(); f_t upper_bound = get_upper_bound(); f_t gap = upper_bound - lower_bound; f_t obj = compute_user_objective(original_lp_, upper_bound); diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 8299dc02c4..02311a1687 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -48,6 +48,9 @@ enum class mip_status_t { template void upper_bound_callback(f_t upper_bound); +// A min-heap for storing the starting nodes for the dives. +// This has a maximum size of 8192, such that the container +// will discard the least promising node if the queue is full. template class dive_queue_t { private: @@ -79,8 +82,6 @@ class dive_queue_t { void clear() { buffer.clear(); } }; -// Note that floating point atomics are only supported in C++20. So, we -// are using omp atomic operations instead. template class branch_and_bound_t { public: @@ -162,8 +163,10 @@ class branch_and_bound_t { // Global status mip_status_t status_; + // Count the number of subtrees that are currently being explored. i_t active_subtrees_; + // Queue for storing the promising node for performing dives. std::mutex mutex_dive_queue_; dive_queue_t dive_queue_; @@ -177,14 +180,29 @@ class branch_and_bound_t { // Repairs low-quality solutions from the heuristics, if it is applicable. void repair_heuristic_solutions(); + // Ramp-up phase of the solver, where we greedily expand the tree until + // a certain depth is reached. This is done recursively using OpenMP tasks. + void exploration_ramp_up(search_tree_t* search_tree, + mip_node_t* node, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow, + i_t max_depth); + // Explore the search tree using the best-first search with plunging strategy. void explore_subtree(search_tree_t& search_tree, mip_node_t* start_node, lp_problem_t& leaf_problem, csc_matrix_t& Arow); - void best_first_thread(search_tree_t& search_tree); - void diving_thread(); + // Each "main" thread pop a node from the global heap and then perform a plunge + // (i.e., a shallow dive) into the subtree determined by the node. + void best_first_thread(search_tree_t& search_tree, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow); + + // Each diving thread pop the first node from the dive queue and then perform + // a deep dive into the subtree determined by the node. + void diving_thread(lp_problem_t& leaf_problem, csc_matrix_t& Arow); // Solve the LP relaxation of a leaf node. mip_status_t solve_node_lp(search_tree_t& search_tree, @@ -204,6 +222,7 @@ class branch_and_bound_t { f_t upper_bound, logger_t& log); + // Sort the children based on the Martin's criteria. std::pair*, mip_node_t*> child_selection( mip_node_t* node_ptr); }; From 189b0232b451ce4867adb142a8c1ab19769cf69f Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 23 Sep 2025 12:16:23 +0200 Subject: [PATCH 06/33] fixed invalid memory access --- cpp/src/dual_simplex/branch_and_bound.cpp | 127 +++++++++++----------- cpp/src/dual_simplex/branch_and_bound.hpp | 13 ++- cpp/src/dual_simplex/mip_node.cpp | 2 +- cpp/src/dual_simplex/mip_node.hpp | 11 +- 4 files changed, 80 insertions(+), 73 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 6fa1f23c96..0b7a4bd28b 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -554,12 +554,13 @@ dual::status_t branch_and_bound_t::node_dual_simplex( } template -mip_status_t branch_and_bound_t::solve_node_lp(search_tree_t& search_tree, - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log) +node_status_t branch_and_bound_t::solve_node_lp(search_tree_t& search_tree, + mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow, + f_t upper_bound, + logger_t& log, + char symbol) { logger_t pc_log; pc_log.log = false; @@ -589,19 +590,20 @@ mip_status_t branch_and_bound_t::solve_node_lp(search_tree_t log); if (lp_status == dual::status_t::DUAL_UNBOUNDED) { + // Node was infeasible. Do not branch node_ptr->lower_bound = inf; search_tree.graphviz_node(node_ptr, "infeasible", 0.0); search_tree.update_tree(node_ptr, node_status_t::INFEASIBLE); - - // Node was infeasible. Do not branch + return node_status_t::INFEASIBLE; } else if (lp_status == dual::status_t::CUTOFF) { + // Node was cut off. Do not branch node_ptr->lower_bound = upper_bound; f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); search_tree.graphviz_node(node_ptr, "cut off", leaf_objective); search_tree.update_tree(node_ptr, node_status_t::FATHOMED); + return node_status_t::FATHOMED; - // Node was cut off. Do not branch } else if (lp_status == dual::status_t::OPTIMAL) { // LP was feasible std::vector leaf_fractional; @@ -617,9 +619,11 @@ mip_status_t branch_and_bound_t::solve_node_lp(search_tree_t node_ptr->lower_bound = leaf_objective; if (leaf_num_fractional == 0) { - add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, 'B'); + // Found a integer feasible solution + add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, symbol); search_tree.graphviz_node(node_ptr, "integer feasible", leaf_objective); search_tree.update_tree(node_ptr, node_status_t::INTEGER_FEASIBLE); + return node_status_t::INTEGER_FEASIBLE; } else if (leaf_objective <= upper_bound + abs_fathom_tol && relative_gap(upper_bound, leaf_objective) > rel_fathom_tol) { @@ -633,24 +637,26 @@ mip_status_t branch_and_bound_t::solve_node_lp(search_tree_t search_tree.branch( node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus, original_lp_); node_ptr->status = node_status_t::HAS_CHILDREN; + return node_status_t::HAS_CHILDREN; } else { search_tree.graphviz_node(node_ptr, "fathomed", leaf_objective); search_tree.update_tree(node_ptr, node_status_t::FATHOMED); + return node_status_t::FATHOMED; } } else if (lp_status == dual::status_t::TIME_LIMIT) { search_tree.graphviz_node(node_ptr, "timeout", 0.0); - return mip_status_t::TIME_LIMIT; + search_tree.update_tree(node_ptr, node_status_t::TIME_LIMIT); + return node_status_t::TIME_LIMIT; } else { search_tree.graphviz_node(node_ptr, "numerical", 0.0); log.printf("Encountered LP status %d on node %d. This indicates a numerical issue.\n", lp_status, node_ptr->node_id); - return mip_status_t::NUMERICAL; + search_tree.update_tree(node_ptr, node_status_t::NUMERICAL); + return node_status_t::NUMERICAL; } - - return mip_status_t::RUNNING; } template @@ -707,25 +713,17 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* if (toc(stats_.start_time) > settings_.time_limit) { #pragma omp atomic write status_ = mip_status_t::TIME_LIMIT; - -#pragma omp single nowait - settings_.log.printf("Time limit reached. Stopping the solver...\n"); return; } + node_status_t node_status = + solve_node_lp(*search_tree, node, leaf_problem, Arow, upper_bound, settings_.log, 'B'); - // We are just checking for numerical issues or timeout during the LP solve, otherwise - // lp_status is set to RUNNING. - mip_status_t lp_status = - solve_node_lp(*search_tree, node, leaf_problem, Arow, upper_bound, settings_.log); - if (lp_status == mip_status_t::TIME_LIMIT) { + if (node_status == node_status_t::TIME_LIMIT) { #pragma omp atomic write - status_ = lp_status; -#pragma omp single nowait - settings_.log.printf("Time limit reached. Stopping the solver...\n"); + status_ = mip_status_t::TIME_LIMIT; return; - } - if (node->status == node_status_t::HAS_CHILDREN) { + } else if (node_status == node_status_t::HAS_CHILDREN) { #pragma omp atomic update stats_.nodes_unexplored += 2; @@ -814,35 +812,33 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear if (toc(stats_.start_time) > settings_.time_limit) { #pragma omp atomic write status_ = mip_status_t::TIME_LIMIT; -#pragma omp single nowait - settings_.log.printf("Time limit reached. Stopping the solver...\n"); break; } // We are just checking for numerical issues or timeout during the LP solve, otherwise // lp_status is set to RUNNING. - mip_status_t lp_status = - solve_node_lp(search_tree, node_ptr, leaf_problem, Arow, upper_bound, settings_.log); - if (lp_status == mip_status_t::TIME_LIMIT) { + node_status_t node_status = + solve_node_lp(search_tree, node_ptr, leaf_problem, Arow, upper_bound, settings_.log, 'B'); + + if (node_status == node_status_t::TIME_LIMIT) { #pragma omp atomic write - status_ = lp_status; -#pragma omp single nowait - settings_.log.printf("Time limit reached. Stopping the solver...\n"); - break; - } + status_ = mip_status_t::TIME_LIMIT; + return; - if (node_ptr->status == node_status_t::HAS_CHILDREN) { + } else if (node_status == node_status_t::HAS_CHILDREN) { if (stack.size() > 0) { mip_node_t* node = stack.back(); stack.pop_back(); + if (get_heap_size() > settings_.num_bfs_threads) { + mutex_dive_queue_.lock(); + dive_queue_.push(node->detach_copy()); + mutex_dive_queue_.unlock(); + } + mutex_heap_.lock(); heap_.push(node); mutex_heap_.unlock(); - - mutex_dive_queue_.lock(); - dive_queue_.push(node->detach_copy()); - mutex_dive_queue_.unlock(); } #pragma omp atomic update @@ -943,11 +939,13 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr continue; } - mip_status_t lp_status = - solve_node_lp(subtree, node_ptr, leaf_problem, Arow, upper_bound, log); - if (lp_status == mip_status_t::TIME_LIMIT) { break; } + node_status_t node_status = + solve_node_lp(subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); + + if (node_status == node_status_t::TIME_LIMIT) { + break; - if (node_ptr->status == node_status_t::HAS_CHILDREN) { + } else if (node_status == node_status_t::HAS_CHILDREN) { auto [first, second] = child_selection(node_ptr); if (dive_queue_.size() < 4 * settings_.num_diving_threads) { @@ -1023,7 +1021,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::UNBOUNDED; } if (root_status == lp_status_t::TIME_LIMIT) { - settings_.log.printf("Hit time limit\n"); + settings_.log.printf("Time limit reached. Stopping the solver...\n"); return mip_status_t::TIME_LIMIT; } @@ -1085,7 +1083,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut pc_); if (toc(stats_.start_time) > settings_.time_limit) { - settings_.log.printf("Hit time limit\n"); + settings_.log.printf("Time limit reached. Stopping the solver...\n"); return mip_status_t::TIME_LIMIT; } @@ -1098,8 +1096,11 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut search_tree.branch( &search_tree.root, branch_var, root_relax_soln_.x[branch_var], root_vstatus_, original_lp_); - settings_.log.printf("Exploring the B&B tree using %d threads and %d diving threads\n", - settings_.num_bfs_threads, + heap_.push(search_tree.root.get_down_child()); + heap_.push(search_tree.root.get_up_child()); + + settings_.log.printf("Exploring the B&B tree using %d threads (%d are diving threads)\n", + settings_.num_threads, settings_.num_threads - settings_.num_bfs_threads); settings_.log.printf( "| Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap | " @@ -1115,20 +1116,20 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut csc_matrix_t Arow(1, 1, 1); leaf_problem.A.transpose(Arow); -#pragma omp master - { - i_t max_depth = 4 * std::ceil(log2((double)settings_.num_bfs_threads)); + // #pragma omp master + // { + // i_t max_depth = std::ceil(std::log2(settings_.num_threads)); -#pragma omp task - exploration_ramp_up( - &search_tree, search_tree.root.get_down_child(), leaf_problem, Arow, max_depth); + // #pragma omp task + // exploration_ramp_up( + // &search_tree, search_tree.root.get_down_child(), leaf_problem, Arow, max_depth); -#pragma omp task - exploration_ramp_up( - &search_tree, search_tree.root.get_up_child(), leaf_problem, Arow, max_depth); - } + // #pragma omp task + // exploration_ramp_up( + // &search_tree, search_tree.root.get_up_child(), leaf_problem, Arow, max_depth); + // } -#pragma omp barrier + // #pragma omp barrier if (omp_get_thread_num() < settings_.num_bfs_threads) { best_first_thread(search_tree, leaf_problem, Arow); @@ -1137,6 +1138,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } } + if (status_ == mip_status_t::TIME_LIMIT) { + settings_.log.printf("Time limit reached. Stopping the solver...\n"); + } + f_t lower_bound = heap_.size() > 0 ? heap_.top()->lower_bound : search_tree.get_lower_bound(); f_t upper_bound = get_upper_bound(); f_t gap = upper_bound - lower_bound; diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 02311a1687..61de2839c9 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -205,12 +205,13 @@ class branch_and_bound_t { void diving_thread(lp_problem_t& leaf_problem, csc_matrix_t& Arow); // Solve the LP relaxation of a leaf node. - mip_status_t solve_node_lp(search_tree_t& search_tree, - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log); + node_status_t solve_node_lp(search_tree_t& search_tree, + mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow, + f_t upper_bound, + logger_t& log, + char symbol); // Solve the LP relaxation of a leaf node using the dual simplex method. dual::status_t node_dual_simplex(i_t leaf_id, diff --git a/cpp/src/dual_simplex/mip_node.cpp b/cpp/src/dual_simplex/mip_node.cpp index f0b275d0a4..2a907c4022 100644 --- a/cpp/src/dual_simplex/mip_node.cpp +++ b/cpp/src/dual_simplex/mip_node.cpp @@ -22,7 +22,7 @@ namespace cuopt::linear_programming::dual_simplex { bool inactive_status(node_status_t status) { return (status == node_status_t::FATHOMED || status == node_status_t::INTEGER_FEASIBLE || - status == node_status_t::INFEASIBLE); + status == node_status_t::INFEASIBLE || status == node_status_t::NUMERICAL); } } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index d59bf82f6e..ee76223f4d 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -30,11 +30,12 @@ namespace cuopt::linear_programming::dual_simplex { enum class node_status_t : int { ACTIVE = 0, // Node still in the tree - IN_PROGRESS = 1, // Node is currently being solved - INTEGER_FEASIBLE = 2, // Node has an integer feasible solution - INFEASIBLE = 3, // Node is infeasible - FATHOMED = 4, // Node objective is greater than the upper bound - HAS_CHILDREN = 5, // Node has children to explore + INTEGER_FEASIBLE = 1, // Node has an integer feasible solution + INFEASIBLE = 2, // Node is infeasible + FATHOMED = 3, // Node objective is greater than the upper bound + HAS_CHILDREN = 4, // Node has children to explore + NUMERICAL = 5, // Encountered numerical issue when solving the LP relaxation + TIME_LIMIT = 6 // Time out during the LP relaxation }; bool inactive_status(node_status_t status); From 065f1124df05d1454c1744f54bc8a080845deea7 Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 23 Sep 2025 14:15:59 +0200 Subject: [PATCH 07/33] re-enabled ramp up phase --- cpp/src/dual_simplex/branch_and_bound.cpp | 25 ++++++++++------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 0b7a4bd28b..c759d4a184 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -1096,9 +1096,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut search_tree.branch( &search_tree.root, branch_var, root_relax_soln_.x[branch_var], root_vstatus_, original_lp_); - heap_.push(search_tree.root.get_down_child()); - heap_.push(search_tree.root.get_up_child()); - settings_.log.printf("Exploring the B&B tree using %d threads (%d are diving threads)\n", settings_.num_threads, settings_.num_threads - settings_.num_bfs_threads); @@ -1116,20 +1113,20 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut csc_matrix_t Arow(1, 1, 1); leaf_problem.A.transpose(Arow); - // #pragma omp master - // { - // i_t max_depth = std::ceil(std::log2(settings_.num_threads)); +#pragma omp master + { + i_t max_depth = std::ceil(std::log2(settings_.num_threads)); - // #pragma omp task - // exploration_ramp_up( - // &search_tree, search_tree.root.get_down_child(), leaf_problem, Arow, max_depth); +#pragma omp task + exploration_ramp_up( + &search_tree, search_tree.root.get_down_child(), leaf_problem, Arow, max_depth); - // #pragma omp task - // exploration_ramp_up( - // &search_tree, search_tree.root.get_up_child(), leaf_problem, Arow, max_depth); - // } +#pragma omp task + exploration_ramp_up( + &search_tree, search_tree.root.get_up_child(), leaf_problem, Arow, max_depth); + } - // #pragma omp barrier +#pragma omp barrier if (omp_get_thread_num() < settings_.num_bfs_threads) { best_first_thread(search_tree, leaf_problem, Arow); From 0e14fe2fb2f2fc2f31a9849e6859110b3326c98b Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 23 Sep 2025 17:45:23 +0200 Subject: [PATCH 08/33] fixed initialization order --- .../dual_simplex/simplex_solver_settings.hpp | 23 +++++++++---------- 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 3fad1d8ad1..e861c55fd1 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -97,15 +97,17 @@ struct simplex_solver_settings_t { bool use_bound_flip_ratio; // true if using the bound flip ratio test bool scale_columns; // true to scale the columns of A bool relaxation; // true to only solve the LP relaxation of a MIP - bool - use_left_looking_lu; // true to use left looking LU factorization, false to use right looking - bool eliminate_singletons; // true to eliminate singletons from the basis - bool print_presolve_stats; // true to print presolve stats - i_t refactor_frequency; // number of basis updates before refactorization - i_t iteration_log_frequency; // number of iterations between log updates - i_t first_iteration_log; // number of iterations to log at beginning of solve - i_t num_threads; // number of threads to use - i_t random_seed; // random seed + bool use_left_looking_lu; // true to use left looking LU factorization, + // false to use right looking + bool eliminate_singletons; // true to eliminate singletons from the basis + bool print_presolve_stats; // true to print presolve stats + i_t refactor_frequency; // number of basis updates before refactorization + i_t iteration_log_frequency; // number of iterations between log updates + i_t first_iteration_log; // number of iterations to log at beginning of solve + i_t num_threads; // number of threads to use + i_t num_bfs_threads; // number of threads dedicated to the best-first search + i_t num_diving_threads; // number of threads dedicated to diving + i_t random_seed; // random seed i_t inside_mip; // 0 if outside MIP, 1 if inside MIP at root node, 2 if inside MIP at leaf node std::function&, f_t)> solution_callback; std::function heuristic_preemption_callback; @@ -113,9 +115,6 @@ struct simplex_solver_settings_t { mutable logger_t log; std::atomic* concurrent_halt; // if nullptr ignored, if !nullptr, 0 if solver should // continue, 1 if solver should halt - - i_t num_bfs_threads; - i_t num_diving_threads; }; } // namespace cuopt::linear_programming::dual_simplex From c41af70db9e4f2967613143a2b9d8db0922a4ba0 Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 23 Sep 2025 21:38:07 +0200 Subject: [PATCH 09/33] fixed incorrect lower bounds Signed-off-by: nicolas --- cpp/src/dual_simplex/branch_and_bound.cpp | 18 +++++++++++++++++- cpp/src/dual_simplex/branch_and_bound.hpp | 2 ++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index c759d4a184..2d5badf9fe 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -15,6 +15,8 @@ * limitations under the License. */ +#include +#include #include #include #include @@ -222,10 +224,18 @@ f_t branch_and_bound_t::get_upper_bound() template f_t branch_and_bound_t::get_lower_bound() { - f_t lower_bound = root_objective_; + f_t lower_bound = inf; mutex_heap_.lock(); if (heap_.size() > 0) { lower_bound = heap_.top()->lower_bound; } mutex_heap_.unlock(); + + for (i_t i = 0; i < lower_bounds_.size(); ++i) { + f_t lb; +#pragma omp atomic read + lb = lower_bounds_[i]; + lower_bound = std::min(lb, lower_bound); + } + return lower_bound; } @@ -759,6 +769,7 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear f_t last_log = 0; i_t nodes_explored = 0; i_t nodes_unexplored = 0; + i_t tid = omp_get_thread_num(); while (stack.size() > 0 && get_status() == mip_status_t::RUNNING) { repair_heuristic_solutions(); @@ -770,6 +781,9 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear upper_bound = get_upper_bound(); gap = upper_bound - lower_bound; +#pragma omp atomic write + lower_bounds_[tid] = lower_bound; + #pragma omp atomic capture nodes_explored = stats_.nodes_explored++; @@ -1028,6 +1042,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut set_uninitialized_steepest_edge_norms(edge_norms_); root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + lower_bounds_.assign(settings_.num_bfs_threads, root_objective_); + if (settings_.set_simplex_solution_callback != nullptr) { std::vector original_x; uncrush_primal_solution(original_problem_, original_lp_, root_relax_soln_.x, original_x); diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 61de2839c9..dbb79a6248 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -124,6 +124,8 @@ class branch_and_bound_t { std::vector new_slacks_; std::vector var_types_; + std::vector lower_bounds_; + // Mutex for upper bound std::mutex mutex_upper_; From d83148e76953af0b48014304ec53692f74900261 Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 24 Sep 2025 18:44:33 +0200 Subject: [PATCH 10/33] added a wrapper class for omp_lock and omp_atomic. --- cpp/src/dual_simplex/branch_and_bound.cpp | 131 ++++++--------------- cpp/src/dual_simplex/branch_and_bound.hpp | 56 +++++---- cpp/src/dual_simplex/mip_node.hpp | 43 ++++--- cpp/src/utilities/omp_helpers.hpp | 136 ++++++++++++++++++++++ 4 files changed, 219 insertions(+), 147 deletions(-) create mode 100644 cpp/src/utilities/omp_helpers.hpp diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 2d5badf9fe..fbe78f66c8 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -170,7 +170,6 @@ f_t relative_gap(f_t obj_value, f_t lower_bound) f_t user_mip_gap = obj_value == 0.0 ? (lower_bound == 0.0 ? 0.0 : std::numeric_limits::infinity()) : std::abs(obj_value - lower_bound) / std::abs(obj_value); - // Handle NaNs (i.e., NaN != NaN) if (std::isnan(user_mip_gap)) { return std::numeric_limits::infinity(); } return user_mip_gap; } @@ -230,26 +229,12 @@ f_t branch_and_bound_t::get_lower_bound() mutex_heap_.unlock(); for (i_t i = 0; i < lower_bounds_.size(); ++i) { - f_t lb; -#pragma omp atomic read - lb = lower_bounds_[i]; - lower_bound = std::min(lb, lower_bound); + lower_bound = std::min(lower_bounds_[i].load(), lower_bound); } return lower_bound; } -template -mip_status_t branch_and_bound_t::get_status() -{ - mip_status_t status; - -#pragma omp atomic read - status = status_; - - return status; -} - template i_t branch_and_bound_t::get_heap_size() { @@ -259,15 +244,6 @@ i_t branch_and_bound_t::get_heap_size() return size; } -template -i_t branch_and_bound_t::get_active_subtrees() -{ - i_t active_subtrees; -#pragma omp atomic read - active_subtrees = active_subtrees_; - return active_subtrees; -} - template void branch_and_bound_t::set_new_solution(const std::vector& solution) { @@ -306,7 +282,7 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu mutex_upper_.unlock(); if (is_feasible) { - if (get_status() == mip_status_t::RUNNING) { + if (status_ == mip_status_t::RUNNING) { f_t user_obj = compute_user_objective(original_lp_, obj); f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); std::string gap = user_mip_gap(user_obj, user_lower); @@ -450,15 +426,9 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, i_t leaf_depth, char symbol) { - bool send_solution = false; - i_t nodes_explored; - i_t nodes_unexplored; - -#pragma omp atomic read - nodes_explored = stats_.nodes_explored; - -#pragma omp atomic read - nodes_unexplored = stats_.nodes_unexplored; + bool send_solution = false; + i_t nodes_explored = stats_.nodes_explored; + i_t nodes_unexplored = stats_.nodes_unexplored; mutex_upper_.lock(); if (leaf_objective < upper_bound_) { @@ -554,12 +524,8 @@ dual::status_t branch_and_bound_t::node_dual_simplex( } } -#pragma omp atomic update stats_.total_lp_solve_time += toc(lp_start_time); - -#pragma omp atomic update stats_.total_lp_iters += node_iter; - return lp_status; } @@ -676,7 +642,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* csc_matrix_t& Arow, i_t max_depth) { - if (get_status() != mip_status_t::RUNNING) { return; } + if (status_ != mip_status_t::RUNNING) { return; } repair_heuristic_solutions(); f_t lower_bound = node->lower_bound; @@ -685,11 +651,8 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* i_t nodes_explored = 0; i_t nodes_unexplored = 0; -#pragma omp atomic capture - nodes_explored = stats_.nodes_explored++; - -#pragma omp atomic capture - nodes_unexplored = stats_.nodes_unexplored--; + nodes_explored = (stats_.nodes_explored++); + nodes_unexplored = (stats_.nodes_unexplored--); if (upper_bound < node->lower_bound || relative_gap(upper_bound, lower_bound) < settings_.relative_mip_gap_tol) { @@ -700,28 +663,24 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* f_t now = toc(stats_.start_time); -#pragma omp single nowait - { - if (nodes_explored % 1000 == 0 || gap < 10 * settings_.absolute_mip_gap_tol || - nodes_explored < 1000) { - f_t obj = compute_user_objective(original_lp_, upper_bound); - f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); - std::string gap_user = user_mip_gap(obj, user_lower); - - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - nodes_explored, - nodes_unexplored, - obj, - user_lower, - node->depth, - nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, - gap_user.c_str(), - now); - } + if (nodes_explored % 1000 == 0 || gap < 10 * settings_.absolute_mip_gap_tol || + nodes_explored < 1000) { + f_t obj = compute_user_objective(original_lp_, upper_bound); + f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); + std::string gap_user = user_mip_gap(obj, user_lower); + + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + nodes_explored, + nodes_unexplored, + obj, + user_lower, + node->depth, + nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, + gap_user.c_str(), + now); } if (toc(stats_.start_time) > settings_.time_limit) { -#pragma omp atomic write status_ = mip_status_t::TIME_LIMIT; return; } @@ -729,12 +688,10 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* solve_node_lp(*search_tree, node, leaf_problem, Arow, upper_bound, settings_.log, 'B'); if (node_status == node_status_t::TIME_LIMIT) { -#pragma omp atomic write status_ = mip_status_t::TIME_LIMIT; return; } else if (node_status == node_status_t::HAS_CHILDREN) { -#pragma omp atomic update stats_.nodes_unexplored += 2; if (node->depth < max_depth) { @@ -771,7 +728,7 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear i_t nodes_unexplored = 0; i_t tid = omp_get_thread_num(); - while (stack.size() > 0 && get_status() == mip_status_t::RUNNING) { + while (stack.size() > 0 && status_ == mip_status_t::RUNNING) { repair_heuristic_solutions(); mip_node_t* node_ptr = stack.front(); @@ -781,14 +738,9 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear upper_bound = get_upper_bound(); gap = upper_bound - lower_bound; -#pragma omp atomic write lower_bounds_[tid] = lower_bound; - -#pragma omp atomic capture - nodes_explored = stats_.nodes_explored++; - -#pragma omp atomic capture - nodes_unexplored = stats_.nodes_unexplored--; + nodes_explored = stats_.nodes_explored++; + nodes_unexplored = stats_.nodes_unexplored--; if (upper_bound < node_ptr->lower_bound || relative_gap(upper_bound, lower_bound) < settings_.relative_mip_gap_tol) { @@ -824,7 +776,6 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear } if (toc(stats_.start_time) > settings_.time_limit) { -#pragma omp atomic write status_ = mip_status_t::TIME_LIMIT; break; } @@ -835,7 +786,6 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear solve_node_lp(search_tree, node_ptr, leaf_problem, Arow, upper_bound, settings_.log, 'B'); if (node_status == node_status_t::TIME_LIMIT) { -#pragma omp atomic write status_ = mip_status_t::TIME_LIMIT; return; @@ -855,7 +805,6 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear mutex_heap_.unlock(); } -#pragma omp atomic update stats_.nodes_unexplored += 2; auto [first, second] = child_selection(node_ptr); @@ -882,8 +831,6 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se if (heap_.size() > 0) { node_ptr = heap_.top(); heap_.pop(); - -#pragma omp atomic update active_subtrees_++; } mutex_heap_.unlock(); @@ -894,16 +841,12 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se // current upper bound search_tree.graphviz_node(node_ptr, "cutoff", node_ptr->lower_bound); search_tree.update_tree(node_ptr, node_status_t::FATHOMED); - -#pragma omp atomic update active_subtrees_--; continue; } // Best-first search with plunging explore_subtree(search_tree, node_ptr, leaf_problem, Arow); - -#pragma omp atomic update active_subtrees_--; } @@ -911,15 +854,12 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se upper_bound = get_upper_bound(); gap = upper_bound - lower_bound; - } while (get_status() == mip_status_t::RUNNING && gap > settings_.absolute_mip_gap_tol && + } while (status_ == mip_status_t::RUNNING && gap > settings_.absolute_mip_gap_tol && relative_gap(upper_bound, lower_bound) > settings_.relative_mip_gap_tol && - (get_active_subtrees() > 0 || get_heap_size() > 0)); + (active_subtrees_ > 0 || get_heap_size() > 0)); // Check if the solver exited naturally, or was due to a timeout or numerical error. - if (status_ == mip_status_t::RUNNING) { -#pragma omp atomic write - status_ = mip_status_t::FINISHED; - } + if (status_ == mip_status_t::RUNNING) { status_ = mip_status_t::COMPLETED; } } template @@ -943,7 +883,7 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr std::deque*> stack; stack.push_front(&subtree.root); - while (stack.size() > 0 && get_status() == mip_status_t::RUNNING) { + while (stack.size() > 0 && status_ == mip_status_t::RUNNING) { mip_node_t* node_ptr = stack.front(); stack.pop_front(); f_t upper_bound = get_upper_bound(); @@ -975,8 +915,7 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr } } - } while (get_status() == mip_status_t::RUNNING && - (get_active_subtrees() > 0 || get_heap_size() > 0)); + } while (status_ == mip_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)); } template @@ -1116,10 +1055,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.num_threads, settings_.num_threads - settings_.num_bfs_threads); settings_.log.printf( - "| Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap | " - " Time \n"); + "| Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap " + "| " + " Time \n"); -#pragma omp atomic write status_ = mip_status_t::RUNNING; #pragma omp parallel num_threads(settings_.num_threads) @@ -1160,7 +1099,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t gap = upper_bound - lower_bound; f_t obj = compute_user_objective(original_lp_, upper_bound); f_t user_lower = compute_user_objective(original_lp_, lower_bound); - f_t gap_rel = relative_gap(get_upper_bound(), lower_bound); + f_t gap_rel = relative_gap(upper_bound, lower_bound); settings_.log.printf( "Explored %d nodes in %.2fs.\n", stats_.nodes_explored, toc(stats_.start_time)); diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index dbb79a6248..9773cd01b6 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -25,24 +25,24 @@ #include #include #include +#include #include -#include #include #include namespace cuopt::linear_programming::dual_simplex { enum class mip_status_t { - OPTIMAL = 0, - UNBOUNDED = 1, - INFEASIBLE = 2, - TIME_LIMIT = 3, - NODE_LIMIT = 4, - NUMERICAL = 5, - UNSET = 6, - RUNNING = 7, - FINISHED = 8 + OPTIMAL = 0, // The optimal integer solution was found + UNBOUNDED = 1, // The problem is unbounded + INFEASIBLE = 2, // The problem is infeasible + TIME_LIMIT = 3, // The solver reached a time limit + NODE_LIMIT = 4, // The maximum number of nodes was reached (not implemented) + NUMERICAL = 5, // The solver encountered a numerical error + RUNNING = 6, // The solver is currently exploring the tree + COMPLETED = 7, // The solver finished exploring the tree + UNSET = 8, // The status is not set }; template @@ -55,7 +55,6 @@ template class dive_queue_t { private: std::vector> buffer; - std::mutex mutex; static constexpr i_t max_size_ = 8192; public: @@ -63,7 +62,7 @@ class dive_queue_t { void push(mip_node_t&& node) { - buffer.push_back(std::forward&&>(node)); + buffer.push_back(std::move(node)); std::push_heap(buffer.begin(), buffer.end(), node_compare_t()); if (buffer.size() > max_size()) { buffer.pop_back(); } } @@ -105,9 +104,7 @@ class branch_and_bound_t { f_t get_upper_bound(); f_t get_lower_bound(); - mip_status_t get_status(); i_t get_heap_size(); - i_t get_active_subtrees(); // The main entry routine. Returns the solver status and populates solution with the incumbent. mip_status_t solve(mip_solution_t& solution); @@ -124,10 +121,11 @@ class branch_and_bound_t { std::vector new_slacks_; std::vector var_types_; - std::vector lower_bounds_; + // Local lower bounds for each thread + std::vector> lower_bounds_; // Mutex for upper bound - std::mutex mutex_upper_; + omp_mutex_t mutex_upper_; // Global variable for upper bound f_t upper_bound_; @@ -137,15 +135,15 @@ class branch_and_bound_t { // Structure with the general info of the solver. struct stats_t { - f_t start_time = 0.0; - f_t total_lp_solve_time = 0.0; - i_t nodes_explored = 0; - i_t nodes_unexplored = 0; - f_t total_lp_iters = 0; + f_t start_time = 0.0; + omp_atomic_t total_lp_solve_time = 0.0; + omp_atomic_t nodes_explored = 0; + omp_atomic_t nodes_unexplored = 0; + omp_atomic_t total_lp_iters = 0; } stats_; // Mutex for repair - std::mutex mutex_repair_; + omp_mutex_t mutex_repair_; std::vector> repair_queue_; // Variables for the root node in the search tree. @@ -156,22 +154,22 @@ class branch_and_bound_t { // Pseudocosts pseudo_costs_t pc_; - std::mutex mutex_pc_; + omp_mutex_t mutex_pc_; // Heap storing the nodes to be explored. - std::mutex mutex_heap_; + omp_mutex_t mutex_heap_; mip_node_heap_t*> heap_; - // Global status - mip_status_t status_; - // Count the number of subtrees that are currently being explored. - i_t active_subtrees_; + omp_atomic_t active_subtrees_; // Queue for storing the promising node for performing dives. - std::mutex mutex_dive_queue_; + omp_mutex_t mutex_dive_queue_; dive_queue_t dive_queue_; + // Global status + omp_atomic_t status_; + // Update the incumbent solution with the new feasible solution. // found during branch and bound. void add_feasible_solution(f_t leaf_objective, diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index ee76223f4d..12f01d3ed1 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -19,11 +19,11 @@ #include #include +#include #include #include #include -#include #include namespace cuopt::linear_programming::dual_simplex { @@ -202,6 +202,10 @@ class mip_node_t { } } + // This method creates a copy of the current node + // with its parent set to `nullptr`, `node_id = 0` + // and `depth = 0` such that it is the root + // of a separated tree. mip_node_t detach_copy() const { mip_node_t copy(lower_bound, vstatus); @@ -242,13 +246,13 @@ void remove_fathomed_nodes(std::vector*>& stack) template class node_compare_t { public: - bool operator()(mip_node_t& a, mip_node_t& b) const + bool operator()(const mip_node_t& a, const mip_node_t& b) const { return a.lower_bound > b.lower_bound; // True if a comes before b, elements that come before are output last } - bool operator()(mip_node_t* a, mip_node_t* b) const + bool operator()(const mip_node_t* a, const mip_node_t* b) const { return a->lower_bound > b->lower_bound; // True if a comes before b, elements that come before are output last @@ -264,7 +268,7 @@ class search_tree_t { } search_tree_t(mip_node_t&& node, logger_t& log) - : root(std::forward&&>(node)), num_nodes(0), log(log) + : root(std::move(node)), num_nodes(0), log(log) { } @@ -272,25 +276,20 @@ class search_tree_t { void update_tree(mip_node_t* node_ptr, node_status_t status) { - std::lock_guard lock(mutex); + mutex.lock(); std::vector*> stack; node_ptr->set_status(status, stack); remove_fathomed_nodes(stack); + mutex.unlock(); } void branch(mip_node_t* parent_node, - i_t branch_var, - f_t branch_var_val, + const i_t branch_var, + const f_t branch_var_val, const std::vector& parent_vstatus, const lp_problem_t& original_lp) { - i_t id; - -#pragma omp atomic capture - { - id = num_nodes; - num_nodes += 2; - } + i_t id = num_nodes.fetch_add(2); // down child auto down_child = std::make_unique>( @@ -309,18 +308,18 @@ class search_tree_t { std::move(up_child)); // child pointers moved into the tree } - void graphviz_node(mip_node_t* node_ptr, std::string label, f_t val) + void graphviz_node(const mip_node_t* node_ptr, const std::string label, const f_t val) { if (write_graphviz) { log.printf("Node%d [label=\"%s %.16e\"]\n", node_ptr->node_id, label.c_str(), val); } } - void graphviz_edge(mip_node_t* origin_ptr, - mip_node_t* dest_ptr, - i_t branch_var, - i_t branch_dir, - f_t bound) + void graphviz_edge(const mip_node_t* origin_ptr, + const mip_node_t* dest_ptr, + const i_t branch_var, + const i_t branch_dir, + const f_t bound) { if (write_graphviz) { log.printf("Node%d -> Node%d [label=\"x%d %s %e\"]\n", @@ -333,8 +332,8 @@ class search_tree_t { } mip_node_t root; - std::mutex mutex; - i_t num_nodes; + omp_mutex_t mutex; + omp_atomic_t num_nodes; logger_t log; static constexpr int write_graphviz = false; diff --git a/cpp/src/utilities/omp_helpers.hpp b/cpp/src/utilities/omp_helpers.hpp new file mode 100644 index 0000000000..5ba50e78ff --- /dev/null +++ b/cpp/src/utilities/omp_helpers.hpp @@ -0,0 +1,136 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +#ifdef _OPENMP + +#include +#include + +namespace cuopt { + +// Wrapper of omp_lock_t. Optionally, you can provide a hint as defined in +// https://www.openmp.org/spec-html/5.1/openmpse39.html#x224-2570003.9 +class omp_mutex_t { + public: + omp_mutex_t() { omp_init_lock(&mutex); } + virtual ~omp_mutex_t() { omp_destroy_lock(&mutex); } + void lock() { omp_set_lock(&mutex); } + void unlock() { omp_unset_lock(&mutex); } + bool try_lock() { return omp_test_lock(&mutex); } + + private: + omp_lock_t mutex; +}; + +// Wrapper for omp atomic operations. See +// https://www.openmp.org/spec-html/5.1/openmpsu105.html. +template +class omp_atomic_t { + public: + omp_atomic_t() = default; + omp_atomic_t(T val) : val(val) {} + + T operator=(T new_val) + { + store(new_val); + return new_val; + } + + operator T() { return load(); } + T operator+=(T inc) { return fetch_add(inc) + inc; } + T operator-=(T inc) { return fetch_sub(inc) - inc; } + + // In theory, this should be enabled only for integers, + // but it works for any numerical types. + T operator++() { return fetch_add(T(1)) + 1; } + T operator++(int) { return fetch_add(T(1)); } + T operator--() { return fetch_add(T(1)) - 1; } + T operator--(int) { return fetch_add(T(1)); } + + T load() + { + T res; +#pragma omp atomic read + res = val; + return res; + } + + void store(T new_val) + { +#pragma omp atomic write + val = new_val; + } + + T exchange(T other) + { + T old; +#pragma omp atomic capture + { + old = val; + val = other; + } + return old; + } + + T fetch_add(T inc) + { + T old; +#pragma omp atomic capture + { + old = val; + val += inc; + } + return old; + } + + T fetch_sub(T inc) { return fetch_add(-inc); } + +// OpenMP only supports atomic CAS operations in v5.1 +#if _OPENMP == 202011 + + T fetch_min(T other) + { + T old; +#pragma omp atomic compare capture + { + old = val; + val = other < val ? other : val; + } + return old; + } + + T fetch_max(T other) + { + T old; +#pragma omp atomic compare capture + { + old = val; + val = other > val ? other : val; + } + return old; + } +#endif + + private: + T val; +}; + +#endif + +} // namespace cuopt From bbdee5e9823352c32a1396425e8553d0a236247d Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 25 Sep 2025 17:48:05 +0200 Subject: [PATCH 11/33] fixed incorrect convergence check (#417) --- cpp/src/dual_simplex/branch_and_bound.cpp | 60 ++++++++++++++--------- cpp/src/dual_simplex/branch_and_bound.hpp | 3 ++ cpp/src/utilities/omp_helpers.hpp | 4 +- 3 files changed, 41 insertions(+), 26 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index fbe78f66c8..581679c817 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -174,6 +174,26 @@ f_t relative_gap(f_t obj_value, f_t lower_bound) return user_mip_gap; } +template +f_t user_gap(const lp_problem_t& lp, f_t obj_value, f_t lower_bound) +{ + f_t user_obj = compute_user_objective(lp, obj_value); + f_t user_lower_bound = compute_user_objective(lp, lower_bound); + return user_obj - user_lower_bound; +} + +template +f_t user_relative_gap(const lp_problem_t& lp, f_t obj_value, f_t lower_bound) +{ + f_t user_obj = compute_user_objective(lp, obj_value); + f_t user_lower_bound = compute_user_objective(lp, lower_bound); + f_t user_mip_gap = user_obj == 0.0 + ? (user_lower_bound == 0.0 ? 0.0 : std::numeric_limits::infinity()) + : std::abs(user_obj - user_lower_bound) / std::abs(user_obj); + if (std::isnan(user_mip_gap)) { return std::numeric_limits::infinity(); } + return user_mip_gap; +} + template std::string user_mip_gap(f_t obj_value, f_t lower_bound) { @@ -244,6 +264,14 @@ i_t branch_and_bound_t::get_heap_size() return size; } +template +bool branch_and_bound_t::check_gap_convergence(f_t lower_bound, f_t upper_bound) +{ + f_t gap_rel = user_relative_gap(original_lp_, upper_bound, lower_bound); + f_t gap_abs = user_gap(original_lp_, upper_bound, lower_bound); + return gap_rel < settings_.relative_mip_gap_tol || gap_abs < settings_.absolute_mip_gap_tol; +} + template void branch_and_bound_t::set_new_solution(const std::vector& solution) { @@ -542,7 +570,6 @@ node_status_t branch_and_bound_t::solve_node_lp(search_tree_t& leaf_vstatus = node_ptr->vstatus; lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); @@ -601,8 +628,7 @@ node_status_t branch_and_bound_t::solve_node_lp(search_tree_t rel_fathom_tol) { + } else if (leaf_objective <= upper_bound + abs_fathom_tol) { // Choose fractional variable to branch on mutex_pc_.lock(); const i_t branch_var = pc_.variable_selection( @@ -654,8 +680,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* nodes_explored = (stats_.nodes_explored++); nodes_unexplored = (stats_.nodes_unexplored--); - if (upper_bound < node->lower_bound || - relative_gap(upper_bound, lower_bound) < settings_.relative_mip_gap_tol) { + if (check_gap_convergence(lower_bound, upper_bound)) { search_tree->graphviz_node(node, "cutoff", node->lower_bound); search_tree->update_tree(node, node_status_t::FATHOMED); return; @@ -742,8 +767,7 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear nodes_explored = stats_.nodes_explored++; nodes_unexplored = stats_.nodes_unexplored--; - if (upper_bound < node_ptr->lower_bound || - relative_gap(upper_bound, lower_bound) < settings_.relative_mip_gap_tol) { + if (check_gap_convergence(lower_bound, upper_bound)) { search_tree.graphviz_node(node_ptr, "cutoff", node_ptr->lower_bound); search_tree.update_tree(node_ptr, node_status_t::FATHOMED); continue; @@ -819,10 +843,6 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se lp_problem_t& leaf_problem, csc_matrix_t& Arow) { - f_t lower_bound = -inf; - f_t upper_bound = inf; - f_t gap = inf; - do { mip_node_t* node_ptr = nullptr; @@ -850,12 +870,8 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se active_subtrees_--; } - lower_bound = get_lower_bound(); - upper_bound = get_upper_bound(); - gap = upper_bound - lower_bound; - - } while (status_ == mip_status_t::RUNNING && gap > settings_.absolute_mip_gap_tol && - relative_gap(upper_bound, lower_bound) > settings_.relative_mip_gap_tol && + } while (status_ == mip_status_t::RUNNING && + !check_gap_convergence(get_lower_bound(), get_upper_bound()) && (active_subtrees_ > 0 || get_heap_size() > 0)); // Check if the solver exited naturally, or was due to a timeout or numerical error. @@ -888,10 +904,7 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr stack.pop_front(); f_t upper_bound = get_upper_bound(); - if (upper_bound < node_ptr->lower_bound || - relative_gap(upper_bound, node_ptr->lower_bound) < settings_.relative_mip_gap_tol) { - continue; - } + if (check_gap_convergence(node_ptr->lower_bound, upper_bound)) { continue; } node_status_t node_status = solve_node_lp(subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); @@ -1056,8 +1069,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.num_threads - settings_.num_bfs_threads); settings_.log.printf( "| Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap " - "| " - " Time \n"); + "| Time \n"); status_ = mip_status_t::RUNNING; @@ -1099,7 +1111,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t gap = upper_bound - lower_bound; f_t obj = compute_user_objective(original_lp_, upper_bound); f_t user_lower = compute_user_objective(original_lp_, lower_bound); - f_t gap_rel = relative_gap(upper_bound, lower_bound); + f_t gap_rel = user_relative_gap(original_lp_, upper_bound, lower_bound); settings_.log.printf( "Explored %d nodes in %.2fs.\n", stats_.nodes_explored, toc(stats_.start_time)); diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 9773cd01b6..eb26833c81 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -106,6 +106,9 @@ class branch_and_bound_t { f_t get_lower_bound(); i_t get_heap_size(); + // Check if the gap and relative gap are below the tolerance. + bool check_gap_convergence(f_t lower_bound, f_t upper_bound); + // The main entry routine. Returns the solver status and populates solution with the incumbent. mip_status_t solve(mip_solution_t& solution); diff --git a/cpp/src/utilities/omp_helpers.hpp b/cpp/src/utilities/omp_helpers.hpp index 5ba50e78ff..f5f45efb2d 100644 --- a/cpp/src/utilities/omp_helpers.hpp +++ b/cpp/src/utilities/omp_helpers.hpp @@ -60,8 +60,8 @@ class omp_atomic_t { // but it works for any numerical types. T operator++() { return fetch_add(T(1)) + 1; } T operator++(int) { return fetch_add(T(1)); } - T operator--() { return fetch_add(T(1)) - 1; } - T operator--(int) { return fetch_add(T(1)); } + T operator--() { return fetch_sub(T(1)) - 1; } + T operator--(int) { return fetch_sub(T(1)); } T load() { From 7f59c73fd8a5882ca14e8f131772aca5945e867e Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 25 Sep 2025 18:21:29 +0200 Subject: [PATCH 12/33] set the default number of threads --- benchmarks/linear_programming/cuopt/run_mip.cpp | 3 +++ cpp/src/dual_simplex/phase2.cpp | 4 +--- cpp/src/dual_simplex/simplex_solver_settings.hpp | 8 ++++---- cpp/src/mip/diversity/recombiners/sub_mip.cuh | 8 +++++--- cpp/src/mip/solver.cu | 10 +++++++++- 5 files changed, 22 insertions(+), 11 deletions(-) diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index 64b1bb4bb1..a4f52cb4ed 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -36,6 +36,7 @@ #include #include +#include #include #include #include @@ -382,6 +383,8 @@ int main(int argc, char* argv[]) double memory_limit = program.get("--memory-limit"); bool track_allocations = program.get("--track-allocations")[0] == 't'; + if (num_cpu_threads < 0) { num_cpu_threads = omp_get_max_threads() / n_gpus; } + if (program.is_used("--out-dir")) { out_dir = program.get("--out-dir"); result_file = out_dir + "/final_result.csv"; diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 370badf332..01d995ba6b 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -29,9 +29,7 @@ #include #include #include -#include -#include -#include +#include namespace cuopt::linear_programming::dual_simplex { diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index e861c55fd1..9ff96ebbd0 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -20,11 +20,11 @@ #include #include +#include #include #include #include #include -#include namespace cuopt::linear_programming::dual_simplex { @@ -59,9 +59,9 @@ struct simplex_solver_settings_t { refactor_frequency(100), iteration_log_frequency(1000), first_iteration_log(2), - num_threads(std::thread::hardware_concurrency() - 1), - num_bfs_threads(std::min(num_threads, 4)), - num_diving_threads(num_threads - num_bfs_threads), + num_threads(omp_get_max_threads() - 1), + num_bfs_threads(std::min(num_threads / 4, 1)), + num_diving_threads(std::min(num_threads - num_bfs_threads, 1)), random_seed(0), inside_mip(0), solution_callback(nullptr), diff --git a/cpp/src/mip/diversity/recombiners/sub_mip.cuh b/cpp/src/mip/diversity/recombiners/sub_mip.cuh index de6d917e7c..fa5968c6f5 100644 --- a/cpp/src/mip/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip/diversity/recombiners/sub_mip.cuh @@ -111,9 +111,11 @@ class sub_mip_recombiner_t : public recombiner_t { branch_and_bound_settings.print_presolve_stats = false; branch_and_bound_settings.absolute_mip_gap_tol = context.settings.tolerances.absolute_mip_gap; branch_and_bound_settings.relative_mip_gap_tol = context.settings.tolerances.relative_mip_gap; - branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance; - branch_and_bound_settings.num_threads = 1; - branch_and_bound_settings.solution_callback = [this](std::vector& solution, + branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance; + branch_and_bound_settings.num_threads = 2; + branch_and_bound_settings.num_bfs_threads = 1; + branch_and_bound_settings.num_diving_threads = 1; + branch_and_bound_settings.solution_callback = [this](std::vector& solution, f_t objective) { this->solution_callback(solution, objective); }; diff --git a/cpp/src/mip/solver.cu b/cpp/src/mip/solver.cu index ea32724191..95eaacc184 100644 --- a/cpp/src/mip/solver.cu +++ b/cpp/src/mip/solver.cu @@ -172,10 +172,18 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.relative_mip_gap_tol = context.settings.tolerances.relative_mip_gap; branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance; - if (context.settings.num_cpu_threads != -1) { + if (context.settings.num_cpu_threads < 0) { + branch_and_bound_settings.num_threads = omp_get_max_threads() - 1; + } else { branch_and_bound_settings.num_threads = std::max(1, context.settings.num_cpu_threads); } + i_t num_threads = branch_and_bound_settings.num_threads; + i_t num_bfs_threads = std::max(1, num_threads / 4); + i_t num_diving_threads = std::max(1, num_threads - num_bfs_threads); + branch_and_bound_settings.num_bfs_threads = num_bfs_threads; + branch_and_bound_settings.num_diving_threads = num_diving_threads; + // Set the branch and bound -> primal heuristics callback branch_and_bound_settings.solution_callback = std::bind(&branch_and_bound_solution_helper_t::solution_callback, From 424331b038e77dc642b01bc2a83e843f1e37f9db Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 25 Sep 2025 18:45:57 +0200 Subject: [PATCH 13/33] fixed log spacing --- cpp/src/dual_simplex/branch_and_bound.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 581679c817..17619dd4a6 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -37,6 +37,7 @@ #include #include #include +#include "cuopt/logger_macros.hpp" namespace cuopt::linear_programming::dual_simplex { @@ -316,7 +317,7 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu std::string gap = user_mip_gap(user_obj, user_lower); settings_.log.printf( - "H %+13.6e %+10.6e %s %9.2f\n", + "H %+13.6e %+10.6e %s %9.2f\n", user_obj, user_lower, gap.c_str(), @@ -465,7 +466,7 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, f_t lower_bound = get_lower_bound(); f_t obj = compute_user_objective(original_lp_, upper_bound_); f_t lower = compute_user_objective(original_lp_, lower_bound); - settings_.log.printf("%c%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + settings_.log.printf("%c%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", symbol, nodes_explored, nodes_unexplored, @@ -694,7 +695,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); std::string gap_user = user_mip_gap(obj, user_lower); - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", nodes_explored, nodes_unexplored, obj, @@ -785,8 +786,7 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear f_t obj = compute_user_objective(original_lp_, upper_bound); f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); std::string gap_user = user_mip_gap(obj, user_lower); - - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", nodes_explored, nodes_unexplored, obj, @@ -1064,12 +1064,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut search_tree.branch( &search_tree.root, branch_var, root_relax_soln_.x[branch_var], root_vstatus_, original_lp_); - settings_.log.printf("Exploring the B&B tree using %d threads (%d are diving threads)\n", + settings_.log.printf("Exploring the B&B tree using %d threads (%d diving threads)\n", settings_.num_threads, - settings_.num_threads - settings_.num_bfs_threads); + settings_.num_diving_threads); settings_.log.printf( - "| Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap " - "| Time \n"); + " | Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap " + "| Time |\n"); status_ = mip_status_t::RUNNING; From a8f3e933a2394f33a61fa2f18d72366e0abf99d0 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 26 Sep 2025 14:07:50 +0200 Subject: [PATCH 14/33] fixed incorrect gap termination when solving a maximation problem --- cpp/src/dual_simplex/branch_and_bound.cpp | 78 +++++++++-------------- cpp/src/dual_simplex/branch_and_bound.hpp | 2 +- 2 files changed, 31 insertions(+), 49 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 17619dd4a6..f4463da09f 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -37,7 +37,6 @@ #include #include #include -#include "cuopt/logger_macros.hpp" namespace cuopt::linear_programming::dual_simplex { @@ -175,14 +174,6 @@ f_t relative_gap(f_t obj_value, f_t lower_bound) return user_mip_gap; } -template -f_t user_gap(const lp_problem_t& lp, f_t obj_value, f_t lower_bound) -{ - f_t user_obj = compute_user_objective(lp, obj_value); - f_t user_lower_bound = compute_user_objective(lp, lower_bound); - return user_obj - user_lower_bound; -} - template f_t user_relative_gap(const lp_problem_t& lp, f_t obj_value, f_t lower_bound) { @@ -266,10 +257,10 @@ i_t branch_and_bound_t::get_heap_size() } template -bool branch_and_bound_t::check_gap_convergence(f_t lower_bound, f_t upper_bound) +bool branch_and_bound_t::check_gap(f_t lower_bound, f_t upper_bound) { f_t gap_rel = user_relative_gap(original_lp_, upper_bound, lower_bound); - f_t gap_abs = user_gap(original_lp_, upper_bound, lower_bound); + f_t gap_abs = upper_bound - lower_bound; return gap_rel < settings_.relative_mip_gap_tol || gap_abs < settings_.absolute_mip_gap_tol; } @@ -674,14 +665,13 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* f_t lower_bound = node->lower_bound; f_t upper_bound = get_upper_bound(); - f_t gap = upper_bound - lower_bound; i_t nodes_explored = 0; i_t nodes_unexplored = 0; nodes_explored = (stats_.nodes_explored++); nodes_unexplored = (stats_.nodes_unexplored--); - if (check_gap_convergence(lower_bound, upper_bound)) { + if (check_gap(lower_bound, upper_bound)) { search_tree->graphviz_node(node, "cutoff", node->lower_bound); search_tree->update_tree(node, node_status_t::FATHOMED); return; @@ -689,8 +679,8 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* f_t now = toc(stats_.start_time); - if (nodes_explored % 1000 == 0 || gap < 10 * settings_.absolute_mip_gap_tol || - nodes_explored < 1000) { + if (nodes_explored % 1000 == 0 || + (upper_bound - lower_bound) < 10 * settings_.absolute_mip_gap_tol || nodes_explored < 1000) { f_t obj = compute_user_objective(original_lp_, upper_bound); f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); std::string gap_user = user_mip_gap(obj, user_lower); @@ -746,9 +736,6 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear stack.push_front(start_node); // It contains the lower bound of the parent for now - f_t lower_bound = start_node->lower_bound; - f_t upper_bound = get_upper_bound(); - f_t gap = upper_bound - lower_bound; f_t last_log = 0; i_t nodes_explored = 0; i_t nodes_unexplored = 0; @@ -760,43 +747,39 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear mip_node_t* node_ptr = stack.front(); stack.pop_front(); - lower_bound = node_ptr->lower_bound; - upper_bound = get_upper_bound(); - gap = upper_bound - lower_bound; + f_t lower_bound = node_ptr->lower_bound; + f_t upper_bound = get_upper_bound(); lower_bounds_[tid] = lower_bound; nodes_explored = stats_.nodes_explored++; nodes_unexplored = stats_.nodes_unexplored--; - if (check_gap_convergence(lower_bound, upper_bound)) { + if (check_gap(lower_bound, upper_bound)) { search_tree.graphviz_node(node_ptr, "cutoff", node_ptr->lower_bound); search_tree.update_tree(node_ptr, node_status_t::FATHOMED); continue; } - f_t now = toc(stats_.start_time); - -#pragma omp single nowait - { - f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); - if (((nodes_explored % 1000 == 0 || gap < 10 * settings_.absolute_mip_gap_tol || - nodes_explored < 1000) && - (time_since_log >= 1)) || - (time_since_log > 60) || now > settings_.time_limit) { - f_t obj = compute_user_objective(original_lp_, upper_bound); - f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); - std::string gap_user = user_mip_gap(obj, user_lower); - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - nodes_explored, - nodes_unexplored, - obj, - user_lower, - node_ptr->depth, - nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, - gap_user.c_str(), - now); - last_log = tic(); - } + f_t now = toc(stats_.start_time); + f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); + if (((nodes_explored % 1000 == 0 || + (upper_bound - lower_bound) < 10 * settings_.absolute_mip_gap_tol || + nodes_explored < 1000) && + (time_since_log >= 1)) || + (time_since_log > 60) || now > settings_.time_limit) { + f_t obj = compute_user_objective(original_lp_, upper_bound); + f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); + std::string gap_user = user_mip_gap(obj, user_lower); + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + nodes_explored, + nodes_unexplored, + obj, + user_lower, + node_ptr->depth, + nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, + gap_user.c_str(), + now); + last_log = tic(); } if (toc(stats_.start_time) > settings_.time_limit) { @@ -870,8 +853,7 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se active_subtrees_--; } - } while (status_ == mip_status_t::RUNNING && - !check_gap_convergence(get_lower_bound(), get_upper_bound()) && + } while (status_ == mip_status_t::RUNNING && !check_gap(get_lower_bound(), get_upper_bound()) && (active_subtrees_ > 0 || get_heap_size() > 0)); // Check if the solver exited naturally, or was due to a timeout or numerical error. @@ -904,7 +886,7 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr stack.pop_front(); f_t upper_bound = get_upper_bound(); - if (check_gap_convergence(node_ptr->lower_bound, upper_bound)) { continue; } + if (check_gap(node_ptr->lower_bound, upper_bound)) { continue; } node_status_t node_status = solve_node_lp(subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index eb26833c81..82c0d9477a 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -107,7 +107,7 @@ class branch_and_bound_t { i_t get_heap_size(); // Check if the gap and relative gap are below the tolerance. - bool check_gap_convergence(f_t lower_bound, f_t upper_bound); + bool check_gap(f_t lower_bound, f_t upper_bound); // The main entry routine. Returns the solver status and populates solution with the incumbent. mip_status_t solve(mip_solution_t& solution); From 19afd5c8c8c7d2e5c844f2c4e4995f4cc97f7214 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 26 Sep 2025 15:10:46 +0200 Subject: [PATCH 15/33] fixed missing final report after timeout (#320) --- cpp/src/dual_simplex/branch_and_bound.cpp | 109 ++++++++++++---------- cpp/src/dual_simplex/branch_and_bound.hpp | 3 + 2 files changed, 64 insertions(+), 48 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index f4463da09f..640a7d5627 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -440,6 +440,55 @@ void branch_and_bound_t::repair_heuristic_solutions() } } +template +void branch_and_bound_t::set_final_solution(mip_solution_t& solution, + f_t lower_bound) +{ + if (status_ == mip_status_t::TIME_LIMIT) { + settings_.log.printf("Time limit reached. Stopping the solver...\n"); + } + + f_t upper_bound = get_upper_bound(); + f_t gap = upper_bound - lower_bound; + f_t obj = compute_user_objective(original_lp_, upper_bound); + f_t user_lower = compute_user_objective(original_lp_, lower_bound); + f_t gap_rel = user_relative_gap(original_lp_, upper_bound, lower_bound); + + settings_.log.printf( + "Explored %d nodes in %.2fs.\n", stats_.nodes_explored, toc(stats_.start_time)); + settings_.log.printf("Absolute Gap %e Objective %.16e Lower Bound %.16e\n", gap, obj, user_lower); + + if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) { + status_ = mip_status_t::OPTIMAL; + if (gap > 0 && gap <= settings_.absolute_mip_gap_tol) { + settings_.log.printf("Optimal solution found within absolute MIP gap tolerance (%.1e)\n", + settings_.absolute_mip_gap_tol); + } else if (gap > 0 && gap_rel <= settings_.relative_mip_gap_tol) { + settings_.log.printf("Optimal solution found within relative MIP gap tolerance (%.1e)\n", + settings_.relative_mip_gap_tol); + } else { + settings_.log.printf("Optimal solution found.\n"); + } + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); + } + } + + if (stats_.nodes_explored > 0 && stats_.nodes_unexplored == 0 && upper_bound == inf) { + settings_.log.printf("Integer infeasible.\n"); + status_ = mip_status_t::INFEASIBLE; + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); + } + } + + uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); + solution.objective = incumbent_.objective; + solution.lower_bound = lower_bound; + solution.nodes_explored = stats_.nodes_explored; + solution.simplex_iterations = stats_.total_lp_iters; +} + template void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, const std::vector& leaf_solution, @@ -922,7 +971,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut stats_.total_lp_iters = 0; stats_.nodes_explored = 0; - stats_.nodes_unexplored = 2; + stats_.nodes_unexplored = 0; active_subtrees_ = 0; if (guess_.size() != 0) { @@ -949,6 +998,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut lp_settings.inside_mip = 1; lp_status_t root_status = solve_linear_program_advanced( original_lp_, stats_.start_time, lp_settings, root_relax_soln_, root_vstatus_, edge_norms_); + stats_.total_lp_iters = root_relax_soln_.iterations; stats_.total_lp_solve_time = toc(stats_.start_time); assert(root_vstatus_.size() == original_lp_.num_cols); if (root_status == lp_status_t::INFEASIBLE) { @@ -968,9 +1018,11 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } return mip_status_t::UNBOUNDED; } + if (root_status == lp_status_t::TIME_LIMIT) { - settings_.log.printf("Time limit reached. Stopping the solver...\n"); - return mip_status_t::TIME_LIMIT; + status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, -inf); + return status_; } set_uninitialized_steepest_edge_norms(edge_norms_); @@ -1011,6 +1063,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.log.printf("Optimal solution found at root node. Objective %.16e. Time %.2f.\n", compute_user_objective(original_lp_, root_objective_), toc(stats_.start_time)); + if (settings_.solution_callback != nullptr) { settings_.solution_callback(solution.x, solution.objective); } @@ -1033,8 +1086,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut pc_); if (toc(stats_.start_time) > settings_.time_limit) { - settings_.log.printf("Time limit reached. Stopping the solver...\n"); - return mip_status_t::TIME_LIMIT; + status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, root_objective_); + return status_; } // Choose variable to branch on @@ -1045,6 +1099,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut search_tree.graphviz_node(&search_tree.root, "lower bound", root_objective_); search_tree.branch( &search_tree.root, branch_var, root_relax_soln_.x[branch_var], root_vstatus_, original_lp_); + stats_.nodes_unexplored = 2; settings_.log.printf("Exploring the B&B tree using %d threads (%d diving threads)\n", settings_.num_threads, @@ -1084,50 +1139,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } } - if (status_ == mip_status_t::TIME_LIMIT) { - settings_.log.printf("Time limit reached. Stopping the solver...\n"); - } - f_t lower_bound = heap_.size() > 0 ? heap_.top()->lower_bound : search_tree.get_lower_bound(); - f_t upper_bound = get_upper_bound(); - f_t gap = upper_bound - lower_bound; - f_t obj = compute_user_objective(original_lp_, upper_bound); - f_t user_lower = compute_user_objective(original_lp_, lower_bound); - f_t gap_rel = user_relative_gap(original_lp_, upper_bound, lower_bound); - - settings_.log.printf( - "Explored %d nodes in %.2fs.\n", stats_.nodes_explored, toc(stats_.start_time)); - settings_.log.printf("Absolute Gap %e Objective %.16e Lower Bound %.16e\n", gap, obj, user_lower); - - if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) { - status_ = mip_status_t::OPTIMAL; - if (gap > 0 && gap <= settings_.absolute_mip_gap_tol) { - settings_.log.printf("Optimal solution found within absolute MIP gap tolerance (%.1e)\n", - settings_.absolute_mip_gap_tol); - } else if (gap > 0 && gap_rel <= settings_.relative_mip_gap_tol) { - settings_.log.printf("Optimal solution found within relative MIP gap tolerance (%.1e)\n", - settings_.relative_mip_gap_tol); - } else { - settings_.log.printf("Optimal solution found.\n"); - } - if (settings_.heuristic_preemption_callback != nullptr) { - settings_.heuristic_preemption_callback(); - } - } - - if (get_heap_size() == 0 && upper_bound == inf) { - settings_.log.printf("Integer infeasible.\n"); - status_ = mip_status_t::INFEASIBLE; - if (settings_.heuristic_preemption_callback != nullptr) { - settings_.heuristic_preemption_callback(); - } - } - - uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); - solution.objective = incumbent_.objective; - solution.lower_bound = lower_bound; - solution.nodes_explored = stats_.nodes_explored; - solution.simplex_iterations = stats_.total_lp_iters; + set_final_solution(solution, lower_bound); return status_; } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 82c0d9477a..a44cdb48e9 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -173,6 +173,9 @@ class branch_and_bound_t { // Global status omp_atomic_t status_; + // Set the final solution. + void set_final_solution(mip_solution_t& solution, f_t lower_bound); + // Update the incumbent solution with the new feasible solution. // found during branch and bound. void add_feasible_solution(f_t leaf_objective, From 3489edfcbb17ca5bb8918b8d21dd17b88c999179 Mon Sep 17 00:00:00 2001 From: nicolas Date: Sat, 27 Sep 2025 15:37:28 +0200 Subject: [PATCH 16/33] fixed styling Signed-off-by: nicolas --- build.sh | 10 +++++++--- cpp/src/dual_simplex/branch_and_bound.cpp | 5 ++--- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/build.sh b/build.sh index a5f57f9c9a..228290fae1 100755 --- a/build.sh +++ b/build.sh @@ -206,7 +206,7 @@ if hasArg -v; then VERBOSE_FLAG="-v" fi if hasArg -g; then - BUILD_TYPE=Debug + BUILD_TYPE=RelWithDebInfo DEFINE_ASSERT=true fi if hasArg -a; then @@ -326,7 +326,8 @@ if buildAll || hasArg libmps_parser; then cmake -DDEFINE_ASSERT=${DEFINE_ASSERT} \ -DCMAKE_INSTALL_PREFIX="${INSTALL_PREFIX}" \ "${CACHE_ARGS[@]}" \ - "${REPODIR}"/cpp/libmps_parser/ + -DCMAKE_EXPORT_COMPILE_COMMANDS=ON \ + "${REPODIR}"/cpp/libmps_parser/ if hasArg -n; then cmake --build "${LIBMPS_PARSER_BUILD_DIR}" ${VERBOSE_FLAG} @@ -348,13 +349,16 @@ if buildAll || hasArg libcuopt; then -DCMAKE_CUDA_ARCHITECTURES=${CUOPT_CMAKE_CUDA_ARCHITECTURES} \ -DDISABLE_DEPRECATION_WARNING=${BUILD_DISABLE_DEPRECATION_WARNING} \ -DCMAKE_BUILD_TYPE=${BUILD_TYPE} \ - -DFETCH_RAPIDS=${FETCH_RAPIDS} \ + -DCMAKE_EXPORT_COMPILE_COMMANDS=ON \ + -DFETCH_RAPIDS=${FETCH_RAPIDS} \ -DBUILD_LP_ONLY=${BUILD_LP_ONLY} \ -DBUILD_SANITIZER=${BUILD_SANITIZER} \ -DSKIP_C_PYTHON_ADAPTERS=${SKIP_C_PYTHON_ADAPTERS} \ -DBUILD_TESTS=$((1 - ${SKIP_TESTS_BUILD})) \ -DSKIP_ROUTING_BUILD=${SKIP_ROUTING_BUILD} \ -DWRITE_FATBIN=${WRITE_FATBIN} \ + -DHOST_LINEINFO=${HOST_LINEINFO} \ + -DBUILD_MIP_BENCHMARKS=ON \ "${CACHE_ARGS[@]}" \ "${EXTRA_CMAKE_ARGS[@]}" \ "${REPODIR}"/cpp diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 9fbbcbefab..c96d39044b 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -791,7 +791,6 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear i_t tid = omp_get_thread_num(); while (stack.size() > 0 && status_ == mip_status_t::RUNNING) { - repair_heuristic_solutions(); mip_node_t* node_ptr = stack.front(); @@ -946,7 +945,7 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr } else if (node_status == node_status_t::HAS_CHILDREN) { auto [first, second] = child_selection(node_ptr); - + if (dive_queue_.size() < 4 * settings_.num_diving_threads) { mutex_dive_queue_.lock(); dive_queue_.push(second->detach_copy()); @@ -1135,7 +1134,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut if (omp_get_thread_num() < settings_.num_bfs_threads) { best_first_thread(search_tree, leaf_problem, Arow); - + } else { diving_thread(leaf_problem, Arow); } From 2b62d05ee9e2d3bb6f52bbdee857c1276553783e Mon Sep 17 00:00:00 2001 From: nicolas Date: Sat, 27 Sep 2025 15:41:45 +0200 Subject: [PATCH 17/33] undo changes in the build script --- build.sh | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/build.sh b/build.sh index 228290fae1..a5f57f9c9a 100755 --- a/build.sh +++ b/build.sh @@ -206,7 +206,7 @@ if hasArg -v; then VERBOSE_FLAG="-v" fi if hasArg -g; then - BUILD_TYPE=RelWithDebInfo + BUILD_TYPE=Debug DEFINE_ASSERT=true fi if hasArg -a; then @@ -326,8 +326,7 @@ if buildAll || hasArg libmps_parser; then cmake -DDEFINE_ASSERT=${DEFINE_ASSERT} \ -DCMAKE_INSTALL_PREFIX="${INSTALL_PREFIX}" \ "${CACHE_ARGS[@]}" \ - -DCMAKE_EXPORT_COMPILE_COMMANDS=ON \ - "${REPODIR}"/cpp/libmps_parser/ + "${REPODIR}"/cpp/libmps_parser/ if hasArg -n; then cmake --build "${LIBMPS_PARSER_BUILD_DIR}" ${VERBOSE_FLAG} @@ -349,16 +348,13 @@ if buildAll || hasArg libcuopt; then -DCMAKE_CUDA_ARCHITECTURES=${CUOPT_CMAKE_CUDA_ARCHITECTURES} \ -DDISABLE_DEPRECATION_WARNING=${BUILD_DISABLE_DEPRECATION_WARNING} \ -DCMAKE_BUILD_TYPE=${BUILD_TYPE} \ - -DCMAKE_EXPORT_COMPILE_COMMANDS=ON \ - -DFETCH_RAPIDS=${FETCH_RAPIDS} \ + -DFETCH_RAPIDS=${FETCH_RAPIDS} \ -DBUILD_LP_ONLY=${BUILD_LP_ONLY} \ -DBUILD_SANITIZER=${BUILD_SANITIZER} \ -DSKIP_C_PYTHON_ADAPTERS=${SKIP_C_PYTHON_ADAPTERS} \ -DBUILD_TESTS=$((1 - ${SKIP_TESTS_BUILD})) \ -DSKIP_ROUTING_BUILD=${SKIP_ROUTING_BUILD} \ -DWRITE_FATBIN=${WRITE_FATBIN} \ - -DHOST_LINEINFO=${HOST_LINEINFO} \ - -DBUILD_MIP_BENCHMARKS=ON \ "${CACHE_ARGS[@]}" \ "${EXTRA_CMAKE_ARGS[@]}" \ "${REPODIR}"/cpp From b580e2cd14aa3b4cf860cd455af208f9090c3a30 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 29 Sep 2025 22:26:40 +0200 Subject: [PATCH 18/33] fix incorrect node order in the logs in the ramp up phase --- cpp/src/dual_simplex/branch_and_bound.cpp | 31 +++++++++++++---------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index c96d39044b..fd5ea3def9 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -728,21 +728,24 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* f_t now = toc(stats_.start_time); - if (nodes_explored % 1000 == 0 || - (upper_bound - lower_bound) < 10 * settings_.absolute_mip_gap_tol || nodes_explored < 1000) { - f_t obj = compute_user_objective(original_lp_, upper_bound); - f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); - std::string gap_user = user_mip_gap(obj, user_lower); + if (omp_get_thread_num() == 0) { + if (nodes_explored % 1000 == 0 || + (upper_bound - lower_bound) < 10 * settings_.absolute_mip_gap_tol || + nodes_explored < 1000) { + f_t obj = compute_user_objective(original_lp_, upper_bound); + f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); + std::string gap_user = user_mip_gap(obj, user_lower); - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - nodes_explored, - nodes_unexplored, - obj, - user_lower, - node->depth, - nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, - gap_user.c_str(), - now); + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + nodes_explored, + nodes_unexplored, + obj, + user_lower, + node->depth, + nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, + gap_user.c_str(), + now); + } } if (toc(stats_.start_time) > settings_.time_limit) { From 0cd7cd69119f177ad6c22418439c95f45100e2e6 Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 30 Sep 2025 11:09:40 +0200 Subject: [PATCH 19/33] fixed incomplete tree after an early stop --- cpp/src/dual_simplex/branch_and_bound.cpp | 57 +++++++++++++++-------- cpp/src/dual_simplex/branch_and_bound.hpp | 3 -- 2 files changed, 38 insertions(+), 22 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index fd5ea3def9..ac82c41b4e 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -256,14 +256,6 @@ i_t branch_and_bound_t::get_heap_size() return size; } -template -bool branch_and_bound_t::check_gap(f_t lower_bound, f_t upper_bound) -{ - f_t gap_rel = user_relative_gap(original_lp_, upper_bound, lower_bound); - f_t gap_abs = upper_bound - lower_bound; - return gap_rel < settings_.relative_mip_gap_tol || gap_abs < settings_.absolute_mip_gap_tol; -} - template void branch_and_bound_t::set_new_solution(const std::vector& solution) { @@ -714,13 +706,15 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* f_t lower_bound = node->lower_bound; f_t upper_bound = get_upper_bound(); + f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); + f_t abs_gap = upper_bound - lower_bound; i_t nodes_explored = 0; i_t nodes_unexplored = 0; nodes_explored = (stats_.nodes_explored++); nodes_unexplored = (stats_.nodes_unexplored--); - if (check_gap(lower_bound, upper_bound)) { + if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree->graphviz_node(node, "cutoff", node->lower_bound); search_tree->update_tree(node, node_status_t::FATHOMED); return; @@ -729,8 +723,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* f_t now = toc(stats_.start_time); if (omp_get_thread_num() == 0) { - if (nodes_explored % 1000 == 0 || - (upper_bound - lower_bound) < 10 * settings_.absolute_mip_gap_tol || + if (nodes_explored % 1000 == 0 || abs_gap < 10 * settings_.absolute_mip_gap_tol || nodes_explored < 1000) { f_t obj = compute_user_objective(original_lp_, upper_bound); f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); @@ -801,12 +794,14 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear f_t lower_bound = node_ptr->lower_bound; f_t upper_bound = get_upper_bound(); + f_t abs_gap = upper_bound - lower_bound; + f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); lower_bounds_[tid] = lower_bound; nodes_explored = stats_.nodes_explored++; nodes_unexplored = stats_.nodes_unexplored--; - if (check_gap(lower_bound, upper_bound)) { + if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree.graphviz_node(node_ptr, "cutoff", node_ptr->lower_bound); search_tree.update_tree(node_ptr, node_status_t::FATHOMED); continue; @@ -814,8 +809,7 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear f_t now = toc(stats_.start_time); f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); - if (((nodes_explored % 1000 == 0 || - (upper_bound - lower_bound) < 10 * settings_.absolute_mip_gap_tol || + if (((nodes_explored % 1000 == 0 || abs_gap < 10 * settings_.absolute_mip_gap_tol || nodes_explored < 1000) && (time_since_log >= 1)) || (time_since_log > 60) || now > settings_.time_limit) { @@ -846,7 +840,7 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_status_t::TIME_LIMIT; - return; + break; } else if (node_status == node_status_t::HAS_CHILDREN) { if (stack.size() > 0) { @@ -871,6 +865,14 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear stack.push_front(first); } } + + if (stack.size() != 0) + printf("%d : stack size = %ld, status = %d, node id = %d, node lb = %f\n", + omp_get_thread_num(), + stack.size(), + (int)status_.load(), + stack.front()->node_id, + stack.front()->lower_bound); } template @@ -878,6 +880,11 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se lp_problem_t& leaf_problem, csc_matrix_t& Arow) { + f_t lower_bound = -inf; + f_t upper_bound = inf; + f_t abs_gap = inf; + f_t rel_gap = inf; + do { mip_node_t* node_ptr = nullptr; @@ -905,11 +912,20 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se active_subtrees_--; } - } while (status_ == mip_status_t::RUNNING && !check_gap(get_lower_bound(), get_upper_bound()) && + lower_bound = get_lower_bound(); + upper_bound = get_upper_bound(); + abs_gap = upper_bound - lower_bound; + rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); + + } while (status_ == mip_status_t::RUNNING && abs_gap > settings_.absolute_mip_gap_tol && + rel_gap > settings_.relative_mip_gap_tol && (active_subtrees_ > 0 || get_heap_size() > 0)); - // Check if the solver exited naturally, or was due to a timeout or numerical error. - if (status_ == mip_status_t::RUNNING) { status_ = mip_status_t::COMPLETED; } + // Check if it is the last thread that exited the loop and no + // timeout or numerical error has happen. + if (status_ == mip_status_t::RUNNING && active_subtrees_ == 0) { + status_ = mip_status_t::COMPLETED; + } } template @@ -937,8 +953,11 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr mip_node_t* node_ptr = stack.front(); stack.pop_front(); f_t upper_bound = get_upper_bound(); + f_t rel_gap = user_relative_gap(original_lp_, upper_bound, node_ptr->lower_bound); - if (check_gap(node_ptr->lower_bound, upper_bound)) { continue; } + if (node_ptr->lower_bound > upper_bound || rel_gap > settings_.relative_mip_gap_tol) { + continue; + } node_status_t node_status = solve_node_lp(subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index a44cdb48e9..8887dfa6da 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -106,9 +106,6 @@ class branch_and_bound_t { f_t get_lower_bound(); i_t get_heap_size(); - // Check if the gap and relative gap are below the tolerance. - bool check_gap(f_t lower_bound, f_t upper_bound); - // The main entry routine. Returns the solver status and populates solution with the incumbent. mip_status_t solve(mip_solution_t& solution); From 0ffbf21d88f8203d8070a91978c411cc8a46ee13 Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 30 Sep 2025 11:40:56 +0200 Subject: [PATCH 20/33] renamed variable --- cpp/src/dual_simplex/branch_and_bound.cpp | 12 ++++++------ cpp/src/dual_simplex/branch_and_bound.hpp | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index ac82c41b4e..3dafa41971 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -240,8 +240,8 @@ f_t branch_and_bound_t::get_lower_bound() if (heap_.size() > 0) { lower_bound = heap_.top()->lower_bound; } mutex_heap_.unlock(); - for (i_t i = 0; i < lower_bounds_.size(); ++i) { - lower_bound = std::min(lower_bounds_[i].load(), lower_bound); + for (i_t i = 0; i < local_lower_bounds_.size(); ++i) { + lower_bound = std::min(local_lower_bounds_[i].load(), lower_bound); } return lower_bound; @@ -797,9 +797,9 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear f_t abs_gap = upper_bound - lower_bound; f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); - lower_bounds_[tid] = lower_bound; - nodes_explored = stats_.nodes_explored++; - nodes_unexplored = stats_.nodes_unexplored--; + local_lower_bounds_[tid] = lower_bound; + nodes_explored = stats_.nodes_explored++; + nodes_unexplored = stats_.nodes_unexplored--; if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree.graphviz_node(node_ptr, "cutoff", node_ptr->lower_bound); @@ -1050,7 +1050,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut set_uninitialized_steepest_edge_norms(edge_norms_); root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - lower_bounds_.assign(settings_.num_bfs_threads, root_objective_); + local_lower_bounds_.assign(settings_.num_bfs_threads, root_objective_); if (settings_.set_simplex_solution_callback != nullptr) { std::vector original_x; diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 8887dfa6da..d4b7a29671 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -122,7 +122,7 @@ class branch_and_bound_t { std::vector var_types_; // Local lower bounds for each thread - std::vector> lower_bounds_; + std::vector> local_lower_bounds_; // Mutex for upper bound omp_mutex_t mutex_upper_; From ff9942037f05b89f46a087468e38aa15d4c20fc8 Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 30 Sep 2025 15:45:05 +0200 Subject: [PATCH 21/33] removed debug info --- cpp/src/dual_simplex/branch_and_bound.cpp | 8 -------- 1 file changed, 8 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 3dafa41971..1929a147fa 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -865,14 +865,6 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear stack.push_front(first); } } - - if (stack.size() != 0) - printf("%d : stack size = %ld, status = %d, node id = %d, node lb = %f\n", - omp_get_thread_num(), - stack.size(), - (int)status_.load(), - stack.front()->node_id, - stack.front()->lower_bound); } template From a9fe075df9f2331c649d0818bdc0bf06f5109237 Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 30 Sep 2025 20:38:52 +0200 Subject: [PATCH 22/33] changed ramp up termination criteria. restrict log to the thread 0. --- cpp/src/dual_simplex/branch_and_bound.cpp | 104 +++++++++++++--------- cpp/src/dual_simplex/branch_and_bound.hpp | 2 +- 2 files changed, 64 insertions(+), 42 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 1929a147fa..5bf734c3da 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -699,7 +699,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* mip_node_t* node, lp_problem_t& leaf_problem, csc_matrix_t& Arow, - i_t max_depth) + i_t initial_heap_size) { if (status_ != mip_status_t::RUNNING) { return; } repair_heuristic_solutions(); @@ -755,12 +755,13 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* } else if (node_status == node_status_t::HAS_CHILDREN) { stats_.nodes_unexplored += 2; - if (node->depth < max_depth) { + if (stats_.nodes_unexplored < initial_heap_size) { #pragma omp task - exploration_ramp_up(search_tree, node->get_down_child(), leaf_problem, Arow, max_depth); + exploration_ramp_up( + search_tree, node->get_down_child(), leaf_problem, Arow, initial_heap_size); #pragma omp task - exploration_ramp_up(search_tree, node->get_up_child(), leaf_problem, Arow, max_depth); + exploration_ramp_up(search_tree, node->get_up_child(), leaf_problem, Arow, initial_heap_size); } else { mutex_heap_.lock(); @@ -781,10 +782,11 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear stack.push_front(start_node); // It contains the lower bound of the parent for now - f_t last_log = 0; - i_t nodes_explored = 0; - i_t nodes_unexplored = 0; - i_t tid = omp_get_thread_num(); + f_t last_log = 0; + i_t nodes_explored = 0; + i_t nodes_unexplored = 0; + i_t tid = omp_get_thread_num(); + i_t nodes_since_last_log = 0; while (stack.size() > 0 && status_ == mip_status_t::RUNNING) { repair_heuristic_solutions(); @@ -797,9 +799,16 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear f_t abs_gap = upper_bound - lower_bound; f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); + // This is based on three assumptions: + // - The stack only contains sibling nodes, i.e., the current node and its siblings, if it + // exists + // - The current node and its siblings uses the lower bound of the parent before solving the LP + // relaxation + // - The lower bound of the parent is lower or equal to its children local_lower_bounds_[tid] = lower_bound; nodes_explored = stats_.nodes_explored++; nodes_unexplored = stats_.nodes_unexplored--; + nodes_since_last_log++; if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree.graphviz_node(node_ptr, "cutoff", node_ptr->lower_bound); @@ -807,25 +816,30 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear continue; } - f_t now = toc(stats_.start_time); - f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); - if (((nodes_explored % 1000 == 0 || abs_gap < 10 * settings_.absolute_mip_gap_tol || - nodes_explored < 1000) && - (time_since_log >= 1)) || - (time_since_log > 60) || now > settings_.time_limit) { - f_t obj = compute_user_objective(original_lp_, upper_bound); - f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); - std::string gap_user = user_mip_gap(obj, user_lower); - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - nodes_explored, - nodes_unexplored, - obj, - user_lower, - node_ptr->depth, - nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, - gap_user.c_str(), - now); - last_log = tic(); + f_t now = toc(stats_.start_time); + + if (tid == 0) { + f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); + + if (((nodes_since_last_log % 1000 == 0 || abs_gap < 10 * settings_.absolute_mip_gap_tol || + nodes_since_last_log < 1000) && + (time_since_log >= 1)) || + (time_since_log > 60) || now > settings_.time_limit) { + f_t obj = compute_user_objective(original_lp_, upper_bound); + f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); + std::string gap_user = user_mip_gap(obj, user_lower); + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + nodes_explored, + nodes_unexplored, + obj, + user_lower, + node_ptr->depth, + nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, + gap_user.c_str(), + now); + last_log = tic(); + nodes_since_last_log = 0; + } } if (toc(stats_.start_time) > settings_.time_limit) { @@ -847,6 +861,12 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear mip_node_t* node = stack.back(); stack.pop_back(); + // The order here matters. We want to create a copy of the node + // before adding to the global heap. Otherwise, + // some thread may consume the node (possibly fathoming it) + // before we had the chance to add to the diving queue. + // This lead to a SIGSEGV. Although, in this case, it + // would be better if we discard the node instead. if (get_heap_size() > settings_.num_bfs_threads) { mutex_dive_queue_.lock(); dive_queue_.push(node->detach_copy()); @@ -959,16 +979,16 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr } else if (node_status == node_status_t::HAS_CHILDREN) { auto [first, second] = child_selection(node_ptr); + stack.push_front(second); + stack.push_front(first); if (dive_queue_.size() < 4 * settings_.num_diving_threads) { mutex_dive_queue_.lock(); - dive_queue_.push(second->detach_copy()); + mip_node_t* new_node = stack.back(); + stack.pop_back(); + dive_queue_.push(new_node->detach_copy()); mutex_dive_queue_.unlock(); - } else { - stack.push_front(second); } - - stack.push_front(first); } } } @@ -1133,24 +1153,26 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut #pragma omp master { - i_t max_depth = std::ceil(std::log2(settings_.num_threads)); + auto down_child = search_tree.root.get_down_child(); + auto up_child = search_tree.root.get_up_child(); + i_t initial_size = 2 * settings_.num_threads; #pragma omp task - exploration_ramp_up( - &search_tree, search_tree.root.get_down_child(), leaf_problem, Arow, max_depth); + exploration_ramp_up(&search_tree, down_child, leaf_problem, Arow, initial_size); #pragma omp task - exploration_ramp_up( - &search_tree, search_tree.root.get_up_child(), leaf_problem, Arow, max_depth); + exploration_ramp_up(&search_tree, up_child, leaf_problem, Arow, initial_size); } #pragma omp barrier - if (omp_get_thread_num() < settings_.num_bfs_threads) { - best_first_thread(search_tree, leaf_problem, Arow); + if (status_ == mip_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { + if (omp_get_thread_num() < settings_.num_bfs_threads) { + best_first_thread(search_tree, leaf_problem, Arow); - } else { - diving_thread(leaf_problem, Arow); + } else { + diving_thread(leaf_problem, Arow); + } } } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index d4b7a29671..3708bf2ced 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -189,7 +189,7 @@ class branch_and_bound_t { mip_node_t* node, lp_problem_t& leaf_problem, csc_matrix_t& Arow, - i_t max_depth); + i_t initial_heap_size); // Explore the search tree using the best-first search with plunging strategy. void explore_subtree(search_tree_t& search_tree, From 82186e6a1faa26c7815836c4884a4f693d47facf Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 1 Oct 2025 17:52:13 +0200 Subject: [PATCH 23/33] address Chris' feedback --- cpp/src/dual_simplex/branch_and_bound.cpp | 276 ++++++++++------------ cpp/src/dual_simplex/branch_and_bound.hpp | 58 ++--- cpp/src/dual_simplex/mip_node.hpp | 27 +-- cpp/src/dual_simplex/pseudo_costs.cpp | 13 +- cpp/src/dual_simplex/pseudo_costs.hpp | 11 +- 5 files changed, 185 insertions(+), 200 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 5bf734c3da..3d4a206e04 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -212,7 +212,7 @@ branch_and_bound_t::branch_and_bound_t( incumbent_(1), root_relax_soln_(1, 1), pc_(1), - status_(mip_status_t::UNSET) + status_(mip_exploration_status_t::UNSET) { stats_.start_time = tic(); convert_user_problem(original_problem_, settings_, original_lp_, new_slacks_); @@ -294,7 +294,7 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu mutex_upper_.unlock(); if (is_feasible) { - if (status_ == mip_status_t::RUNNING) { + if (status_ == mip_exploration_status_t::RUNNING) { f_t user_obj = compute_user_objective(original_lp_, obj); f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); std::string gap = user_mip_gap(user_obj, user_lower); @@ -433,11 +433,14 @@ void branch_and_bound_t::repair_heuristic_solutions() } template -void branch_and_bound_t::set_final_solution(mip_solution_t& solution, - f_t lower_bound) +mip_status_t branch_and_bound_t::set_final_solution(mip_solution_t& solution, + f_t lower_bound) { - if (status_ == mip_status_t::TIME_LIMIT) { + mip_status_t mip_status = mip_status_t::UNSET; + + if (status_ == mip_exploration_status_t::TIME_LIMIT) { settings_.log.printf("Time limit reached. Stopping the solver...\n"); + mip_status = mip_status_t::TIME_LIMIT; } f_t upper_bound = get_upper_bound(); @@ -451,7 +454,7 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& settings_.log.printf("Absolute Gap %e Objective %.16e Lower Bound %.16e\n", gap, obj, user_lower); if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) { - status_ = mip_status_t::OPTIMAL; + mip_status = mip_status_t::OPTIMAL; if (gap > 0 && gap <= settings_.absolute_mip_gap_tol) { settings_.log.printf("Optimal solution found within absolute MIP gap tolerance (%.1e)\n", settings_.absolute_mip_gap_tol); @@ -468,7 +471,7 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& if (stats_.nodes_explored > 0 && stats_.nodes_unexplored == 0 && upper_bound == inf) { settings_.log.printf("Integer infeasible.\n"); - status_ = mip_status_t::INFEASIBLE; + mip_status = mip_status_t::INFEASIBLE; if (settings_.heuristic_preemption_callback != nullptr) { settings_.heuristic_preemption_callback(); } @@ -479,6 +482,8 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& solution.lower_bound = lower_bound; solution.nodes_explored = stats_.nodes_explored; solution.simplex_iterations = stats_.total_lp_iters; + + return mip_status; } template @@ -540,20 +545,35 @@ branch_and_bound_t::child_selection(mip_node_t* node_ptr) } template -dual::status_t branch_and_bound_t::node_dual_simplex( - i_t leaf_id, +node_status_t branch_and_bound_t::solve_node_lp_and_update_tree( + search_tree_t& search_tree, + mip_node_t* node_ptr, lp_problem_t& leaf_problem, - std::vector& leaf_vstatus, - lp_solution_t& leaf_solution, - std::vector& bounds_changed, csc_matrix_t& Arow, f_t upper_bound, - logger_t& log) + logger_t& log, + char symbol) { - i_t node_iter = 0; + f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; + + std::vector& leaf_vstatus = node_ptr->vstatus; + lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); assert(leaf_vstatus.size() == leaf_problem.num_cols); - f_t lp_start_time = tic(); - std::vector leaf_edge_norms = edge_norms_; // = node.steepest_edge_norms; + + // Set the correct bounds for the leaf problem + leaf_problem.lower = original_lp_.lower; + leaf_problem.upper = original_lp_.upper; + + std::vector bounds_changed(leaf_problem.num_cols, false); + // Technically, we can get the already strengthened bounds from the node/parent instead of + // getting it from the original problem and re-strengthening. But this requires storing + // two vectors at each node and potentially cause memory issues + node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); + + i_t node_iter = 0; + f_t lp_start_time = tic(); + std::vector leaf_edge_norms = edge_norms_; // = node.steepest_edge_norms; + simplex_solver_settings_t lp_settings = settings_; lp_settings.set_log(false); lp_settings.cut_off = upper_bound + settings_.dual_tol; @@ -578,7 +598,7 @@ dual::status_t branch_and_bound_t::node_dual_simplex( leaf_edge_norms); if (lp_status == dual::status_t::NUMERICAL) { - log.printf("Numerical issue node %d. Resolving from scratch.\n", leaf_id); + log.printf("Numerical issue node %d. Resolving from scratch.\n", node_ptr->node_id); lp_status_t second_status = solve_linear_program_advanced( leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms); lp_status = convert_lp_status_to_dual_status(second_status); @@ -587,48 +607,11 @@ dual::status_t branch_and_bound_t::node_dual_simplex( stats_.total_lp_solve_time += toc(lp_start_time); stats_.total_lp_iters += node_iter; - return lp_status; -} - -template -node_status_t branch_and_bound_t::solve_node_lp(search_tree_t& search_tree, - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log, - char symbol) -{ - logger_t pc_log; - pc_log.log = false; - - f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; - - std::vector& leaf_vstatus = node_ptr->vstatus; - lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); - - // Set the correct bounds for the leaf problem - leaf_problem.lower = original_lp_.lower; - leaf_problem.upper = original_lp_.upper; - std::vector bounds_changed(leaf_problem.num_cols, false); - // Technically, we can get the already strengthened bounds from the node/parent instead of - // getting it from the original problem and re-strengthening. But this requires storing - // two vectors at each node and potentially cause memory issues - node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); - - dual::status_t lp_status = node_dual_simplex(node_ptr->node_id, - leaf_problem, - leaf_vstatus, - leaf_solution, - bounds_changed, - Arow, - upper_bound, - log); if (lp_status == dual::status_t::DUAL_UNBOUNDED) { // Node was infeasible. Do not branch node_ptr->lower_bound = inf; - search_tree.graphviz_node(node_ptr, "infeasible", 0.0); + search_tree.graphviz_node(log, node_ptr, "infeasible", 0.0); search_tree.update_tree(node_ptr, node_status_t::INFEASIBLE); return node_status_t::INFEASIBLE; @@ -636,7 +619,7 @@ node_status_t branch_and_bound_t::solve_node_lp(search_tree_tlower_bound = upper_bound; f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); - search_tree.graphviz_node(node_ptr, "cut off", leaf_objective); + search_tree.graphviz_node(log, node_ptr, "cut off", leaf_objective); search_tree.update_tree(node_ptr, node_status_t::FATHOMED); return node_status_t::FATHOMED; @@ -646,46 +629,41 @@ node_status_t branch_and_bound_t::solve_node_lp(search_tree_tlower_bound = leaf_objective; if (leaf_num_fractional == 0) { // Found a integer feasible solution add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, symbol); - search_tree.graphviz_node(node_ptr, "integer feasible", leaf_objective); + search_tree.graphviz_node(log, node_ptr, "integer feasible", leaf_objective); search_tree.update_tree(node_ptr, node_status_t::INTEGER_FEASIBLE); return node_status_t::INTEGER_FEASIBLE; } else if (leaf_objective <= upper_bound + abs_fathom_tol) { // Choose fractional variable to branch on - mutex_pc_.lock(); - const i_t branch_var = pc_.variable_selection( - leaf_fractional, leaf_solution.x, leaf_problem.lower, leaf_problem.upper, pc_log); - mutex_pc_.unlock(); + const i_t branch_var = + pc_.variable_selection(leaf_fractional, leaf_solution.x, lp_settings.log); assert(leaf_vstatus.size() == leaf_problem.num_cols); search_tree.branch( - node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus, original_lp_); + node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus, original_lp_, log); node_ptr->status = node_status_t::HAS_CHILDREN; return node_status_t::HAS_CHILDREN; } else { - search_tree.graphviz_node(node_ptr, "fathomed", leaf_objective); + search_tree.graphviz_node(log, node_ptr, "fathomed", leaf_objective); search_tree.update_tree(node_ptr, node_status_t::FATHOMED); return node_status_t::FATHOMED; } } else if (lp_status == dual::status_t::TIME_LIMIT) { - search_tree.graphviz_node(node_ptr, "timeout", 0.0); + search_tree.graphviz_node(log, node_ptr, "timeout", 0.0); search_tree.update_tree(node_ptr, node_status_t::TIME_LIMIT); return node_status_t::TIME_LIMIT; } else { - search_tree.graphviz_node(node_ptr, "numerical", 0.0); + search_tree.graphviz_node(log, node_ptr, "numerical", 0.0); log.printf("Encountered LP status %d on node %d. This indicates a numerical issue.\n", lp_status, node_ptr->node_id); @@ -701,7 +679,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* csc_matrix_t& Arow, i_t initial_heap_size) { - if (status_ != mip_status_t::RUNNING) { return; } + if (status_ != mip_exploration_status_t::RUNNING) { return; } repair_heuristic_solutions(); f_t lower_bound = node->lower_bound; @@ -715,7 +693,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* nodes_unexplored = (stats_.nodes_unexplored--); if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { - search_tree->graphviz_node(node, "cutoff", node->lower_bound); + search_tree->graphviz_node(settings_.log, node, "cutoff", node->lower_bound); search_tree->update_tree(node, node_status_t::FATHOMED); return; } @@ -723,10 +701,14 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* f_t now = toc(stats_.start_time); if (omp_get_thread_num() == 0) { - if (nodes_explored % 1000 == 0 || abs_gap < 10 * settings_.absolute_mip_gap_tol || - nodes_explored < 1000) { + stats_.nodes_since_last_log++; + f_t time_since_last_log = stats_.last_log == 0 ? 1.0 : toc(stats_.last_log); + + if (((stats_.nodes_since_last_log > 10 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && + (time_since_last_log >= 1)) || + (time_since_last_log > 30) || now > settings_.time_limit) { f_t obj = compute_user_objective(original_lp_, upper_bound); - f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); + f_t user_lower = compute_user_objective(original_lp_, root_objective_); std::string gap_user = user_mip_gap(obj, user_lower); settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", @@ -742,14 +724,14 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* } if (toc(stats_.start_time) > settings_.time_limit) { - status_ = mip_status_t::TIME_LIMIT; + status_ = mip_exploration_status_t::TIME_LIMIT; return; } - node_status_t node_status = - solve_node_lp(*search_tree, node, leaf_problem, Arow, upper_bound, settings_.log, 'B'); + node_status_t node_status = solve_node_lp_and_update_tree( + *search_tree, node, leaf_problem, Arow, upper_bound, settings_.log, 'B'); if (node_status == node_status_t::TIME_LIMIT) { - status_ = mip_status_t::TIME_LIMIT; + status_ = mip_exploration_status_t::TIME_LIMIT; return; } else if (node_status == node_status_t::HAS_CHILDREN) { @@ -773,7 +755,8 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* } template -void branch_and_bound_t::explore_subtree(search_tree_t& search_tree, +void branch_and_bound_t::explore_subtree(i_t id, + search_tree_t& search_tree, mip_node_t* start_node, lp_problem_t& leaf_problem, csc_matrix_t& Arow) @@ -781,14 +764,7 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear std::deque*> stack; stack.push_front(start_node); - // It contains the lower bound of the parent for now - f_t last_log = 0; - i_t nodes_explored = 0; - i_t nodes_unexplored = 0; - i_t tid = omp_get_thread_num(); - i_t nodes_since_last_log = 0; - - while (stack.size() > 0 && status_ == mip_status_t::RUNNING) { + while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { repair_heuristic_solutions(); mip_node_t* node_ptr = stack.front(); @@ -805,26 +781,25 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear // - The current node and its siblings uses the lower bound of the parent before solving the LP // relaxation // - The lower bound of the parent is lower or equal to its children - local_lower_bounds_[tid] = lower_bound; - nodes_explored = stats_.nodes_explored++; - nodes_unexplored = stats_.nodes_unexplored--; - nodes_since_last_log++; + local_lower_bounds_[id] = lower_bound; + i_t nodes_explored = stats_.nodes_explored++; + i_t nodes_unexplored = stats_.nodes_unexplored--; if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { - search_tree.graphviz_node(node_ptr, "cutoff", node_ptr->lower_bound); + search_tree.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound); search_tree.update_tree(node_ptr, node_status_t::FATHOMED); continue; } f_t now = toc(stats_.start_time); - if (tid == 0) { - f_t time_since_log = last_log == 0 ? 1.0 : toc(last_log); + if (id == 0) { + stats_.nodes_since_last_log++; + f_t time_since_last_log = stats_.last_log == 0 ? 1.0 : toc(stats_.last_log); - if (((nodes_since_last_log % 1000 == 0 || abs_gap < 10 * settings_.absolute_mip_gap_tol || - nodes_since_last_log < 1000) && - (time_since_log >= 1)) || - (time_since_log > 60) || now > settings_.time_limit) { + if (((stats_.nodes_since_last_log >= 1000 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && + time_since_last_log >= 1) || + (time_since_last_log > 30) || now > settings_.time_limit) { f_t obj = compute_user_objective(original_lp_, upper_bound); f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); std::string gap_user = user_mip_gap(obj, user_lower); @@ -837,26 +812,28 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, gap_user.c_str(), now); - last_log = tic(); - nodes_since_last_log = 0; + stats_.last_log = tic(); + stats_.nodes_since_last_log = 0; } } if (toc(stats_.start_time) > settings_.time_limit) { - status_ = mip_status_t::TIME_LIMIT; + status_ = mip_exploration_status_t::TIME_LIMIT; break; } - // We are just checking for numerical issues or timeout during the LP solve, otherwise - // lp_status is set to RUNNING. - node_status_t node_status = - solve_node_lp(search_tree, node_ptr, leaf_problem, Arow, upper_bound, settings_.log, 'B'); + node_status_t node_status = solve_node_lp_and_update_tree( + search_tree, node_ptr, leaf_problem, Arow, upper_bound, settings_.log, 'B'); if (node_status == node_status_t::TIME_LIMIT) { - status_ = mip_status_t::TIME_LIMIT; + status_ = mip_exploration_status_t::TIME_LIMIT; break; } else if (node_status == node_status_t::HAS_CHILDREN) { + // The stack should only contain the children of the current parent. + // If the stack size is greater than 0, + // we pop the current node from the stack and place it in the global heap, + // since we are about to add the two children to the stack if (stack.size() > 0) { mip_node_t* node = stack.back(); stack.pop_back(); @@ -888,7 +865,8 @@ void branch_and_bound_t::explore_subtree(search_tree_t& sear } template -void branch_and_bound_t::best_first_thread(search_tree_t& search_tree, +void branch_and_bound_t::best_first_thread(i_t id, + search_tree_t& search_tree, lp_problem_t& leaf_problem, csc_matrix_t& Arow) { @@ -897,7 +875,9 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se f_t abs_gap = inf; f_t rel_gap = inf; - do { + while (status_ == mip_exploration_status_t::RUNNING && abs_gap > settings_.absolute_mip_gap_tol && + rel_gap > settings_.relative_mip_gap_tol && + (active_subtrees_ > 0 || get_heap_size() > 0)) { mip_node_t* node_ptr = nullptr; // If there any node left in the heap, we pop the top node and explore it. @@ -913,14 +893,14 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se if (get_upper_bound() < node_ptr->lower_bound) { // This node was put on the heap earlier but its lower bound is now greater than the // current upper bound - search_tree.graphviz_node(node_ptr, "cutoff", node_ptr->lower_bound); + search_tree.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound); search_tree.update_tree(node_ptr, node_status_t::FATHOMED); active_subtrees_--; continue; } // Best-first search with plunging - explore_subtree(search_tree, node_ptr, leaf_problem, Arow); + explore_subtree(id, search_tree, node_ptr, leaf_problem, Arow); active_subtrees_--; } @@ -928,15 +908,12 @@ void branch_and_bound_t::best_first_thread(search_tree_t& se upper_bound = get_upper_bound(); abs_gap = upper_bound - lower_bound; rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); - - } while (status_ == mip_status_t::RUNNING && abs_gap > settings_.absolute_mip_gap_tol && - rel_gap > settings_.relative_mip_gap_tol && - (active_subtrees_ > 0 || get_heap_size() > 0)); + } // Check if it is the last thread that exited the loop and no // timeout or numerical error has happen. - if (status_ == mip_status_t::RUNNING && active_subtrees_ == 0) { - status_ = mip_status_t::COMPLETED; + if (status_ == mip_exploration_status_t::RUNNING && active_subtrees_ == 0) { + status_ = mip_exploration_status_t::COMPLETED; } } @@ -947,7 +924,8 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr logger_t log; log.log = false; - do { + while (status_ == mip_exploration_status_t::RUNNING && + (active_subtrees_ > 0 || get_heap_size() > 0)) { std::optional> start_node; mutex_dive_queue_.lock(); @@ -957,11 +935,11 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr if (start_node) { if (get_upper_bound() < start_node->lower_bound) { continue; } - search_tree_t subtree(std::move(start_node.value()), settings_.log); + search_tree_t subtree(std::move(start_node.value())); std::deque*> stack; stack.push_front(&subtree.root); - while (stack.size() > 0 && status_ == mip_status_t::RUNNING) { + while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { mip_node_t* node_ptr = stack.front(); stack.pop_front(); f_t upper_bound = get_upper_bound(); @@ -971,8 +949,8 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr continue; } - node_status_t node_status = - solve_node_lp(subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); + node_status_t node_status = solve_node_lp_and_update_tree( + subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); if (node_status == node_status_t::TIME_LIMIT) { break; @@ -992,8 +970,7 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr } } } - - } while (status_ == mip_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)); + } } template @@ -1001,7 +978,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut { logger_t log; log.log = false; - status_ = mip_status_t::UNSET; + status_ = mip_exploration_status_t::UNSET; stats_.total_lp_iters = 0; stats_.nodes_explored = 0; @@ -1054,9 +1031,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (root_status == lp_status_t::TIME_LIMIT) { - status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, -inf); - return status_; + status_ = mip_exploration_status_t::TIME_LIMIT; + return set_final_solution(solution, -inf); } set_uninitialized_steepest_edge_norms(edge_norms_); @@ -1120,19 +1096,21 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut pc_); if (toc(stats_.start_time) > settings_.time_limit) { - status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, root_objective_); - return status_; + status_ = mip_exploration_status_t::TIME_LIMIT; + return set_final_solution(solution, root_objective_); } // Choose variable to branch on - i_t branch_var = pc_.variable_selection( - fractional, root_relax_soln_.x, original_lp_.lower, original_lp_.upper, log); - - search_tree_t search_tree(root_objective_, root_vstatus_, settings_.log); - search_tree.graphviz_node(&search_tree.root, "lower bound", root_objective_); - search_tree.branch( - &search_tree.root, branch_var, root_relax_soln_.x[branch_var], root_vstatus_, original_lp_); + i_t branch_var = pc_.variable_selection(fractional, root_relax_soln_.x, log); + + search_tree_t search_tree(root_objective_, root_vstatus_); + search_tree.graphviz_node(settings_.log, &search_tree.root, "lower bound", root_objective_); + search_tree.branch(&search_tree.root, + branch_var, + root_relax_soln_.x[branch_var], + root_vstatus_, + original_lp_, + log); stats_.nodes_unexplored = 2; settings_.log.printf("Exploring the B&B tree using %d threads (%d diving threads)\n", @@ -1142,7 +1120,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut " | Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap " "| Time |\n"); - status_ = mip_status_t::RUNNING; + status_ = mip_exploration_status_t::RUNNING; #pragma omp parallel num_threads(settings_.num_threads) { @@ -1166,19 +1144,25 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut #pragma omp barrier - if (status_ == mip_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { - if (omp_get_thread_num() < settings_.num_bfs_threads) { - best_first_thread(search_tree, leaf_problem, Arow); +#pragma omp master + { + if (status_ == mip_exploration_status_t::RUNNING && + (active_subtrees_ > 0 || get_heap_size() > 0)) { + for (i_t i = 0; i < settings_.num_bfs_threads; i++) { +#pragma omp task + best_first_thread(i, search_tree, leaf_problem, Arow); + } - } else { - diving_thread(leaf_problem, Arow); + for (i_t i = 0; i < settings_.num_diving_threads; i++) { +#pragma omp task + diving_thread(leaf_problem, Arow); + } } } } - f_t lower_bound = heap_.size() > 0 ? heap_.top()->lower_bound : search_tree.get_lower_bound(); - set_final_solution(solution, lower_bound); - return status_; + f_t lower_bound = heap_.size() > 0 ? heap_.top()->lower_bound : search_tree.root.lower_bound; + return set_final_solution(solution, lower_bound); } #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 3708bf2ced..f5d9323741 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -40,9 +40,16 @@ enum class mip_status_t { TIME_LIMIT = 3, // The solver reached a time limit NODE_LIMIT = 4, // The maximum number of nodes was reached (not implemented) NUMERICAL = 5, // The solver encountered a numerical error - RUNNING = 6, // The solver is currently exploring the tree - COMPLETED = 7, // The solver finished exploring the tree - UNSET = 8, // The status is not set + UNSET = 6, // The status is not set +}; + +enum class mip_exploration_status_t { + UNSET = 0, // The status is not set + TIME_LIMIT = 1, // The solver reached a time limit + NODE_LIMIT = 2, // The maximum number of nodes was reached (not implemented) + NUMERICAL = 3, // The solver encountered a numerical error + RUNNING = 4, // The solver is currently exploring the tree + COMPLETED = 5, // The solver finished exploring the tree }; template @@ -140,6 +147,10 @@ class branch_and_bound_t { omp_atomic_t nodes_explored = 0; omp_atomic_t nodes_unexplored = 0; omp_atomic_t total_lp_iters = 0; + + // This should only be used by the main thread + f_t last_log = 0.0; + i_t nodes_since_last_log = 0; } stats_; // Mutex for repair @@ -154,7 +165,6 @@ class branch_and_bound_t { // Pseudocosts pseudo_costs_t pc_; - omp_mutex_t mutex_pc_; // Heap storing the nodes to be explored. omp_mutex_t mutex_heap_; @@ -168,10 +178,10 @@ class branch_and_bound_t { dive_queue_t dive_queue_; // Global status - omp_atomic_t status_; + omp_atomic_t status_; // Set the final solution. - void set_final_solution(mip_solution_t& solution, f_t lower_bound); + mip_status_t set_final_solution(mip_solution_t& solution, f_t lower_bound); // Update the incumbent solution with the new feasible solution. // found during branch and bound. @@ -184,7 +194,7 @@ class branch_and_bound_t { void repair_heuristic_solutions(); // Ramp-up phase of the solver, where we greedily expand the tree until - // a certain depth is reached. This is done recursively using OpenMP tasks. + // there is enough unexplored nodes. This is done recursively using OpenMP tasks. void exploration_ramp_up(search_tree_t* search_tree, mip_node_t* node, lp_problem_t& leaf_problem, @@ -192,39 +202,31 @@ class branch_and_bound_t { i_t initial_heap_size); // Explore the search tree using the best-first search with plunging strategy. - void explore_subtree(search_tree_t& search_tree, + void explore_subtree(i_t id, + search_tree_t& search_tree, mip_node_t* start_node, lp_problem_t& leaf_problem, csc_matrix_t& Arow); - // Each "main" thread pop a node from the global heap and then perform a plunge + // Each "main" thread pops a node from the global heap and then performs a plunge // (i.e., a shallow dive) into the subtree determined by the node. - void best_first_thread(search_tree_t& search_tree, + void best_first_thread(i_t id, + search_tree_t& search_tree, lp_problem_t& leaf_problem, csc_matrix_t& Arow); - // Each diving thread pop the first node from the dive queue and then perform + // Each diving thread pops the first node from the dive queue and then performs // a deep dive into the subtree determined by the node. void diving_thread(lp_problem_t& leaf_problem, csc_matrix_t& Arow); // Solve the LP relaxation of a leaf node. - node_status_t solve_node_lp(search_tree_t& search_tree, - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log, - char symbol); - - // Solve the LP relaxation of a leaf node using the dual simplex method. - dual::status_t node_dual_simplex(i_t leaf_id, - lp_problem_t& leaf_problem, - std::vector& leaf_vstatus, - lp_solution_t& leaf_solution, - std::vector& bounds_changed, - csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log); + node_status_t solve_node_lp_and_update_tree(search_tree_t& search_tree, + mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + csc_matrix_t& Arow, + f_t upper_bound, + logger_t& log, + char symbol); // Sort the children based on the Martin's criteria. std::pair*, mip_node_t*> child_selection( diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 12f01d3ed1..dccfa4480b 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -262,17 +262,12 @@ class node_compare_t { template class search_tree_t { public: - search_tree_t(f_t root_lower_bound, const std::vector& basis, logger_t& log) - : root(root_lower_bound, basis), num_nodes(0), log(log) + search_tree_t(f_t root_lower_bound, const std::vector& basis) + : root(root_lower_bound, basis), num_nodes(0) { } - search_tree_t(mip_node_t&& node, logger_t& log) - : root(std::move(node)), num_nodes(0), log(log) - { - } - - f_t get_lower_bound() const { return root.lower_bound; } + search_tree_t(mip_node_t&& node) : root(std::move(node)), num_nodes(0) {} void update_tree(mip_node_t* node_ptr, node_status_t status) { @@ -287,7 +282,8 @@ class search_tree_t { const i_t branch_var, const f_t branch_var_val, const std::vector& parent_vstatus, - const lp_problem_t& original_lp) + const lp_problem_t& original_lp, + logger_t& log) { i_t id = num_nodes.fetch_add(2); @@ -295,27 +291,31 @@ class search_tree_t { auto down_child = std::make_unique>( original_lp, parent_node, ++id, branch_var, 0, branch_var_val, parent_vstatus); - graphviz_edge(parent_node, down_child.get(), branch_var, 0, std::floor(branch_var_val)); + graphviz_edge(log, parent_node, down_child.get(), branch_var, 0, std::floor(branch_var_val)); // up child auto up_child = std::make_unique>( original_lp, parent_node, ++id, branch_var, 1, branch_var_val, parent_vstatus); - graphviz_edge(parent_node, up_child.get(), branch_var, 1, std::ceil(branch_var_val)); + graphviz_edge(log, parent_node, up_child.get(), branch_var, 1, std::ceil(branch_var_val)); assert(parent_vstatus.size() == original_lp.num_cols); parent_node->add_children(std::move(down_child), std::move(up_child)); // child pointers moved into the tree } - void graphviz_node(const mip_node_t* node_ptr, const std::string label, const f_t val) + void graphviz_node(logger_t& log, + const mip_node_t* node_ptr, + const std::string label, + const f_t val) { if (write_graphviz) { log.printf("Node%d [label=\"%s %.16e\"]\n", node_ptr->node_id, label.c_str(), val); } } - void graphviz_edge(const mip_node_t* origin_ptr, + void graphviz_edge(logger_t& log, + const mip_node_t* origin_ptr, const mip_node_t* dest_ptr, const i_t branch_var, const i_t branch_dir, @@ -334,7 +334,6 @@ class search_tree_t { mip_node_t root; omp_mutex_t mutex; omp_atomic_t num_nodes; - logger_t log; static constexpr int write_graphviz = false; }; diff --git a/cpp/src/dual_simplex/pseudo_costs.cpp b/cpp/src/dual_simplex/pseudo_costs.cpp index 4d3db6f06b..4bd9590e16 100644 --- a/cpp/src/dual_simplex/pseudo_costs.cpp +++ b/cpp/src/dual_simplex/pseudo_costs.cpp @@ -79,7 +79,6 @@ void strong_branch_helper(i_t start, solution, iter, child_edge_norms); - // const f_t lp_solve_time = toc(lp_start_time); f_t obj = std::numeric_limits::infinity(); if (status == dual::status_t::DUAL_UNBOUNDED) { @@ -126,9 +125,7 @@ void strong_branch_helper(i_t start, } if (toc(start_time) > settings.time_limit) { break; } - pc.strong_branches_completed.lock(); const i_t completed = pc.num_strong_branches_completed++; - pc.strong_branches_completed.unlock(); if (thread_id == 0 && toc(last_log) > 10) { last_log = tic(); @@ -211,6 +208,7 @@ template void pseudo_costs_t::update_pseudo_costs(mip_node_t* node_ptr, f_t leaf_objective) { + mutex.lock(); const f_t change_in_obj = leaf_objective - node_ptr->lower_bound; const f_t frac = node_ptr->branch_dir == 0 ? node_ptr->fractional_val - std::floor(node_ptr->fractional_val) @@ -222,6 +220,7 @@ void pseudo_costs_t::update_pseudo_costs(mip_node_t* node_pt pseudo_cost_sum_up[node_ptr->branch_var] += change_in_obj / frac; pseudo_cost_num_up[node_ptr->branch_var]++; } + mutex.unlock(); } template @@ -260,10 +259,10 @@ void pseudo_costs_t::initialized(i_t& num_initialized_down, template i_t pseudo_costs_t::variable_selection(const std::vector& fractional, const std::vector& solution, - const std::vector& lower, - const std::vector& upper, - logger_t& log) const + logger_t& log) { + mutex.lock(); + const i_t num_fractional = fractional.size(); std::vector pseudo_cost_up(num_fractional); std::vector pseudo_cost_down(num_fractional); @@ -316,6 +315,8 @@ i_t pseudo_costs_t::variable_selection(const std::vector& fractio log.printf( "pc branching on %d. Value %e. Score %e\n", branch_var, solution[branch_var], score[select]); + mutex.unlock(); + return branch_var; } diff --git a/cpp/src/dual_simplex/pseudo_costs.hpp b/cpp/src/dual_simplex/pseudo_costs.hpp index d26fe8d489..5bd03e3fcd 100644 --- a/cpp/src/dual_simplex/pseudo_costs.hpp +++ b/cpp/src/dual_simplex/pseudo_costs.hpp @@ -21,8 +21,9 @@ #include #include #include +#include -#include +#include namespace cuopt::linear_programming::dual_simplex { @@ -54,9 +55,7 @@ class pseudo_costs_t { i_t variable_selection(const std::vector& fractional, const std::vector& solution, - const std::vector& lower, - const std::vector& upper, - logger_t& log) const; + logger_t& log); void update_pseudo_costs_from_strong_branching(const std::vector& fractional, const std::vector& root_soln); @@ -67,8 +66,8 @@ class pseudo_costs_t { std::vector strong_branch_down; std::vector strong_branch_up; - std::mutex strong_branches_completed; - i_t num_strong_branches_completed = 0; + omp_mutex_t mutex; + omp_atomic_t num_strong_branches_completed = 0; }; template From ca507393b99dde3d1758bdf5f0b43775b6a66050 Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 1 Oct 2025 18:02:31 +0200 Subject: [PATCH 24/33] added numerical status --- cpp/src/dual_simplex/branch_and_bound.cpp | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 3d4a206e04..970ee19d4c 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -438,6 +438,11 @@ mip_status_t branch_and_bound_t::set_final_solution(mip_solution_t::solve(mip_solution_t& solut log.log = false; status_ = mip_exploration_status_t::UNSET; - stats_.total_lp_iters = 0; - stats_.nodes_explored = 0; - stats_.nodes_unexplored = 0; - active_subtrees_ = 0; + stats_.total_lp_iters = 0; + stats_.nodes_explored = 0; + stats_.nodes_unexplored = 0; + stats_.nodes_since_last_log = 0; + stats_.last_log = 0.0; + active_subtrees_ = 0; if (guess_.size() != 0) { std::vector crushed_guess; From 7954284b826528b93a9b8c2ce54b084c5cabbacf Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 1 Oct 2025 18:19:44 +0200 Subject: [PATCH 25/33] adjusted report frequency --- cpp/src/dual_simplex/branch_and_bound.cpp | 8 +++++--- cpp/src/dual_simplex/branch_and_bound.hpp | 4 ++-- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 970ee19d4c..b247c8e844 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -696,6 +696,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* nodes_explored = (stats_.nodes_explored++); nodes_unexplored = (stats_.nodes_unexplored--); + stats_.nodes_since_last_log++; if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree->graphviz_node(settings_.log, node, "cutoff", node->lower_bound); @@ -706,10 +707,9 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* f_t now = toc(stats_.start_time); if (omp_get_thread_num() == 0) { - stats_.nodes_since_last_log++; f_t time_since_last_log = stats_.last_log == 0 ? 1.0 : toc(stats_.last_log); - if (((stats_.nodes_since_last_log > 10 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && + if (((stats_.nodes_since_last_log >= 10 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && (time_since_last_log >= 1)) || (time_since_last_log > 30) || now > settings_.time_limit) { f_t obj = compute_user_objective(original_lp_, upper_bound); @@ -725,6 +725,8 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, gap_user.c_str(), now); + + stats_.nodes_since_last_log = 0; } } @@ -789,6 +791,7 @@ void branch_and_bound_t::explore_subtree(i_t id, local_lower_bounds_[id] = lower_bound; i_t nodes_explored = stats_.nodes_explored++; i_t nodes_unexplored = stats_.nodes_unexplored--; + stats_.nodes_since_last_log++; if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound); @@ -799,7 +802,6 @@ void branch_and_bound_t::explore_subtree(i_t id, f_t now = toc(stats_.start_time); if (id == 0) { - stats_.nodes_since_last_log++; f_t time_since_last_log = stats_.last_log == 0 ? 1.0 : toc(stats_.last_log); if (((stats_.nodes_since_last_log >= 1000 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index f5d9323741..d1ed1c495d 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -149,8 +149,8 @@ class branch_and_bound_t { omp_atomic_t total_lp_iters = 0; // This should only be used by the main thread - f_t last_log = 0.0; - i_t nodes_since_last_log = 0; + f_t last_log = 0.0; + omp_atomic_t nodes_since_last_log = 0; } stats_; // Mutex for repair From 25e5fda7edb46ba5d53a81e6a1eb3c4c9559235c Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 2 Oct 2025 11:20:24 +0200 Subject: [PATCH 26/33] Arow is now shared by all threads --- cpp/src/dual_simplex/branch_and_bound.cpp | 67 +++++++++++------------ cpp/src/dual_simplex/branch_and_bound.hpp | 26 ++++----- cpp/src/dual_simplex/mip_node.hpp | 2 +- cpp/src/linear_programming/utils.cuh | 2 +- 4 files changed, 48 insertions(+), 49 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index b247c8e844..6492690360 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -550,14 +550,13 @@ branch_and_bound_t::child_selection(mip_node_t* node_ptr) } template -node_status_t branch_and_bound_t::solve_node_lp_and_update_tree( - search_tree_t& search_tree, - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log, - char symbol) +node_status_t branch_and_bound_t::solve_node(search_tree_t& search_tree, + mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + const csc_matrix_t& Arow, + f_t upper_bound, + logger_t& log, + char symbol) { f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; @@ -633,11 +632,11 @@ node_status_t branch_and_bound_t::solve_node_lp_and_update_tree( std::vector leaf_fractional; i_t leaf_num_fractional = fractional_variables(settings_, leaf_solution.x, var_types_, leaf_fractional); - f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); - search_tree.graphviz_node(log, node_ptr, "lower bound", leaf_objective); - pc_.update_pseudo_costs(node_ptr, leaf_objective); + f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); node_ptr->lower_bound = leaf_objective; + search_tree.graphviz_node(log, node_ptr, "lower bound", leaf_objective); + pc_.update_pseudo_costs(node_ptr, leaf_objective); if (leaf_num_fractional == 0) { // Found a integer feasible solution @@ -681,7 +680,7 @@ template void branch_and_bound_t::exploration_ramp_up(search_tree_t* search_tree, mip_node_t* node, lp_problem_t& leaf_problem, - csc_matrix_t& Arow, + const csc_matrix_t& Arow, i_t initial_heap_size) { if (status_ != mip_exploration_status_t::RUNNING) { return; } @@ -734,8 +733,8 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* status_ = mip_exploration_status_t::TIME_LIMIT; return; } - node_status_t node_status = solve_node_lp_and_update_tree( - *search_tree, node, leaf_problem, Arow, upper_bound, settings_.log, 'B'); + node_status_t node_status = + solve_node(*search_tree, node, leaf_problem, Arow, upper_bound, settings_.log, 'B'); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -766,7 +765,7 @@ void branch_and_bound_t::explore_subtree(i_t id, search_tree_t& search_tree, mip_node_t* start_node, lp_problem_t& leaf_problem, - csc_matrix_t& Arow) + const csc_matrix_t& Arow) { std::deque*> stack; stack.push_front(start_node); @@ -829,8 +828,8 @@ void branch_and_bound_t::explore_subtree(i_t id, break; } - node_status_t node_status = solve_node_lp_and_update_tree( - search_tree, node_ptr, leaf_problem, Arow, upper_bound, settings_.log, 'B'); + node_status_t node_status = + solve_node(search_tree, node_ptr, leaf_problem, Arow, upper_bound, settings_.log, 'B'); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -875,7 +874,7 @@ template void branch_and_bound_t::best_first_thread(i_t id, search_tree_t& search_tree, lp_problem_t& leaf_problem, - csc_matrix_t& Arow) + const csc_matrix_t& Arow) { f_t lower_bound = -inf; f_t upper_bound = inf; @@ -926,7 +925,7 @@ void branch_and_bound_t::best_first_thread(i_t id, template void branch_and_bound_t::diving_thread(lp_problem_t& leaf_problem, - csc_matrix_t& Arow) + const csc_matrix_t& Arow) { logger_t log; log.log = false; @@ -956,8 +955,8 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr continue; } - node_status_t node_status = solve_node_lp_and_update_tree( - subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); + node_status_t node_status = + solve_node(subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); if (node_status == node_status_t::TIME_LIMIT) { break; @@ -984,15 +983,10 @@ template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { logger_t log; - log.log = false; - status_ = mip_exploration_status_t::UNSET; - - stats_.total_lp_iters = 0; - stats_.nodes_explored = 0; - stats_.nodes_unexplored = 0; - stats_.nodes_since_last_log = 0; - stats_.last_log = 0.0; - active_subtrees_ = 0; + log.log = false; + status_ = mip_exploration_status_t::UNSET; + stats_.nodes_unexplored = 0; + stats_.nodes_explored = 0; if (guess_.size() != 0) { std::vector crushed_guess; @@ -1120,7 +1114,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_vstatus_, original_lp_, log); - stats_.nodes_unexplored = 2; settings_.log.printf("Exploring the B&B tree using %d threads (%d diving threads)\n", settings_.num_threads, @@ -1129,14 +1122,20 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut " | Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap " "| Time |\n"); - status_ = mip_exploration_status_t::RUNNING; + csc_matrix_t Arow(1, 1, 1); + original_lp_.A.transpose(Arow); + + stats_.nodes_explored = 0; + stats_.nodes_unexplored = 2; + stats_.nodes_since_last_log = 0; + stats_.last_log = 0.0; + active_subtrees_ = 0; + status_ = mip_exploration_status_t::RUNNING; #pragma omp parallel num_threads(settings_.num_threads) { // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; - csc_matrix_t Arow(1, 1, 1); - leaf_problem.A.transpose(Arow); #pragma omp master { diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index d1ed1c495d..d5f65d2138 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -198,7 +198,7 @@ class branch_and_bound_t { void exploration_ramp_up(search_tree_t* search_tree, mip_node_t* node, lp_problem_t& leaf_problem, - csc_matrix_t& Arow, + const csc_matrix_t& Arow, i_t initial_heap_size); // Explore the search tree using the best-first search with plunging strategy. @@ -206,27 +206,27 @@ class branch_and_bound_t { search_tree_t& search_tree, mip_node_t* start_node, lp_problem_t& leaf_problem, - csc_matrix_t& Arow); + const csc_matrix_t& Arow); // Each "main" thread pops a node from the global heap and then performs a plunge // (i.e., a shallow dive) into the subtree determined by the node. void best_first_thread(i_t id, search_tree_t& search_tree, lp_problem_t& leaf_problem, - csc_matrix_t& Arow); + const csc_matrix_t& Arow); // Each diving thread pops the first node from the dive queue and then performs // a deep dive into the subtree determined by the node. - void diving_thread(lp_problem_t& leaf_problem, csc_matrix_t& Arow); - - // Solve the LP relaxation of a leaf node. - node_status_t solve_node_lp_and_update_tree(search_tree_t& search_tree, - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log, - char symbol); + void diving_thread(lp_problem_t& leaf_problem, const csc_matrix_t& Arow); + + // Solve the LP relaxation of a leaf node and update the tree. + node_status_t solve_node(search_tree_t& search_tree, + mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + const csc_matrix_t& Arow, + f_t upper_bound, + logger_t& log, + char symbol); // Sort the children based on the Martin's criteria. std::pair*, mip_node_t*> child_selection( diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index dccfa4480b..f3476d5ba8 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -335,7 +335,7 @@ class search_tree_t { omp_mutex_t mutex; omp_atomic_t num_nodes; - static constexpr int write_graphviz = false; + static constexpr bool write_graphviz = false; }; } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/linear_programming/utils.cuh b/cpp/src/linear_programming/utils.cuh index ca77146f6d..2b2211b279 100644 --- a/cpp/src/linear_programming/utils.cuh +++ b/cpp/src/linear_programming/utils.cuh @@ -202,7 +202,7 @@ void inline compute_sum_bounds(const rmm::device_uvector& constraint_lower_ rmm::cuda_stream_view stream_view) { rmm::device_buffer d_temp_storage; - size_t bytes; + size_t bytes = 0; auto main_op = [] HD(const thrust::tuple t) { const f_t lower = thrust::get<0>(t); const f_t upper = thrust::get<1>(t); From de24ca3376d9e8cd26fb7bed9b26053e73d1c1f6 Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 2 Oct 2025 14:19:16 +0200 Subject: [PATCH 27/33] minor changes and variable renaming --- cpp/src/dual_simplex/branch_and_bound.cpp | 12 +++++++----- cpp/src/dual_simplex/branch_and_bound.hpp | 1 + cpp/src/dual_simplex/mip_node.hpp | 10 +++++----- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 6492690360..7563e2c8eb 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -540,8 +540,9 @@ branch_and_bound_t::child_selection(mip_node_t* node_ptr) const f_t up_val = std::ceil(root_relax_soln_.x[branch_var]); const f_t down_dist = branch_var_val - down_val; const f_t up_dist = up_val - branch_var_val; + constexpr f_t eps = 1e-6; - if (down_dist < up_dist) { + if (down_dist < up_dist + eps) { return std::make_pair(node_ptr->get_down_child(), node_ptr->get_up_child()); } else { @@ -782,8 +783,8 @@ void branch_and_bound_t::explore_subtree(i_t id, f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); // This is based on three assumptions: - // - The stack only contains sibling nodes, i.e., the current node and its siblings, if it - // exists + // - The stack only contains sibling nodes, i.e., the current node and it sibling (if + // applicable) // - The current node and its siblings uses the lower bound of the parent before solving the LP // relaxation // - The lower bound of the parent is lower or equal to its children @@ -938,7 +939,7 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr if (dive_queue_.size() > 0) { start_node = dive_queue_.pop(); } mutex_dive_queue_.unlock(); - if (start_node) { + if (start_node.has_value()) { if (get_upper_bound() < start_node->lower_bound) { continue; } search_tree_t subtree(std::move(start_node.value())); @@ -966,7 +967,7 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr stack.push_front(second); stack.push_front(first); - if (dive_queue_.size() < 4 * settings_.num_diving_threads) { + if (dive_queue_.size() < min_diving_queue_size_) { mutex_dive_queue_.lock(); mip_node_t* new_node = stack.back(); stack.pop_back(); @@ -1130,6 +1131,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut stats_.nodes_since_last_log = 0; stats_.last_log = 0.0; active_subtrees_ = 0; + min_diving_queue_size_ = 4 * settings_.num_diving_threads; status_ = mip_exploration_status_t::RUNNING; #pragma omp parallel num_threads(settings_.num_threads) diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index d5f65d2138..3e6fa9f445 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -176,6 +176,7 @@ class branch_and_bound_t { // Queue for storing the promising node for performing dives. omp_mutex_t mutex_dive_queue_; dive_queue_t dive_queue_; + i_t min_diving_queue_size_; // Global status omp_atomic_t status_; diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index f3476d5ba8..3f262d393d 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -280,7 +280,7 @@ class search_tree_t { void branch(mip_node_t* parent_node, const i_t branch_var, - const f_t branch_var_val, + const f_t fractional_val, const std::vector& parent_vstatus, const lp_problem_t& original_lp, logger_t& log) @@ -289,15 +289,15 @@ class search_tree_t { // down child auto down_child = std::make_unique>( - original_lp, parent_node, ++id, branch_var, 0, branch_var_val, parent_vstatus); + original_lp, parent_node, ++id, branch_var, 0, fractional_val, parent_vstatus); - graphviz_edge(log, parent_node, down_child.get(), branch_var, 0, std::floor(branch_var_val)); + graphviz_edge(log, parent_node, down_child.get(), branch_var, 0, std::floor(fractional_val)); // up child auto up_child = std::make_unique>( - original_lp, parent_node, ++id, branch_var, 1, branch_var_val, parent_vstatus); + original_lp, parent_node, ++id, branch_var, 1, fractional_val, parent_vstatus); - graphviz_edge(log, parent_node, up_child.get(), branch_var, 1, std::ceil(branch_var_val)); + graphviz_edge(log, parent_node, up_child.get(), branch_var, 1, std::ceil(fractional_val)); assert(parent_vstatus.size() == original_lp.num_cols); parent_node->add_children(std::move(down_child), From 518639dbd6aafa7fa4154a52e4117e42ed0fb74e Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 2 Oct 2025 17:16:03 +0200 Subject: [PATCH 28/33] minor change based on feedback. --- cpp/src/dual_simplex/branch_and_bound.cpp | 29 ++++++++++++++++------- cpp/src/dual_simplex/branch_and_bound.hpp | 2 +- 2 files changed, 22 insertions(+), 9 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 7563e2c8eb..13d67a2b79 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -685,7 +685,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* i_t initial_heap_size) { if (status_ != mip_exploration_status_t::RUNNING) { return; } - repair_heuristic_solutions(); + if (omp_get_thread_num() == 0) { repair_heuristic_solutions(); } f_t lower_bound = node->lower_bound; f_t upper_bound = get_upper_bound(); @@ -744,6 +744,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* } else if (node_status == node_status_t::HAS_CHILDREN) { stats_.nodes_unexplored += 2; + // If we haven't generated enough nodes to keep the threads busy, continue the ramp up phase if (stats_.nodes_unexplored < initial_heap_size) { #pragma omp task exploration_ramp_up( @@ -753,6 +754,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* exploration_ramp_up(search_tree, node->get_up_child(), leaf_problem, Arow, initial_heap_size); } else { + // We've generated enough nodes, push further nodes onto the heap mutex_heap_.lock(); heap_.push(node->get_down_child()); heap_.push(node->get_up_child()); @@ -772,7 +774,7 @@ void branch_and_bound_t::explore_subtree(i_t id, stack.push_front(start_node); while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { - repair_heuristic_solutions(); + if (omp_get_thread_num() == 0) { repair_heuristic_solutions(); } mip_node_t* node_ptr = stack.front(); stack.pop_front(); @@ -919,8 +921,12 @@ void branch_and_bound_t::best_first_thread(i_t id, // Check if it is the last thread that exited the loop and no // timeout or numerical error has happen. - if (status_ == mip_exploration_status_t::RUNNING && active_subtrees_ == 0) { - status_ = mip_exploration_status_t::COMPLETED; + if (status_ == mip_exploration_status_t::RUNNING) { + if (active_subtrees_ == 0) { + status_ = mip_exploration_status_t::COMPLETED; + } else { + local_lower_bounds_[id] = inf; + } } } @@ -967,6 +973,10 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr stack.push_front(second); stack.push_front(first); + // If the diving thread is consuming the nodes faster than the + // best first search, then we split the current subtree at the + // lowest possible point and move to the queue, so it can + // be picked by another thread. if (dive_queue_.size() < min_diving_queue_size_) { mutex_dive_queue_.lock(); mip_node_t* new_node = stack.back(); @@ -1116,9 +1126,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut original_lp_, log); - settings_.log.printf("Exploring the B&B tree using %d threads (%d diving threads)\n", - settings_.num_threads, - settings_.num_diving_threads); + settings_.log.printf( + "Exploring the B&B tree using %d best-first threads and %d diving threads (%d threads)\n", + settings_.num_bfs_threads, + settings_.num_diving_threads, + settings_.num_threads); + settings_.log.printf( " | Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap " "| Time |\n"); @@ -1129,7 +1142,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut stats_.nodes_explored = 0; stats_.nodes_unexplored = 2; stats_.nodes_since_last_log = 0; - stats_.last_log = 0.0; + stats_.last_log = tic(); active_subtrees_ = 0; min_diving_queue_size_ = 4 * settings_.num_diving_threads; status_ = mip_exploration_status_t::RUNNING; diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 3e6fa9f445..c03f46856a 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -62,7 +62,7 @@ template class dive_queue_t { private: std::vector> buffer; - static constexpr i_t max_size_ = 8192; + static constexpr i_t max_size_ = 2048; public: dive_queue_t() { buffer.reserve(max_size_); } From 4aa8549262162f3ece88b9495d999bd2764992b5 Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 2 Oct 2025 17:38:23 +0200 Subject: [PATCH 29/33] reduce the verbosity in the default log level --- cpp/src/mip/diversity/diversity_manager.cu | 6 +++--- cpp/src/mip/diversity/population.cu | 12 ++++++------ cpp/src/mip/local_search/rounding/constraint_prop.cu | 1 - .../mip/local_search/rounding/lb_constraint_prop.cu | 1 - 4 files changed, 9 insertions(+), 11 deletions(-) diff --git a/cpp/src/mip/diversity/diversity_manager.cu b/cpp/src/mip/diversity/diversity_manager.cu index 9df0966554..a9cd335a90 100644 --- a/cpp/src/mip/diversity/diversity_manager.cu +++ b/cpp/src/mip/diversity/diversity_manager.cu @@ -312,9 +312,9 @@ void diversity_manager_t::run_fj_alone(solution_t& solution) template void diversity_manager_t::run_fp_alone(solution_t& solution) { - CUOPT_LOG_INFO("Running FP alone!"); + CUOPT_LOG_DEBUG("Running FP alone!"); ls.run_fp(solution, timer, &population); - CUOPT_LOG_INFO("FP alone finished!"); + CUOPT_LOG_DEBUG("FP alone finished!"); } template @@ -543,7 +543,7 @@ void diversity_manager_t::recombine_and_ls_with_all( { raft::common::nvtx::range fun_scope("recombine_and_ls_with_all"); if (solutions.size() > 0) { - CUOPT_LOG_INFO("Running recombiners on B&B solutions with size %lu", solutions.size()); + CUOPT_LOG_DEBUG("Running recombiners on B&B solutions with size %lu", solutions.size()); // add all solutions because time limit might have been consumed and we might have exited before for (auto& sol : solutions) { cuopt_func_call(sol.test_feasibility(true)); diff --git a/cpp/src/mip/diversity/population.cu b/cpp/src/mip/diversity/population.cu index 17c88df3fc..0b506cbcb4 100644 --- a/cpp/src/mip/diversity/population.cu +++ b/cpp/src/mip/diversity/population.cu @@ -168,10 +168,10 @@ void population_t::add_external_solution(const std::vector& solut external_solution_queue_cpufj.erase(external_solution_queue_cpufj.begin() + worst_obj_idx); } - CUOPT_LOG_INFO("%s added a solution to population, solution queue size %lu with objective %g", - solution_origin_to_string(origin), - external_solution_queue.size(), - problem_ptr->get_user_obj_from_solver_obj(objective)); + CUOPT_LOG_DEBUG("%s added a solution to population, solution queue size %lu with objective %g", + solution_origin_to_string(origin), + external_solution_queue.size(), + problem_ptr->get_user_obj_from_solver_obj(objective)); if (external_solution_queue.size() >= 5) { early_exit_primal_generation = true; } } @@ -227,8 +227,8 @@ std::vector> population_t::get_external_solutions } } if (external_solution_queue.size() > 0) { - CUOPT_LOG_INFO("Consuming B&B solutions, solution queue size %lu", - external_solution_queue.size()); + CUOPT_LOG_DEBUG("Consuming B&B solutions, solution queue size %lu", + external_solution_queue.size()); external_solution_queue.clear(); } external_solution_queue_cpufj.clear(); diff --git a/cpp/src/mip/local_search/rounding/constraint_prop.cu b/cpp/src/mip/local_search/rounding/constraint_prop.cu index 4dfd1b216b..28dfa4ccde 100644 --- a/cpp/src/mip/local_search/rounding/constraint_prop.cu +++ b/cpp/src/mip/local_search/rounding/constraint_prop.cu @@ -907,7 +907,6 @@ bool constraint_prop_t::find_integer( CUOPT_LOG_DEBUG("Bounds propagation rounding: unset vars %lu", unset_integer_vars.size()); if (unset_integer_vars.size() == 0) { - CUOPT_LOG_ERROR("No integer variables provided in the bounds prop rounding"); expand_device_copy(orig_sol.assignment, sol.assignment, sol.handle_ptr->get_stream()); cuopt_func_call(orig_sol.test_variable_bounds()); return orig_sol.compute_feasibility(); diff --git a/cpp/src/mip/local_search/rounding/lb_constraint_prop.cu b/cpp/src/mip/local_search/rounding/lb_constraint_prop.cu index b9408755c8..be74611063 100644 --- a/cpp/src/mip/local_search/rounding/lb_constraint_prop.cu +++ b/cpp/src/mip/local_search/rounding/lb_constraint_prop.cu @@ -770,7 +770,6 @@ bool lb_constraint_prop_t::find_integer( { RAFT_CHECK_CUDA(problem.handle_ptr->get_stream()); if (orig_sol.problem_ptr->n_integer_vars == 0) { - CUOPT_LOG_ERROR("No integer variables provided in the bounds prop rounding"); cuopt_func_call(orig_sol.test_variable_bounds()); return orig_sol.compute_feasibility(); } From f6b1c9dc302dc1ce0f575150e7be3eb0518f639f Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 2 Oct 2025 18:49:35 +0200 Subject: [PATCH 30/33] fixed incorrect comparison --- cpp/src/dual_simplex/branch_and_bound.cpp | 2 +- cpp/src/dual_simplex/mip_node.hpp | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 13d67a2b79..3af673c8e2 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -958,7 +958,7 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr f_t upper_bound = get_upper_bound(); f_t rel_gap = user_relative_gap(original_lp_, upper_bound, node_ptr->lower_bound); - if (node_ptr->lower_bound > upper_bound || rel_gap > settings_.relative_mip_gap_tol) { + if (node_ptr->lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { continue; } diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 3f262d393d..a3a34eb81d 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -214,6 +214,7 @@ class mip_node_t { copy.branch_var_lower = branch_var_lower; copy.branch_var_upper = branch_var_upper; copy.fractional_val = fractional_val; + copy.node_id = node_id; return copy; } From 08b2191f858a75c582bdf3d32b0eca670476622d Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 2 Oct 2025 22:02:57 +0200 Subject: [PATCH 31/33] handle nodes with numerical issues --- cpp/src/dual_simplex/branch_and_bound.cpp | 10 ++++++++-- cpp/src/dual_simplex/branch_and_bound.hpp | 6 +++++- cpp/src/utilities/omp_helpers.hpp | 6 ++++-- 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 3af673c8e2..93e63d7563 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -235,7 +235,7 @@ f_t branch_and_bound_t::get_upper_bound() template f_t branch_and_bound_t::get_lower_bound() { - f_t lower_bound = inf; + f_t lower_bound = lower_bound_ceiling_.load(); mutex_heap_.lock(); if (heap_.size() > 0) { lower_bound = heap_.top()->lower_bound; } mutex_heap_.unlock(); @@ -669,10 +669,14 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& } else { search_tree.graphviz_node(log, node_ptr, "numerical", 0.0); - log.printf("Encountered LP status %d on node %d. This indicates a numerical issue.\n", + lower_bound_ceiling_.fetch_min(node_ptr->lower_bound); + log.printf("LP returned status %d on node %d. This indicates a numerical issue.\n", lp_status, node_ptr->node_id); + log.printf("The maximum lower bound is set to %+10.6e.\n", + compute_user_objective(original_lp_, lower_bound_ceiling_.load())); search_tree.update_tree(node_ptr, node_status_t::NUMERICAL); + return node_status_t::NUMERICAL; } } @@ -790,6 +794,7 @@ void branch_and_bound_t::explore_subtree(i_t id, // - The current node and its siblings uses the lower bound of the parent before solving the LP // relaxation // - The lower bound of the parent is lower or equal to its children + assert(id < local_lower_bounds_.size()); local_lower_bounds_[id] = lower_bound; i_t nodes_explored = stats_.nodes_explored++; i_t nodes_unexplored = stats_.nodes_unexplored--; @@ -1146,6 +1151,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut active_subtrees_ = 0; min_diving_queue_size_ = 4 * settings_.num_diving_threads; status_ = mip_exploration_status_t::RUNNING; + lower_bound_ceiling_ = inf; #pragma omp parallel num_threads(settings_.num_threads) { diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index c03f46856a..8a022c0c79 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -178,9 +178,13 @@ class branch_and_bound_t { dive_queue_t dive_queue_; i_t min_diving_queue_size_; - // Global status + // Global status of the solver. omp_atomic_t status_; + // In case, a best-first thread encounters a numerical issue when solving a node, + // its blocks the progression of the lower bound. + omp_atomic_t lower_bound_ceiling_; + // Set the final solution. mip_status_t set_final_solution(mip_solution_t& solution, f_t lower_bound); diff --git a/cpp/src/utilities/omp_helpers.hpp b/cpp/src/utilities/omp_helpers.hpp index f5f45efb2d..d5b91c6f91 100644 --- a/cpp/src/utilities/omp_helpers.hpp +++ b/cpp/src/utilities/omp_helpers.hpp @@ -101,8 +101,10 @@ class omp_atomic_t { T fetch_sub(T inc) { return fetch_add(-inc); } -// OpenMP only supports atomic CAS operations in v5.1 -#if _OPENMP == 202011 +// Atomic CAS are only supported in OpenMP v5.1 +// (gcc 12+ or clang 14+), however, nvcc (or the host compiler) cannot +// parse it correctly yet +#ifndef __NVCC__ T fetch_min(T other) { From 523d405c6185e511d8eb22494cf4d681ca10a16c Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 2 Oct 2025 22:11:08 +0200 Subject: [PATCH 32/33] fixed missing operation --- cpp/src/dual_simplex/branch_and_bound.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 93e63d7563..fd62504ee6 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -237,7 +237,7 @@ f_t branch_and_bound_t::get_lower_bound() { f_t lower_bound = lower_bound_ceiling_.load(); mutex_heap_.lock(); - if (heap_.size() > 0) { lower_bound = heap_.top()->lower_bound; } + if (heap_.size() > 0) { lower_bound = std::min(heap_.top()->lower_bound, lower_bound); } mutex_heap_.unlock(); for (i_t i = 0; i < local_lower_bounds_.size(); ++i) { From f31d5f3fb92cf09e97bdb5d35dcc52c283392346 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 3 Oct 2025 16:43:05 +0200 Subject: [PATCH 33/33] limiting the numerical issue reporting to the bfs thread --- cpp/src/dual_simplex/branch_and_bound.cpp | 26 +++++++++++++---------- cpp/src/dual_simplex/branch_and_bound.hpp | 4 ++-- 2 files changed, 17 insertions(+), 13 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index fd62504ee6..5870a1b46a 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -495,7 +495,7 @@ template void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, const std::vector& leaf_solution, i_t leaf_depth, - char symbol) + char thread_type) { bool send_solution = false; i_t nodes_explored = stats_.nodes_explored; @@ -509,7 +509,7 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, f_t obj = compute_user_objective(original_lp_, upper_bound_); f_t lower = compute_user_objective(original_lp_, lower_bound); settings_.log.printf("%c%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - symbol, + thread_type, nodes_explored, nodes_unexplored, obj, @@ -557,7 +557,7 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& const csc_matrix_t& Arow, f_t upper_bound, logger_t& log, - char symbol) + char thread_type) { f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; @@ -641,7 +641,7 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& if (leaf_num_fractional == 0) { // Found a integer feasible solution - add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, symbol); + add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, thread_type); search_tree.graphviz_node(log, node_ptr, "integer feasible", leaf_objective); search_tree.update_tree(node_ptr, node_status_t::INTEGER_FEASIBLE); return node_status_t::INTEGER_FEASIBLE; @@ -668,15 +668,19 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& return node_status_t::TIME_LIMIT; } else { + if (thread_type == 'B') { + lower_bound_ceiling_.fetch_min(node_ptr->lower_bound); + log.printf( + "LP returned status %d on node %d. This indicates a numerical issue. The best bound is set " + "to " + "%+10.6e.\n", + lp_status, + node_ptr->node_id, + compute_user_objective(original_lp_, lower_bound_ceiling_.load())); + } + search_tree.graphviz_node(log, node_ptr, "numerical", 0.0); - lower_bound_ceiling_.fetch_min(node_ptr->lower_bound); - log.printf("LP returned status %d on node %d. This indicates a numerical issue.\n", - lp_status, - node_ptr->node_id); - log.printf("The maximum lower bound is set to %+10.6e.\n", - compute_user_objective(original_lp_, lower_bound_ceiling_.load())); search_tree.update_tree(node_ptr, node_status_t::NUMERICAL); - return node_status_t::NUMERICAL; } } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 8a022c0c79..7b80f88fa9 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -193,7 +193,7 @@ class branch_and_bound_t { void add_feasible_solution(f_t leaf_objective, const std::vector& leaf_solution, i_t leaf_depth, - char symbol); + char thread_type); // Repairs low-quality solutions from the heuristics, if it is applicable. void repair_heuristic_solutions(); @@ -231,7 +231,7 @@ class branch_and_bound_t { const csc_matrix_t& Arow, f_t upper_bound, logger_t& log, - char symbol); + char thread_type); // Sort the children based on the Martin's criteria. std::pair*, mip_node_t*> child_selection(