diff --git a/cpp/src/mip/CMakeLists.txt b/cpp/src/mip/CMakeLists.txt index 7b4440276a..d8a6ae55b4 100644 --- a/cpp/src/mip/CMakeLists.txt +++ b/cpp/src/mip/CMakeLists.txt @@ -35,6 +35,8 @@ list(PREPEND ${CMAKE_CURRENT_SOURCE_DIR}/local_search/line_segment_search/line_segment_search.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bounds_presolve.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bounds_update_data.cu + ${CMAKE_CURRENT_SOURCE_DIR}/presolve/lb_bounds_update_data.cu + ${CMAKE_CURRENT_SOURCE_DIR}/presolve/lb_multi_probe.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/conditional_bound_strengthening.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/lb_probing_cache.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/load_balanced_bounds_presolve.cu diff --git a/cpp/src/mip/presolve/lb_bounds_update_data.cu b/cpp/src/mip/presolve/lb_bounds_update_data.cu new file mode 100644 index 0000000000..f31adc055e --- /dev/null +++ b/cpp/src/mip/presolve/lb_bounds_update_data.cu @@ -0,0 +1,112 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2022-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. + */ + +#include +#include + +#include +#include "lb_bounds_update_data.cuh" + +namespace cuopt::linear_programming::detail { + +template +lb_bounds_update_data_t::lb_bounds_update_data_t(const raft::handle_t* handle_ptr) + : bounds_changed(handle_ptr->get_stream()), + cnst_slack(0, handle_ptr->get_stream()), + vars_bnd(0, handle_ptr->get_stream()), + tmp_cnst_slack(0, handle_ptr->get_stream()), + tmp_vars_bnd(0, handle_ptr->get_stream()), + var_bounds_changed(0, handle_ptr->get_stream()), + changed_constraints(0, handle_ptr->get_stream()), + next_changed_constraints(0, handle_ptr->get_stream()), + changed_variables(0, handle_ptr->get_stream()) +{ +} + +template +void lb_bounds_update_data_t::resize(const raft::handle_t* handle_ptr, + i_t n_constraints, + i_t n_variables, + i_t num_blocks_heavy_cnst, + i_t num_blocks_heavy_vars) +{ + cnst_slack.resize(2 * n_constraints, handle_ptr->get_stream()); + tmp_cnst_slack.resize(2 * num_blocks_heavy_cnst, handle_ptr->get_stream()); + vars_bnd.resize(2 * n_variables, handle_ptr->get_stream()); + tmp_vars_bnd.resize(2 * num_blocks_heavy_vars, handle_ptr->get_stream()); + + var_bounds_changed.resize(n_variables, handle_ptr->get_stream()); + changed_constraints.resize(n_constraints, handle_ptr->get_stream()); + next_changed_constraints.resize(n_constraints, handle_ptr->get_stream()); + changed_variables.resize(n_variables, handle_ptr->get_stream()); +} + +template +typename lb_bounds_update_data_t::view_t lb_bounds_update_data_t::view() +{ + view_t v; + v.bounds_changed = bounds_changed.data(); + v.cnst_slack = make_span_2(cnst_slack); + v.vars_bnd = make_span_2(vars_bnd); + v.tmp_cnst_slack = make_span_2(cnst_slack); + v.tmp_vars_bnd = make_span_2(vars_bnd); + v.var_bounds_changed = make_span(var_bounds_changed); + v.changed_constraints = make_span(changed_constraints); + v.next_changed_constraints = make_span(next_changed_constraints); + v.changed_variables = make_span(changed_variables); + return v; +} + +template +void lb_bounds_update_data_t::init_changed_constraints(const raft::handle_t* handle_ptr) +{ + thrust::fill( + handle_ptr->get_thrust_policy(), var_bounds_changed.begin(), var_bounds_changed.end(), 0); + thrust::fill( + handle_ptr->get_thrust_policy(), changed_variables.begin(), changed_variables.end(), 1); + thrust::fill( + handle_ptr->get_thrust_policy(), changed_constraints.begin(), changed_constraints.end(), 1); + thrust::fill(handle_ptr->get_thrust_policy(), + next_changed_constraints.begin(), + next_changed_constraints.end(), + 0); +} + +template +void lb_bounds_update_data_t::prepare_for_next_iteration(const raft::handle_t* handle_ptr) +{ + std::swap(changed_constraints, next_changed_constraints); + handle_ptr->sync_stream(); + thrust::fill(handle_ptr->get_thrust_policy(), + next_changed_constraints.begin(), + next_changed_constraints.end(), + 0); + thrust::fill( + handle_ptr->get_thrust_policy(), changed_variables.begin(), changed_variables.end(), 0); + thrust::fill( + handle_ptr->get_thrust_policy(), var_bounds_changed.begin(), var_bounds_changed.end(), 0); +} + +#if MIP_INSTANTIATE_FLOAT +template class lb_bounds_update_data_t; +#endif + +#if MIP_INSTANTIATE_DOUBLE +template class lb_bounds_update_data_t; +#endif + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/lb_bounds_update_data.cuh b/cpp/src/mip/presolve/lb_bounds_update_data.cuh new file mode 100644 index 0000000000..749f2dd625 --- /dev/null +++ b/cpp/src/mip/presolve/lb_bounds_update_data.cuh @@ -0,0 +1,64 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2022-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 + +#include +#include + +#include + +namespace cuopt::linear_programming::detail { + +template +struct lb_bounds_update_data_t { + rmm::device_scalar bounds_changed; + rmm::device_uvector cnst_slack; + rmm::device_uvector vars_bnd; + rmm::device_uvector tmp_cnst_slack; + rmm::device_uvector tmp_vars_bnd; + rmm::device_uvector var_bounds_changed; + rmm::device_uvector changed_constraints; + rmm::device_uvector next_changed_constraints; + rmm::device_uvector changed_variables; + + struct view_t { + using f_t2 = typename type_2::type; + i_t* bounds_changed; + raft::device_span cnst_slack; + raft::device_span vars_bnd; + raft::device_span tmp_cnst_slack; + raft::device_span tmp_vars_bnd; + raft::device_span var_bounds_changed; + raft::device_span changed_constraints; + raft::device_span next_changed_constraints; + raft::device_span changed_variables; + }; + + lb_bounds_update_data_t(const raft::handle_t* handle); + // void resize(load_balanced_bounds_presolve_t& prs); + void resize(const raft::handle_t* handle_ptr, + i_t n_constraints, + i_t n_variables, + i_t num_blocks_heavy_cnst, + i_t num_blocks_heavy_vars); + void init_changed_constraints(const raft::handle_t* handle_ptr); + void prepare_for_next_iteration(const raft::handle_t* handle_ptr); + view_t view(); +}; + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/lb_multi_probe.cu b/cpp/src/mip/presolve/lb_multi_probe.cu new file mode 100644 index 0000000000..aa3e067beb --- /dev/null +++ b/cpp/src/mip/presolve/lb_multi_probe.cu @@ -0,0 +1,487 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2022-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. + */ + +#include +#include "lb_multi_probe.cuh" +#include "lb_multi_probe_helpers.cuh" +#include "load_balanced_bounds_presolve_helpers.cuh" + +namespace cuopt::linear_programming::detail { + +template +lb_multi_probe_t::lb_multi_probe_t(const load_balanced_problem_t& problem_, + mip_solver_context_t& context_, + settings_t in_settings, + i_t max_stream_count_) + : context(context_), + pb(&problem_), + streams(max_stream_count_), + upd_0(problem_.handle_ptr), + upd_1(problem_.handle_ptr), + heavy_cnst_block_segments(0, problem_.handle_ptr->get_stream()), + heavy_vars_block_segments(0, problem_.handle_ptr->get_stream()), + heavy_cnst_vertex_ids(0, problem_.handle_ptr->get_stream()), + heavy_vars_vertex_ids(0, problem_.handle_ptr->get_stream()), + heavy_cnst_pseudo_block_ids(0, problem_.handle_ptr->get_stream()), + heavy_vars_pseudo_block_ids(0, problem_.handle_ptr->get_stream()), + num_blocks_heavy_cnst(0), + num_blocks_heavy_vars(0), + warp_cnst_offsets(0, problem_.handle_ptr->get_stream()), + warp_cnst_id_offsets(0, problem_.handle_ptr->get_stream()), + warp_vars_offsets(0, problem_.handle_ptr->get_stream()), + warp_vars_id_offsets(0, problem_.handle_ptr->get_stream()), + calc_slack_exec(nullptr), + calc_slack_erase_inf_cnst_exec(nullptr), + upd_bnd_exec(nullptr), + calc_slack_erase_inf_cnst_graph_created(false), + calc_slack_graph_created(false), + upd_bnd_graph_created(false), + settings(in_settings) +{ + setup(problem_); +} + +template +lb_multi_probe_t::~lb_multi_probe_t() +{ + if (calc_slack_erase_inf_cnst_graph_created) { + cudaGraphExecDestroy(calc_slack_erase_inf_cnst_exec); + } + if (calc_slack_graph_created) { cudaGraphExecDestroy(calc_slack_exec); } + if (upd_bnd_graph_created) { cudaGraphExecDestroy(upd_bnd_exec); } +} + +template +void lb_multi_probe_t::setup(const load_balanced_problem_t& problem) +{ + auto handle_ptr = pb->handle_ptr; + auto stream = handle_ptr->get_stream(); + + num_blocks_heavy_cnst = create_heavy_item_block_segments(stream, + heavy_cnst_vertex_ids, + heavy_cnst_pseudo_block_ids, + heavy_cnst_block_segments, + heavy_degree_cutoff, + problem.cnst_bin_offsets, + problem.offsets); + + num_blocks_heavy_vars = create_heavy_item_block_segments(stream, + heavy_vars_vertex_ids, + heavy_vars_pseudo_block_ids, + heavy_vars_block_segments, + heavy_degree_cutoff, + problem.vars_bin_offsets, + problem.reverse_offsets); + upd_0.resize( + handle_ptr, pb->n_constraints, pb->n_variables, num_blocks_heavy_cnst, num_blocks_heavy_vars); + upd_1.resize( + handle_ptr, pb->n_constraints, pb->n_variables, num_blocks_heavy_cnst, num_blocks_heavy_vars); + + std::tie(is_cnst_sub_warp_single_bin, cnst_sub_warp_count) = + sub_warp_meta(stream, warp_cnst_offsets, warp_cnst_id_offsets, pb->cnst_bin_offsets, 4); + + std::tie(is_vars_sub_warp_single_bin, vars_sub_warp_count) = + sub_warp_meta(stream, warp_vars_offsets, warp_vars_id_offsets, pb->vars_bin_offsets, 4); + + stream.synchronize(); + streams.sync_all_issued(); + + if (!calc_slack_erase_inf_cnst_graph_created) { + bool erase_inf_cnst = true; + calc_slack_erase_inf_cnst_graph_created = build_graph( + streams, + handle_ptr, + calc_slack_erase_inf_cnst_graph, + calc_slack_erase_inf_cnst_exec, + [erase_inf_cnst, this]() { this->calculate_cnst_slack_graph(erase_inf_cnst, true); }, + [erase_inf_cnst, this]() { this->calculate_cnst_slack_graph(erase_inf_cnst); }); + } + + if (!calc_slack_graph_created) { + bool erase_inf_cnst = false; + calc_slack_graph_created = build_graph( + streams, + handle_ptr, + calc_slack_graph, + calc_slack_exec, + [erase_inf_cnst, this]() { this->calculate_cnst_slack_graph(erase_inf_cnst, true); }, + [erase_inf_cnst, this]() { this->calculate_cnst_slack_graph(erase_inf_cnst); }); + } + + if (!upd_bnd_graph_created) { + upd_bnd_graph_created = build_graph( + streams, + handle_ptr, + upd_bnd_graph, + upd_bnd_exec, + [this]() { this->calculate_bounds_update_graph(true); }, + [this]() { this->calculate_bounds_update_graph(); }); + } +} + +template +typename lb_multi_probe_t::activity_view_t lb_multi_probe_t::get_activity_view( + const load_balanced_problem_t& pb) +{ + lb_multi_probe_t::activity_view_t v; + v.cnst_reorg_ids = make_span(pb.cnst_reorg_ids); + v.coeff = make_span(pb.coefficients); + v.vars = make_span(pb.variables); + v.offsets = make_span(pb.offsets); + v.cnst_bnd = make_span_2(pb.cnst_bounds_data); + v.nnz = pb.nnz; + v.tolerances = pb.tolerances; + return v; +} + +template +typename lb_multi_probe_t::bounds_update_view_t +lb_multi_probe_t::get_bounds_update_view(const load_balanced_problem_t& pb) +{ + lb_multi_probe_t::bounds_update_view_t v; + v.vars_reorg_ids = make_span(pb.vars_reorg_ids); + v.coeff = make_span(pb.reverse_coefficients); + v.cnst = make_span(pb.reverse_constraints); + v.offsets = make_span(pb.reverse_offsets); + v.vars_types = make_span(pb.vars_types); + v.nnz = pb.nnz; + v.tolerances = pb.tolerances; + return v; +} + +template +void lb_multi_probe_t::calculate_cnst_slack(const raft::handle_t* handle_ptr) +{ + cudaGraphLaunch(calc_slack_erase_inf_cnst_exec, handle_ptr->get_stream()); +} + +template +void lb_multi_probe_t::calculate_cnst_slack_graph(bool erase_inf_cnst, bool dry_run) +{ + using f_t2 = typename type_2::type; + + auto activity_view = get_activity_view(*pb); + auto upd_0_v = upd_0.view(); + auto upd_1_v = upd_1.view(); + + calc_activity_heavy_cnst(streams, + activity_view, + upd_0_v, + upd_1_v, + heavy_cnst_vertex_ids, + heavy_cnst_pseudo_block_ids, + heavy_cnst_block_segments, + pb->cnst_bin_offsets, + heavy_degree_cutoff, + num_blocks_heavy_cnst, + erase_inf_cnst, + dry_run); + calc_activity_per_block(streams, + activity_view, + upd_0_v, + upd_1_v, + pb->cnst_bin_offsets, + heavy_degree_cutoff, + erase_inf_cnst, + dry_run); + calc_activity_sub_warp(streams, + activity_view, + upd_0_v, + upd_1_v, + is_cnst_sub_warp_single_bin, + cnst_sub_warp_count, + warp_cnst_offsets, + warp_cnst_id_offsets, + pb->cnst_bin_offsets, + erase_inf_cnst, + dry_run); +} + +template +void lb_multi_probe_t::calculate_bounds_update_graph(bool dry_run) +{ + using f_t2 = typename type_2::type; + + auto bounds_update_view = get_bounds_update_view(*pb); + auto upd_0_v = upd_0.view(); + auto upd_1_v = upd_1.view(); + + upd_bounds_heavy_vars(streams, + bounds_update_view, + upd_0_v, + upd_1_v, + heavy_vars_vertex_ids, + heavy_vars_pseudo_block_ids, + heavy_vars_block_segments, + pb->vars_bin_offsets, + heavy_degree_cutoff, + num_blocks_heavy_vars, + dry_run); + upd_bounds_per_block(streams, + bounds_update_view, + upd_0_v, + upd_1_v, + pb->vars_bin_offsets, + heavy_degree_cutoff, + dry_run); + upd_bounds_sub_warp(streams, + bounds_update_view, + upd_0_v, + upd_1_v, + is_vars_sub_warp_single_bin, + vars_sub_warp_count, + warp_vars_offsets, + warp_vars_id_offsets, + pb->vars_bin_offsets, + dry_run); +} + +template +bool lb_multi_probe_t::calculate_bounds_update(const raft::handle_t* handle_ptr) +{ + constexpr i_t zero = 0; + upd_0.bounds_changed.set_value_async(zero, handle_ptr->get_stream()); + upd_1.bounds_changed.set_value_async(zero, handle_ptr->get_stream()); + cudaGraphLaunch(upd_bnd_exec, handle_ptr->get_stream()); + i_t h_bounds_changed_0 = upd_0.bounds_changed.value(handle_ptr->get_stream()); + i_t h_bounds_changed_1 = upd_1.bounds_changed.value(handle_ptr->get_stream()); + skip_0 = (h_bounds_changed_0 == zero); + skip_1 = (h_bounds_changed_1 == zero); + return !(skip_0 && skip_1); +} + +template +termination_criterion_t lb_multi_probe_t::bound_update_loop( + const raft::handle_t* handle_ptr, timer_t timer) +{ + termination_criterion_t criteria = termination_criterion_t::ITERATION_LIMIT; + i_t iter_0 = 0; + i_t iter_1 = 0; + if (init_changed_constraints) { + // all changed constraints are 1, next are zero + upd_0.init_changed_constraints(handle_ptr); + upd_1.init_changed_constraints(handle_ptr); + } else { + // reset for the next calls on the same object + init_changed_constraints = true; + } + // settings.iteration_limit = 1; + for (i_t iter = 0; iter < settings.iteration_limit; ++iter) { + if (timer.check_time_limit()) { + criteria = termination_criterion_t::TIME_LIMIT; + break; + } + // calculate activity for both probes + calculate_cnst_slack(handle_ptr); + if (!calculate_bounds_update(handle_ptr)) { + if (iter == 0) { + criteria = termination_criterion_t::NO_UPDATE; + } else { + criteria = termination_criterion_t::CONVERGENCE; + } + break; + } + // next_changed are updated, fill current changed with zero and swap + // swap next and current changed constraints + if (!skip_0) { upd_0.prepare_for_next_iteration(handle_ptr); } + if (!skip_1) { upd_1.prepare_for_next_iteration(handle_ptr); } + iter_0 += !skip_0; + iter_1 += !skip_1; + } + + return criteria; +} + +template +void lb_multi_probe_t::copy_problem_into_probing_buffers(const raft::handle_t* handle_ptr) +{ + cuopt_assert(upd_0.vars_bnd.size() == pb->variable_bounds.size(), + "size of variable bounds mismatch"); + raft::copy(upd_0.vars_bnd.data(), + pb->variable_bounds.data(), + upd_0.vars_bnd.size(), + handle_ptr->get_stream()); + + cuopt_assert(upd_1.vars_bnd.size() == pb->variable_bounds.size(), + "size of variable bounds mismatch"); + raft::copy(upd_1.vars_bnd.data(), + pb->variable_bounds.data(), + upd_1.vars_bnd.size(), + handle_ptr->get_stream()); +} + +template +void lb_multi_probe_t::set_interval_bounds( + const std::tuple, std::pair>& var_interval_vals, + const raft::handle_t* handle_ptr) +{ + // TODO : upd.vars_bnd_changed + using f_t2 = typename type_2::type; + const i_t& probe_var = std::get<0>(var_interval_vals); + const std::pair& probe_vals_0 = std::get<1>(var_interval_vals); + const std::pair& probe_vals_1 = std::get<2>(var_interval_vals); + run_device_lambda(handle_ptr->get_stream(), + [probe_var = probe_var, + lb_0 = probe_vals_0.first, + ub_0 = probe_vals_0.second, + lb_1 = probe_vals_1.first, + ub_1 = probe_vals_1.second, + upd_0_v = upd_0.view(), + upd_1_v = upd_1.view()] __device__() { + upd_0_v.vars_bnd[probe_var] = f_t2{lb_0, ub_0}; + upd_1_v.vars_bnd[probe_var] = f_t2{lb_1, ub_1}; + }); + // init changed constraints + auto orig_pb = pb->pb; + i_t var_offset_begin = orig_pb->reverse_offsets.element(probe_var, handle_ptr->get_stream()); + i_t var_offset_end = orig_pb->reverse_offsets.element(probe_var + 1, handle_ptr->get_stream()); + thrust::fill(handle_ptr->get_thrust_policy(), + upd_0.changed_constraints.begin(), + upd_0.changed_constraints.end(), + 0); + thrust::fill(handle_ptr->get_thrust_policy(), + upd_1.changed_constraints.begin(), + upd_1.changed_constraints.end(), + 0); + thrust::fill(handle_ptr->get_thrust_policy(), + upd_0.next_changed_constraints.begin(), + upd_0.next_changed_constraints.end(), + 0); + thrust::fill(handle_ptr->get_thrust_policy(), + upd_1.next_changed_constraints.begin(), + upd_1.next_changed_constraints.end(), + 0); + // set changed constraints from the vars + thrust::for_each(handle_ptr->get_thrust_policy(), + orig_pb->reverse_constraints.begin() + var_offset_begin, + orig_pb->reverse_constraints.begin() + var_offset_end, + [upd_0_v = upd_0.view(), upd_1_v = upd_1.view()] __device__(auto i) { + upd_0_v.changed_constraints[i] = 1; + upd_1_v.changed_constraints[i] = 1; + }); + init_changed_constraints = false; + handle_ptr->sync_stream(); + RAFT_CHECK_CUDA(handle_ptr->get_stream()); +} + +template +termination_criterion_t lb_multi_probe_t::solve_for_interval( + const std::tuple, std::pair>& var_interval_vals, + const raft::handle_t* handle_ptr) +{ + timer_t timer(settings.time_limit); + + copy_problem_into_probing_buffers(handle_ptr); + set_interval_bounds(var_interval_vals, handle_ptr); + + return bound_update_loop(handle_ptr, timer); +} + +template +void lb_multi_probe_t::set_updated_bounds(const raft::handle_t* handle_ptr, + raft::device_span output_bounds, + i_t select_update) +{ + auto& bnds = select_update ? upd_1.vars_bnd : upd_0.vars_bnd; + + cuopt_assert(bnds.size() == output_bounds.size(), "size of variable bound mismatch"); + raft::copy(output_bounds.data(), bnds.data(), bnds.size(), handle_ptr->get_stream()); +} + +template +void lb_multi_probe_t::set_updated_bounds(const raft::handle_t* handle_ptr, + raft::device_span output_lb, + raft::device_span output_ub, + i_t select_update) +{ + auto& bnds = select_update ? upd_1.vars_bnd : upd_0.vars_bnd; + + auto bnd_span = make_span_2(bnds); + cuopt_assert(bnd_span.size() == output_lb.size(), "size of variable lower bound mismatch"); + cuopt_assert(bnd_span.size() == output_ub.size(), "size of variable upper bound mismatch"); + thrust::for_each(handle_ptr->get_thrust_policy(), + thrust::make_counting_iterator(0), + thrust::make_counting_iterator(bnd_span.size()), + [output_lb, output_ub, bnd_span] __device__(auto idx) { + auto bnd = bnd_span[idx]; + output_lb[idx] = bnd.x; + output_ub[idx] = bnd.y; + }); +} + +template +void lb_multi_probe_t::set_updated_bounds(load_balanced_problem_t& problem, + i_t select_update, + const raft::handle_t* handle_ptr) +{ + set_updated_bounds(handle_ptr, make_span(problem.variable_bounds), select_update); +} + +template +void lb_multi_probe_t::set_bounds( + const std::tuple, std::vector, std::vector>& var_probe_vals, + const raft::handle_t* handle_ptr) +{ + const std::vector& probe_vars = std::get<0>(var_probe_vals); + const std::vector& probe_vals_0 = std::get<1>(var_probe_vals); + const std::vector& probe_vals_1 = std::get<2>(var_probe_vals); + auto d_vars = device_copy(probe_vars, handle_ptr->get_stream()); + auto d_vals_0 = device_copy(probe_vals_0, handle_ptr->get_stream()); + auto d_vals_1 = device_copy(probe_vals_1, handle_ptr->get_stream()); + + auto upd_0_v = upd_0.view(); + auto upd_1_v = upd_1.view(); + auto z_iter = thrust::make_zip_iterator( + thrust::make_tuple(d_vars.begin(), d_vals_0.begin(), d_vals_1.begin())); + thrust::for_each( + handle_ptr->get_thrust_policy(), + z_iter, + z_iter + d_vars.size(), + [upd_0_v, upd_1_v] __device__(auto t) { + using f_t2 = typename type_2::type; + upd_0_v.vars_bnd[thrust::get<0>(t)] = f_t2{thrust::get<1>(t), thrust::get<1>(t)}; + upd_1_v.vars_bnd[thrust::get<0>(t)] = f_t2{thrust::get<2>(t), thrust::get<2>(t)}; + // upd_0_v.ub[thrust::get<0>(t)] = thrust::get<1>(t); + // upd_1_v.lb[thrust::get<0>(t)] = thrust::get<2>(t); + // upd_1_v.ub[thrust::get<0>(t)] = thrust::get<2>(t); + }); + handle_ptr->sync_stream(); + RAFT_CHECK_CUDA(handle_ptr->get_stream()); +} + +template +termination_criterion_t lb_multi_probe_t::solve( + const std::tuple, std::vector, std::vector>& var_probe_vals) +{ + timer_t timer(settings.time_limit); + auto& handle_ptr = pb->handle_ptr; + + copy_problem_into_probing_buffers(handle_ptr); + set_bounds(var_probe_vals, handle_ptr); + + return bound_update_loop(handle_ptr, timer); +} + +#if MIP_INSTANTIATE_FLOAT +template class lb_multi_probe_t; +#endif + +#if MIP_INSTANTIATE_DOUBLE +template class lb_multi_probe_t; +#endif + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/lb_multi_probe.cuh b/cpp/src/mip/presolve/lb_multi_probe.cuh new file mode 100644 index 0000000000..4823dcdda3 --- /dev/null +++ b/cpp/src/mip/presolve/lb_multi_probe.cuh @@ -0,0 +1,159 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2022-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 + +#include +#include +#include +#include + +#include + +#include "lb_bounds_update_data.cuh" +#include "load_balanced_bounds_presolve.cuh" +#include "utils.cuh" + +namespace cuopt::linear_programming::detail { + +template +class lb_multi_probe_t { + public: + static constexpr i_t heavy_degree_cutoff = 16 * 1024; + struct settings_t { + f_t time_limit{60.0}; + i_t iteration_limit{std::numeric_limits::max()}; + }; + + struct activity_view_t { + using f_t2 = typename type_2::type; + raft::device_span cnst_reorg_ids; + raft::device_span coeff; + raft::device_span vars; + raft::device_span offsets; + raft::device_span cnst_bnd; // new indexing + i_t nnz; + typename mip_solver_settings_t::tolerances_t tolerances; + }; + + struct bounds_update_view_t { + using f_t2 = typename type_2::type; + raft::device_span vars_reorg_ids; + raft::device_span coeff; + raft::device_span cnst; + raft::device_span offsets; + raft::device_span vars_types; // new indexing + i_t nnz; + typename mip_solver_settings_t::tolerances_t tolerances; + }; + ~lb_multi_probe_t(); + lb_multi_probe_t(lb_multi_probe_t&&) = default; + + lb_multi_probe_t(const load_balanced_problem_t& problem_, + mip_solver_context_t& context_, + settings_t settings = settings_t{}, + i_t max_stream_count_ = 32); + void setup(const load_balanced_problem_t& problem); + + activity_view_t get_activity_view(const load_balanced_problem_t& pb); + bounds_update_view_t get_bounds_update_view(const load_balanced_problem_t& pb); + + void calculate_cnst_slack(const raft::handle_t* handle_ptr); + bool calculate_bounds_update(const raft::handle_t* handle_ptr); + void calculate_bounds_update_graph(bool dry_run = false); + void calculate_cnst_slack_graph(bool erase_inf_cnst, bool dry_run = false); + termination_criterion_t solve( + const std::tuple, std::vector, std::vector>& var_probe_vals); + + termination_criterion_t solve_for_interval( + const std::tuple, std::pair>& var_interval_vals, + const raft::handle_t* handle_ptr); + + void set_updated_bounds(load_balanced_problem_t& problem, + i_t select_update, + const raft::handle_t* handle_ptr); + void set_updated_bounds(const raft::handle_t* handle_ptr, + raft::device_span output_bounds, + i_t select_update); + void set_updated_bounds(const raft::handle_t* handle_ptr, + raft::device_span output_lb, + raft::device_span output_ub, + i_t select_update); + termination_criterion_t bound_update_loop(const raft::handle_t* handle_ptr, timer_t timer); + void set_interval_bounds( + const std::tuple, std::pair>& var_interval_vals, + const raft::handle_t* handle_ptr); + void set_bounds( + const std::tuple, std::vector, std::vector>& var_probe_vals, + const raft::handle_t* handle_ptr); + // void constraint_stats(problem_t& pb, const raft::handle_t* handle_ptr); + void copy_problem_into_probing_buffers(const raft::handle_t* handle_ptr); + + mip_solver_context_t& context; + + const load_balanced_problem_t* pb; + + managed_stream_pool streams; + lb_bounds_update_data_t upd_0; + lb_bounds_update_data_t upd_1; + + // Number of blocks for heavy ids + rmm::device_uvector heavy_cnst_block_segments; + rmm::device_uvector heavy_vars_block_segments; + rmm::device_uvector heavy_cnst_vertex_ids; + rmm::device_uvector heavy_vars_vertex_ids; + rmm::device_uvector heavy_cnst_pseudo_block_ids; + rmm::device_uvector heavy_vars_pseudo_block_ids; + + i_t num_blocks_heavy_cnst; + i_t num_blocks_heavy_vars; + + // sub warp meta data + bool is_cnst_sub_warp_single_bin; + i_t cnst_sub_warp_count; + rmm::device_uvector warp_cnst_offsets; + rmm::device_uvector warp_cnst_id_offsets; + + bool is_vars_sub_warp_single_bin; + i_t vars_sub_warp_count; + rmm::device_uvector warp_vars_offsets; + rmm::device_uvector warp_vars_id_offsets; + + // graphs + bool calc_slack_erase_inf_cnst_graph_created; + bool calc_slack_graph_created; + bool upd_bnd_graph_created; + + cudaGraphExec_t calc_slack_erase_inf_cnst_exec; + cudaGraph_t calc_slack_erase_inf_cnst_graph; + cudaGraphExec_t calc_slack_exec; + cudaGraph_t calc_slack_graph; + cudaGraphExec_t upd_bnd_exec; + cudaGraph_t upd_bnd_graph; + + bool skip_0; + bool skip_1; + settings_t settings; + bool compute_stats = true; + bool init_changed_constraints = true; + i_t infeas_constraints_count_0 = 0; + i_t redund_constraints_count_0 = 0; + i_t infeas_constraints_count_1 = 0; + i_t redund_constraints_count_1 = 0; +}; + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/lb_multi_probe_helpers.cuh b/cpp/src/mip/presolve/lb_multi_probe_helpers.cuh new file mode 100644 index 0000000000..8108a85640 --- /dev/null +++ b/cpp/src/mip/presolve/lb_multi_probe_helpers.cuh @@ -0,0 +1,558 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2022-2025, NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: LicenseRef-NvidiaProprietary + * + * NVIDIA CORPORATION, its affiliates and licensors retain all intellectual + * property and proprietary rights in and to this material, related + * documentation and any modifications thereto. Any use, reproduction, + * disclosure or distribution of this material and related documentation + * without an express license agreement from NVIDIA CORPORATION or + * its affiliates is strictly prohibited. + */ + +#pragma once + +#include "load_balanced_bounds_next_constraint_kernels.cuh" +#include "load_balanced_bounds_presolve_helpers.cuh" +#include "load_balanced_bounds_presolve_kernels.cuh" +#include "load_balanced_partition_helpers.cuh" + +#include +#include +#include +#include +#include +#include +#include + +namespace cuopt::linear_programming::detail { + +/// CALCULATE ACTIVITY + +template +void calc_activity_heavy_cnst(stream_pool_t& streams, + activity_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + const rmm::device_uvector& heavy_cnst_vertex_ids, + const rmm::device_uvector& heavy_cnst_pseudo_block_ids, + const rmm::device_uvector& heavy_cnst_block_segments, + const std::vector& cnst_bin_offsets, + i_t heavy_degree_cutoff, + i_t num_blocks_heavy_cnst, + bool erase_inf_cnst, + bool dry_run = false) +{ + if (num_blocks_heavy_cnst != 0) { + auto heavy_cnst_stream = streams.get_stream(); + // TODO : Check heavy_cnst_block_segments size for profiling + if (!dry_run) { + auto heavy_cnst_beg_id = get_id_offset(cnst_bin_offsets, heavy_degree_cutoff); + lb_calc_act_heavy_kernel + <<>>( + heavy_cnst_beg_id, + make_span(heavy_cnst_vertex_ids), + make_span(heavy_cnst_pseudo_block_ids), + heavy_degree_cutoff, + view, + upd_0, + upd_1); + auto num_heavy_cnst = cnst_bin_offsets.back() - heavy_cnst_beg_id; + if (erase_inf_cnst) { + finalize_calc_act_kernel + <<>>( + heavy_cnst_beg_id, make_span(heavy_cnst_block_segments), view, upd_0, upd_1); + } else { + finalize_calc_act_kernel + <<>>( + heavy_cnst_beg_id, make_span(heavy_cnst_block_segments), view, upd_0, upd_1); + } + } + } +} + +template +void calc_activity_per_block(stream_pool_t& streams, + activity_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + const std::vector& cnst_bin_offsets, + i_t degree_beg, + i_t degree_end, + bool erase_inf_cnst, + bool dry_run) +{ + static_assert(block_dim <= 1024, "Cannot launch kernel with more than 1024 threads"); + + auto [cnst_id_beg, cnst_id_end] = get_id_range(cnst_bin_offsets, degree_beg, degree_end); + + auto block_count = cnst_id_end - cnst_id_beg; + if (block_count > 0) { + auto block_stream = streams.get_stream(); + if (!dry_run) { + if (erase_inf_cnst) { + lb_calc_act_block_kernel + <<>>(cnst_id_beg, view, upd_0, upd_1); + } else { + lb_calc_act_block_kernel + <<>>(cnst_id_beg, view, upd_0, upd_1); + } + } + } +} + +template +void calc_activity_per_block(stream_pool_t& streams, + activity_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + const std::vector& cnst_bin_offsets, + i_t heavy_degree_cutoff, + bool erase_inf_cnst, + bool dry_run = false) +{ + if (view.nnz < 10000) { + calc_activity_per_block( + streams, view, upd_0, upd_1, cnst_bin_offsets, 32, 32, erase_inf_cnst, dry_run); + calc_activity_per_block( + streams, view, upd_0, upd_1, cnst_bin_offsets, 64, 64, erase_inf_cnst, dry_run); + calc_activity_per_block( + streams, view, upd_0, upd_1, cnst_bin_offsets, 128, 128, erase_inf_cnst, dry_run); + calc_activity_per_block( + streams, view, upd_0, upd_1, cnst_bin_offsets, 256, 256, erase_inf_cnst, dry_run); + } else { + //[1024, heavy_degree_cutoff/2] -> 1024 block size + calc_activity_per_block(streams, + view, + upd_0, + upd_1, + cnst_bin_offsets, + 1024, + heavy_degree_cutoff / 2, + erase_inf_cnst, + dry_run); + //[512, 512] -> 128 block size + calc_activity_per_block( + streams, view, upd_0, upd_1, cnst_bin_offsets, 128, 512, erase_inf_cnst, dry_run); + } +} + +template +void calc_activity_sub_warp(stream_pool_t& streams, + activity_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + i_t degree_beg, + i_t degree_end, + const std::vector& cnst_bin_offsets, + bool erase_inf_cnst, + bool dry_run) +{ + constexpr i_t block_dim = 32; + auto cnst_per_block = block_dim / threads_per_constraint; + auto [cnst_id_beg, cnst_id_end] = get_id_range(cnst_bin_offsets, degree_beg, degree_end); + + auto block_count = raft::ceildiv(cnst_id_end - cnst_id_beg, cnst_per_block); + if (block_count != 0) { + auto sub_warp_thread = streams.get_stream(); + if (!dry_run) { + if (erase_inf_cnst) { + lb_calc_act_sub_warp_kernel + <<>>( + cnst_id_beg, cnst_id_end, view, upd_0, upd_1); + } else { + lb_calc_act_sub_warp_kernel + <<>>( + cnst_id_beg, cnst_id_end, view, upd_0, upd_1); + } + } + } +} + +template +void calc_activity_sub_warp(stream_pool_t& streams, + activity_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + i_t degree, + const std::vector& cnst_bin_offsets, + bool erase_inf_cnst, + bool dry_run) +{ + calc_activity_sub_warp( + streams, view, upd_0, upd_1, degree, degree, cnst_bin_offsets, erase_inf_cnst, dry_run); +} + +template +void calc_activity_sub_warp(stream_pool_t& streams, + activity_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + i_t cnst_sub_warp_count, + rmm::device_uvector& warp_cnst_offsets, + rmm::device_uvector& warp_cnst_id_offsets, + bool erase_inf_cnst, + bool dry_run) +{ + constexpr i_t block_dim = 256; + + auto block_count = raft::ceildiv(cnst_sub_warp_count * 32, block_dim); + if (block_count != 0) { + auto sub_warp_stream = streams.get_stream(); + if (!dry_run) { + if (erase_inf_cnst) { + lb_calc_act_sub_warp_kernel + <<>>( + view, upd_0, upd_1, make_span(warp_cnst_offsets), make_span(warp_cnst_id_offsets)); + } else { + lb_calc_act_sub_warp_kernel + <<>>( + view, upd_0, upd_1, make_span(warp_cnst_offsets), make_span(warp_cnst_id_offsets)); + } + } + } +} + +template +void calc_activity_sub_warp(stream_pool_t& streams, + activity_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + bool is_cnst_sub_warp_single_bin, + i_t cnst_sub_warp_count, + rmm::device_uvector& warp_cnst_offsets, + rmm::device_uvector& warp_cnst_id_offsets, + const std::vector& cnst_bin_offsets, + bool erase_inf_cnst, + bool dry_run = false) +{ + if (view.nnz < 10000) { + calc_activity_sub_warp( + streams, view, upd_0, upd_1, 16, cnst_bin_offsets, erase_inf_cnst, dry_run); + calc_activity_sub_warp( + streams, view, upd_0, upd_1, 8, cnst_bin_offsets, erase_inf_cnst, dry_run); + calc_activity_sub_warp( + streams, view, upd_0, upd_1, 4, cnst_bin_offsets, erase_inf_cnst, dry_run); + calc_activity_sub_warp( + streams, view, upd_0, upd_1, 2, cnst_bin_offsets, erase_inf_cnst, dry_run); + calc_activity_sub_warp( + streams, view, upd_0, upd_1, 1, cnst_bin_offsets, erase_inf_cnst, dry_run); + } else { + if (is_cnst_sub_warp_single_bin) { + calc_activity_sub_warp( + streams, view, upd_0, upd_1, 64, cnst_bin_offsets, erase_inf_cnst, dry_run); + calc_activity_sub_warp( + streams, view, upd_0, upd_1, 32, cnst_bin_offsets, erase_inf_cnst, dry_run); + calc_activity_sub_warp( + streams, view, upd_0, upd_1, 16, cnst_bin_offsets, erase_inf_cnst, dry_run); + calc_activity_sub_warp( + streams, view, upd_0, upd_1, 8, cnst_bin_offsets, erase_inf_cnst, dry_run); + calc_activity_sub_warp( + streams, view, upd_0, upd_1, 1, 4, cnst_bin_offsets, erase_inf_cnst, dry_run); + } else { + calc_activity_sub_warp(streams, + view, + upd_0, + upd_1, + cnst_sub_warp_count, + warp_cnst_offsets, + warp_cnst_id_offsets, + erase_inf_cnst, + dry_run); + } + } +} + +// bounds + +template +void upd_bounds_heavy_vars(stream_pool_t& streams, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + const rmm::device_uvector& heavy_vars_vertex_ids, + const rmm::device_uvector& heavy_vars_pseudo_block_ids, + const rmm::device_uvector& heavy_vars_block_segments, + const std::vector& vars_bin_offsets, + i_t heavy_degree_cutoff, + i_t num_blocks_heavy_vars, + bool dry_run = false) +{ + if (num_blocks_heavy_vars != 0) { + auto heavy_vars_stream = streams.get_stream(); + // TODO : Check heavy_vars_block_segments size for profiling + if (!dry_run) { + auto heavy_vars_beg_id = get_id_offset(vars_bin_offsets, heavy_degree_cutoff); + lb_upd_bnd_heavy_kernel + <<>>( + heavy_vars_beg_id, + make_span(heavy_vars_vertex_ids), + make_span(heavy_vars_pseudo_block_ids), + heavy_degree_cutoff, + view, + upd_0, + upd_1); + auto num_heavy_vars = vars_bin_offsets.back() - heavy_vars_beg_id; + finalize_upd_bnd_kernel<<>>( + heavy_vars_beg_id, make_span(heavy_vars_block_segments), view, upd_0, upd_1); + lb_upd_next_constraint_heavy_kernel + <<>>( + heavy_vars_beg_id, + make_span(heavy_vars_vertex_ids), + make_span(heavy_vars_pseudo_block_ids), + heavy_degree_cutoff, + view, + upd_0, + upd_1); + } + } +} + +template +void upd_bounds_per_block(stream_pool_t& streams, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + const std::vector& vars_bin_offsets, + i_t degree_beg, + i_t degree_end, + bool dry_run) +{ + static_assert(block_dim <= 1024, "Cannot launch kernel with more than 1024 threads"); + + auto [vars_id_beg, vars_id_end] = get_id_range(vars_bin_offsets, degree_beg, degree_end); + + auto block_count = vars_id_end - vars_id_beg; + if (block_count > 0) { + auto block_stream = streams.get_stream(); + if (!dry_run) { + lb_upd_bnd_block_kernel + <<>>(vars_id_beg, view, upd_0, upd_1); + lb_upd_next_constraint_block_kernel + <<>>(vars_id_beg, view, upd_0, upd_1); + } + } +} + +template +void upd_bounds_per_block(stream_pool_t& streams, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + const std::vector& vars_bin_offsets, + i_t heavy_degree_cutoff, + bool dry_run = false) +{ + if (view.nnz < 10000) { + upd_bounds_per_block( + streams, view, upd_0, upd_1, vars_bin_offsets, 32, 32, dry_run); + upd_bounds_per_block( + streams, view, upd_0, upd_1, vars_bin_offsets, 64, 64, dry_run); + upd_bounds_per_block( + streams, view, upd_0, upd_1, vars_bin_offsets, 128, 128, dry_run); + upd_bounds_per_block( + streams, view, upd_0, upd_1, vars_bin_offsets, 256, 256, dry_run); + } else { + //[1024, heavy_degree_cutoff/2] -> 128 block size + upd_bounds_per_block( + streams, view, upd_0, upd_1, vars_bin_offsets, 1024, heavy_degree_cutoff / 2, dry_run); + //[64, 512] -> 32 block size + upd_bounds_per_block( + streams, view, upd_0, upd_1, vars_bin_offsets, 128, 512, dry_run); + } +} + +template +void upd_bounds_sub_warp(stream_pool_t& streams, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + i_t degree_beg, + i_t degree_end, + const std::vector& vars_bin_offsets, + bool dry_run) +{ + constexpr i_t block_dim = 32; + auto vars_per_block = block_dim / threads_per_variable; + auto [vars_id_beg, vars_id_end] = get_id_range(vars_bin_offsets, degree_beg, degree_end); + + auto block_count = raft::ceildiv(vars_id_end - vars_id_beg, vars_per_block); + if (block_count != 0) { + auto sub_warp_stream = streams.get_stream(); + if (!dry_run) { + lb_upd_bnd_sub_warp_kernel + <<>>( + vars_id_beg, vars_id_end, view, upd_0, upd_1); + lb_upd_next_constraint_sub_warp_kernel + <<>>( + vars_id_beg, vars_id_end, view, upd_0, upd_1); + } + } +} + +template +void upd_bounds_sub_warp(stream_pool_t& streams, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + i_t vars_sub_warp_count, + rmm::device_uvector& warp_vars_offsets, + rmm::device_uvector& warp_vars_id_offsets, + bool dry_run) +{ + constexpr i_t block_dim = 256; + + auto block_count = raft::ceildiv(vars_sub_warp_count * 32, block_dim); + if (block_count != 0) { + auto sub_warp_stream = streams.get_stream(); + if (!dry_run) { + lb_upd_bnd_sub_warp_kernel + <<>>( + view, upd_0, upd_1, make_span(warp_vars_offsets), make_span(warp_vars_id_offsets)); + lb_upd_next_constraint_sub_warp_kernel + <<>>( + view, upd_0, upd_1, make_span(warp_vars_offsets), make_span(warp_vars_id_offsets)); + } + } +} + +template +void upd_bounds_sub_warp(stream_pool_t& streams, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + i_t degree, + const std::vector& vars_bin_offsets, + bool dry_run) +{ + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, degree, degree, vars_bin_offsets, dry_run); +} + +template +void upd_bounds_sub_warp(stream_pool_t& streams, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + bool is_vars_sub_warp_single_bin, + i_t vars_sub_warp_count, + rmm::device_uvector& warp_vars_offsets, + rmm::device_uvector& warp_vars_id_offsets, + const std::vector& vars_bin_offsets, + bool dry_run = false) +{ + if (view.nnz < 10000) { + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, 16, vars_bin_offsets, dry_run); + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, 8, vars_bin_offsets, dry_run); + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, 4, vars_bin_offsets, dry_run); + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, 2, vars_bin_offsets, dry_run); + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, 1, vars_bin_offsets, dry_run); + } else { + if (is_vars_sub_warp_single_bin) { + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, 64, vars_bin_offsets, dry_run); + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, 32, vars_bin_offsets, dry_run); + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, 16, vars_bin_offsets, dry_run); + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, 8, vars_bin_offsets, dry_run); + upd_bounds_sub_warp( + streams, view, upd_0, upd_1, 1, 4, vars_bin_offsets, dry_run); + } else { + upd_bounds_sub_warp(streams, + view, + upd_0, + upd_1, + vars_sub_warp_count, + warp_vars_offsets, + warp_vars_id_offsets, + dry_run); + } + } +} + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/load_balanced_bounds_common_kernels.cuh b/cpp/src/mip/presolve/load_balanced_bounds_common_kernels.cuh new file mode 100644 index 0000000000..e887b1d160 --- /dev/null +++ b/cpp/src/mip/presolve/load_balanced_bounds_common_kernels.cuh @@ -0,0 +1,44 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2022-2025, NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: LicenseRef-NvidiaProprietary + * + * NVIDIA CORPORATION, its affiliates and licensors retain all intellectual + * property and proprietary rights in and to this material, related + * documentation and any modifications thereto. Any use, reproduction, + * disclosure or distribution of this material and related documentation + * without an express license agreement from NVIDIA CORPORATION or + * its affiliates is strictly prohibited. + */ + +#pragma once + +#include + +namespace cuopt::linear_programming::detail { + +template +__device__ __forceinline__ void detect_range_sub_warp(i_t* id_warp_beg, + i_t* id_range_end, + i_t* threads_per_item, + raft::device_span warp_offsets, + raft::device_span bin_offsets) +{ + i_t warp_id = (blockDim.x * blockIdx.x + threadIdx.x) / raft::WarpSize; + i_t lane_id = threadIdx.x & 31; + bool pred = false; + if (lane_id < warp_offsets.size()) { pred = (warp_id >= warp_offsets[lane_id]); } + unsigned int m = __ballot_sync(0xffffffff, pred); + i_t seg = 31 - __clz(m); + i_t it_per_warp = (1 << (5 - seg)); // item per warp = raft::WarpSize/(2^seg) + if (5 - seg < 0) { + *threads_per_item = 0; + return; + } + i_t beg = bin_offsets[seg] + (warp_id - warp_offsets[seg]) * it_per_warp; + i_t end = bin_offsets[seg + 1]; + *id_warp_beg = beg; + *id_range_end = end; + *threads_per_item = (1 << seg); +} + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/load_balanced_bounds_next_constraint_kernels.cuh b/cpp/src/mip/presolve/load_balanced_bounds_next_constraint_kernels.cuh new file mode 100644 index 0000000000..959c7d0148 --- /dev/null +++ b/cpp/src/mip/presolve/load_balanced_bounds_next_constraint_kernels.cuh @@ -0,0 +1,512 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2022-2025, NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: LicenseRef-NvidiaProprietary + * + * NVIDIA CORPORATION, its affiliates and licensors retain all intellectual + * property and proprietary rights in and to this material, related + * documentation and any modifications thereto. Any use, reproduction, + * disclosure or distribution of this material and related documentation + * without an express license agreement from NVIDIA CORPORATION or + * its affiliates is strictly prohibited. + */ + +#pragma once + +#include +#include "load_balanced_bounds_common_kernels.cuh" + +namespace cuopt::linear_programming::detail { + +/// BOUNDS UPDATE + +template +inline __device__ thrust::pair skip_mark(upd_view_t upd_0, + upd_view_t upd_1, + i_t var_idx) +{ + return thrust::make_pair((upd_0.var_bounds_changed[var_idx] == i_t{0}), + (upd_1.var_bounds_changed[var_idx] == i_t{0})); +} + +template +__device__ void update_next_constraints( + bounds_update_view_t view, upd_view_t upd, i_t tid, i_t beg, i_t end) +{ + for (i_t i = tid + beg; i < end; i += MAX_EDGE_PER_VAR) { + // auto a = view.coeff[i]; + auto cnst_idx = view.cnst[i]; + atomicExch(&upd.next_changed_constraints[cnst_idx], 1); + } +} + +template +__device__ void update_next_constraints( + bounds_update_view_t view, upd_view_t upd_0, upd_view_t upd_1, i_t tid, i_t beg, i_t end) +{ + for (i_t i = tid + beg; i < end; i += MAX_EDGE_PER_VAR) { + // auto a = view.coeff[i]; + auto cnst_idx = view.cnst[i]; + atomicExch(&upd_0.next_changed_constraints[cnst_idx], 1); + atomicExch(&upd_1.next_changed_constraints[cnst_idx], 1); + } +} + +template +__device__ void update_next_constraints(bounds_update_view_t view, i_t tid, i_t beg, i_t end) +{ + for (i_t i = tid + beg; i < end; i += MAX_EDGE_PER_VAR) { + // auto a = view.coeff[i]; + auto cnst_idx = view.cnst[i]; + atomicExch(&view.next_changed_constraints[cnst_idx], 1); + } +} + +template +__global__ void lb_upd_next_constraint_heavy_kernel(i_t id_range_beg, + raft::device_span ids, + raft::device_span pseudo_block_ids, + i_t work_per_block, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1) +{ + auto idx = ids[blockIdx.x] + id_range_beg; + auto pseudo_block_id = pseudo_block_ids[blockIdx.x]; + auto var_idx = view.vars_reorg_ids[idx]; + auto skip_mark_next = skip_mark(upd_0, upd_1, var_idx); + // x is lb, y is ub + // auto old_bounds = view.vars_bnd[var_idx]; + // bool is_int = (view.vars_types[idx] == var_t::INTEGER); + i_t item_off_beg = view.offsets[idx] + work_per_block * pseudo_block_id; + i_t item_off_end = min(item_off_beg + work_per_block, view.offsets[idx + 1]); + + // typedef cub::BlockReduce BlockReduce; + //__shared__ typename BlockReduce::TempStorage temp_storage; + + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + // if (old_bounds.x + view.tolerances.integrality_tolerance >= old_bounds.y) { + // tmp_bnd[blockIdx.x] = old_bounds; + // return; + //} + // auto bounds = + if (thrust::get<1>(skip_mark_next)) { + update_next_constraints(view, upd_0, threadIdx.x, item_off_beg, item_off_end); + } else if (thrust::get<0>(skip_mark_next)) { + update_next_constraints(view, upd_1, threadIdx.x, item_off_beg, item_off_end); + } else { + update_next_constraints( + view, upd_0, upd_1, threadIdx.x, item_off_beg, item_off_end); + } + + // bounds.x = BlockReduce(temp_storage).Reduce(bounds.x, cub::Max()); + //__syncthreads(); + // bounds.y = BlockReduce(temp_storage).Reduce(bounds.y, cub::Min()); + + // if (threadIdx.x == 0) { + // bool changed = write_updated_bounds(&tmp_bnd[blockIdx.x], is_int, view, bounds, old_bounds); + // atomicExch(&view.var_bounds_changed[var_idx], 1); + // } +} + +template +__global__ void lb_upd_next_constraint_heavy_kernel(i_t id_range_beg, + raft::device_span ids, + raft::device_span pseudo_block_ids, + i_t work_per_block, + bounds_update_view_t view) +{ + auto idx = ids[blockIdx.x] + id_range_beg; + auto pseudo_block_id = pseudo_block_ids[blockIdx.x]; + auto var_idx = view.vars_reorg_ids[idx]; + if (view.var_bounds_changed[var_idx] == 0) { return; } + // x is lb, y is ub + // auto old_bounds = view.vars_bnd[var_idx]; + // bool is_int = (view.vars_types[idx] == var_t::INTEGER); + i_t item_off_beg = view.offsets[idx] + work_per_block * pseudo_block_id; + i_t item_off_end = min(item_off_beg + work_per_block, view.offsets[idx + 1]); + + // typedef cub::BlockReduce BlockReduce; + //__shared__ typename BlockReduce::TempStorage temp_storage; + + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + // if (old_bounds.x + view.tolerances.integrality_tolerance >= old_bounds.y) { + // tmp_bnd[blockIdx.x] = old_bounds; + // return; + //} + // auto bounds = + update_next_constraints(view, threadIdx.x, item_off_beg, item_off_end); + + // bounds.x = BlockReduce(temp_storage).Reduce(bounds.x, cub::Max()); + //__syncthreads(); + // bounds.y = BlockReduce(temp_storage).Reduce(bounds.y, cub::Min()); + + // if (threadIdx.x == 0) { + // bool changed = write_updated_bounds(&tmp_bnd[blockIdx.x], is_int, view, bounds, old_bounds); + // atomicExch(&view.var_bounds_changed[var_idx], 1); + // } +} + +// template +//__global__ void finalize_upd_next_constraint_kernel(i_t heavy_vars_beg_id, +// raft::device_span item_offsets, +// raft::device_span tmp_bnd, +// bounds_update_view_t view) +//{ +// using warp_reduce = cub::WarpReduce; +// __shared__ typename warp_reduce::TempStorage temp_storage; +// i_t idx = heavy_vars_beg_id + blockIdx.x; +// i_t var_idx = view.vars_reorg_ids[idx]; +// if (view.var_bounds_changed[var_idx] == 0) { return; } +// +// // assumes cnst_bnd[i].x has ub and cnst_bnd[i].y has lb +// i_t item_off_beg = item_offsets[blockIdx.x]; +// i_t item_off_end = item_offsets[blockIdx.x + 1]; +// f_t2 bounds = f_t2{-std::numeric_limits::infinity(), +// std::numeric_limits::infinity()}; +// // assumes tmp_act[i].x has min activity and tmp_act[i].y has max activity +// for (i_t i = threadIdx.x + item_off_beg; i < item_off_end; i += blockDim.x) { +// auto bnd = tmp_bnd[i]; +// bounds.x = max(bounds.x, bnd.x); +// bounds.y = min(bounds.y, bnd.y); +// } +// bounds.x = warp_reduce(temp_storage).Reduce(bounds.x, cub::Max()); +// __syncwarp(); +// bounds.y = warp_reduce(temp_storage).Reduce(bounds.y, cub::Min()); +// if (threadIdx.x == 0) { view.vars_bnd[var_idx] = bounds; } +// } + +template +__global__ void lb_upd_next_constraint_block_kernel(i_t id_range_beg, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1) +{ + i_t idx = id_range_beg + blockIdx.x; + i_t var_idx = view.vars_reorg_ids[idx]; + auto skip_mark_next = skip_mark(upd_0, upd_1, var_idx); + // x is lb, y is ub + // auto old_bounds = view.vars_bnd[var_idx]; + // bool is_int = (view.vars_types[idx] == var_t::INTEGER); + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + + // typedef cub::BlockReduce BlockReduce; + //__shared__ typename BlockReduce::TempStorage temp_storage; + + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + // if (old_bounds.x + view.tolerances.integrality_tolerance >= old_bounds.y) { return; } + // auto bounds = + if (thrust::get<1>(skip_mark_next)) { + update_next_constraints(view, upd_0, threadIdx.x, item_off_beg, item_off_end); + } else if (thrust::get<0>(skip_mark_next)) { + update_next_constraints(view, upd_1, threadIdx.x, item_off_beg, item_off_end); + } else { + update_next_constraints( + view, upd_0, upd_1, threadIdx.x, item_off_beg, item_off_end); + } + + // bounds.x = BlockReduce(temp_storage).Reduce(bounds.x, cub::Max()); + //__syncthreads(); + // bounds.y = BlockReduce(temp_storage).Reduce(bounds.y, cub::Min()); + + // if (threadIdx.x == 0) { + // bool changed = write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, + // old_bounds); view.var_bounds_changed[var_idx] = changed; + // } +} + +template +__global__ void lb_upd_next_constraint_block_kernel(i_t id_range_beg, bounds_update_view_t view) +{ + i_t idx = id_range_beg + blockIdx.x; + i_t var_idx = view.vars_reorg_ids[idx]; + if (view.var_bounds_changed[var_idx] == 0) { return; } + // x is lb, y is ub + // auto old_bounds = view.vars_bnd[var_idx]; + // bool is_int = (view.vars_types[idx] == var_t::INTEGER); + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + + // typedef cub::BlockReduce BlockReduce; + //__shared__ typename BlockReduce::TempStorage temp_storage; + + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + // if (old_bounds.x + view.tolerances.integrality_tolerance >= old_bounds.y) { return; } + // auto bounds = + update_next_constraints(view, threadIdx.x, item_off_beg, item_off_end); + + // bounds.x = BlockReduce(temp_storage).Reduce(bounds.x, cub::Max()); + //__syncthreads(); + // bounds.y = BlockReduce(temp_storage).Reduce(bounds.y, cub::Min()); + + // if (threadIdx.x == 0) { + // bool changed = write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, + // old_bounds); view.var_bounds_changed[var_idx] = changed; + // } +} + +template +__global__ void lb_upd_next_constraint_sub_warp_kernel( + i_t id_range_beg, i_t id_range_end, activity_view_t view, upd_view_t upd_0, upd_view_t upd_1) +{ + constexpr i_t ids_per_block = BDIM / MAX_EDGE_PER_VAR; + i_t id_beg = blockIdx.x * ids_per_block + id_range_beg; + i_t idx = id_beg + (threadIdx.x / MAX_EDGE_PER_VAR); + i_t var_idx; + // auto old_bounds = + // f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + // auto bounds = old_bounds; + // bool is_int = false; + auto skip_mark_next = thrust::make_pair(!(idx < id_range_end), !(idx < id_range_end)); + if (idx < id_range_end) { + var_idx = view.vars_reorg_ids[idx]; + skip_mark_next = skip_mark(upd_0, upd_1, var_idx); + // old_bounds = view.vars_bnd[var_idx]; + // is_int = (view.vars_types[idx] == var_t::INTEGER); + } + i_t p_tid = threadIdx.x % MAX_EDGE_PER_VAR; + + // i_t head_flag = (p_tid == 0); + + // using warp_reduce = cub::WarpReduce; + //__shared__ typename warp_reduce::TempStorage temp_storage; + + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + if (thrust::get<1>(skip_mark_next)) { + update_next_constraints( + view, upd_0, p_tid, item_off_beg, item_off_end); + } else if (thrust::get<0>(skip_mark_next)) { + update_next_constraints( + view, upd_1, p_tid, item_off_beg, item_off_end); + } else { + update_next_constraints( + view, upd_0, upd_1, p_tid, item_off_beg, item_off_end); + } + + // bounds.x = warp_reduce(temp_storage).Reduce(bounds.x, cub::Max()); + //__syncwarp(); + // bounds.y = warp_reduce(temp_storage).Reduce(bounds.y, cub::Min()); + + // if (head_flag && continue_calc) { + // bool changed = write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, + // old_bounds); view.var_bounds_changed[var_idx] = changed; + // } +} + +template +__global__ void lb_upd_next_constraint_sub_warp_kernel(i_t id_range_beg, + i_t id_range_end, + activity_view_t view) +{ + constexpr i_t ids_per_block = BDIM / MAX_EDGE_PER_VAR; + i_t id_beg = blockIdx.x * ids_per_block + id_range_beg; + i_t idx = id_beg + (threadIdx.x / MAX_EDGE_PER_VAR); + i_t var_idx; + // auto old_bounds = + // f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + // auto bounds = old_bounds; + // bool is_int = false; + bool continue_calc = (idx < id_range_end); + if (continue_calc) { + var_idx = view.vars_reorg_ids[idx]; + if (view.var_bounds_changed[var_idx] == 0) { continue_calc = false; } + // old_bounds = view.vars_bnd[var_idx]; + // is_int = (view.vars_types[idx] == var_t::INTEGER); + } + i_t p_tid = threadIdx.x % MAX_EDGE_PER_VAR; + + // i_t head_flag = (p_tid == 0); + + // using warp_reduce = cub::WarpReduce; + //__shared__ typename warp_reduce::TempStorage temp_storage; + + if (continue_calc) { + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + // if (old_bounds.x + view.tolerances.integrality_tolerance >= old_bounds.y) { head_flag = 0; } + + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + + update_next_constraints(view, p_tid, item_off_beg, item_off_end); + } + + // bounds.x = warp_reduce(temp_storage).Reduce(bounds.x, cub::Max()); + //__syncwarp(); + // bounds.y = warp_reduce(temp_storage).Reduce(bounds.y, cub::Min()); + + // if (head_flag && continue_calc) { + // bool changed = write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, + // old_bounds); view.var_bounds_changed[var_idx] = changed; + // } +} + +template +__device__ void upd_next_constraint_sub_warp(i_t id_warp_beg, + i_t id_range_end, + bounds_update_view_t view) +{ + i_t lane_id = (threadIdx.x & 31); + i_t idx = id_warp_beg + (lane_id / MAX_EDGE_PER_VAR); + i_t var_idx; + // auto old_bounds = + // f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + // auto bounds = old_bounds; + // bool is_int = false; + bool continue_calc = (idx < id_range_end); + if (continue_calc) { + var_idx = view.vars_reorg_ids[idx]; + if (view.var_bounds_changed[var_idx] == 0) { continue_calc = false; } + // old_bounds = view.vars_bnd[var_idx]; + // is_int = (view.vars_types[idx] == var_t::INTEGER); + } + // Equivalent to + // i_t p_tid = threadIdx.x % MAX_EDGE_PER_VAR; + i_t p_tid = lane_id & (MAX_EDGE_PER_VAR - 1); + + // i_t head_flag = (p_tid == 0); + + // using warp_reduce = cub::WarpReduce; + //__shared__ typename warp_reduce::TempStorage temp_storage; + + if (continue_calc) { + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + // if (old_bounds.x + view.tolerances.integrality_tolerance >= old_bounds.y) { head_flag = 0; } + + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + + update_next_constraints(view, p_tid, item_off_beg, item_off_end); + } + + // bounds.x = warp_reduce(temp_storage).Reduce(bounds.x, cub::Max()); + //__syncwarp(); + // bounds.y = warp_reduce(temp_storage).Reduce(bounds.y, cub::Min()); + + // if (head_flag && continue_calc) { + // bool changed = write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, + // old_bounds); view.var_bounds_changed[var_idx] = changed; + // } +} + +template +__device__ void upd_next_constraint_sub_warp( + i_t id_warp_beg, i_t id_range_end, bounds_update_view_t view, upd_view_t upd_0, upd_view_t upd_1) +{ + i_t lane_id = (threadIdx.x & 31); + i_t idx = id_warp_beg + (lane_id / MAX_EDGE_PER_VAR); + i_t var_idx; + // auto old_bounds = + // f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + // auto bounds = old_bounds; + // bool is_int = false; + auto skip_mark_next = thrust::make_pair(!(idx < id_range_end), !(idx < id_range_end)); + if (idx < id_range_end) { + var_idx = view.vars_reorg_ids[idx]; + skip_mark_next = skip_mark(upd_0, upd_1, var_idx); + // old_bounds = view.vars_bnd[var_idx]; + // is_int = (view.vars_types[idx] == var_t::INTEGER); + } + // Equivalent to + // i_t p_tid = threadIdx.x % MAX_EDGE_PER_VAR; + i_t p_tid = lane_id & (MAX_EDGE_PER_VAR - 1); + + // i_t head_flag = (p_tid == 0); + + // using warp_reduce = cub::WarpReduce; + //__shared__ typename warp_reduce::TempStorage temp_storage; + + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + if (thrust::get<1>(skip_mark_next)) { + update_next_constraints( + view, upd_0, p_tid, item_off_beg, item_off_end); + } else if (thrust::get<0>(skip_mark_next)) { + update_next_constraints( + view, upd_1, p_tid, item_off_beg, item_off_end); + } else { + update_next_constraints( + view, upd_0, upd_1, p_tid, item_off_beg, item_off_end); + } + + // bounds.x = warp_reduce(temp_storage).Reduce(bounds.x, cub::Max()); + //__syncwarp(); + // bounds.y = warp_reduce(temp_storage).Reduce(bounds.y, cub::Min()); + + // if (head_flag && continue_calc) { + // bool changed = write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, + // old_bounds); view.var_bounds_changed[var_idx] = changed; + // } +} + +template +__global__ void lb_upd_next_constraint_sub_warp_kernel(bounds_update_view_t view, + raft::device_span warp_vars_offsets, + raft::device_span warp_vars_id_offsets) +{ + i_t id_warp_beg, id_range_end, threads_per_variable; + detect_range_sub_warp( + &id_warp_beg, &id_range_end, &threads_per_variable, warp_vars_offsets, warp_vars_id_offsets); + + if (threads_per_variable == 1) { + upd_next_constraint_sub_warp(id_warp_beg, id_range_end, view); + } else if (threads_per_variable == 2) { + upd_next_constraint_sub_warp(id_warp_beg, id_range_end, view); + } else if (threads_per_variable == 4) { + upd_next_constraint_sub_warp(id_warp_beg, id_range_end, view); + } else if (threads_per_variable == 8) { + upd_next_constraint_sub_warp(id_warp_beg, id_range_end, view); + } else if (threads_per_variable == 16) { + upd_next_constraint_sub_warp(id_warp_beg, id_range_end, view); + } +} + +template +__global__ void lb_upd_next_constraint_sub_warp_kernel(bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + raft::device_span warp_vars_offsets, + raft::device_span warp_vars_id_offsets) +{ + i_t id_warp_beg, id_range_end, threads_per_variable; + detect_range_sub_warp( + &id_warp_beg, &id_range_end, &threads_per_variable, warp_vars_offsets, warp_vars_id_offsets); + + if (threads_per_variable == 1) { + upd_next_constraint_sub_warp(id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_variable == 2) { + upd_next_constraint_sub_warp(id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_variable == 4) { + upd_next_constraint_sub_warp(id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_variable == 8) { + upd_next_constraint_sub_warp(id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_variable == 16) { + upd_next_constraint_sub_warp(id_warp_beg, id_range_end, view, upd_0, upd_1); + } +} + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/load_balanced_bounds_presolve.cu b/cpp/src/mip/presolve/load_balanced_bounds_presolve.cu index 7f98a42cb0..0fff8d66e5 100644 --- a/cpp/src/mip/presolve/load_balanced_bounds_presolve.cu +++ b/cpp/src/mip/presolve/load_balanced_bounds_presolve.cu @@ -56,6 +56,10 @@ load_balanced_bounds_presolve_t::load_balanced_bounds_presolve_t( : streams(max_stream_count_), pb(&problem_), bounds_changed(problem_.handle_ptr->get_stream()), + var_bounds_changed(0, problem_.handle_ptr->get_stream()), + changed_constraints(0, problem_.handle_ptr->get_stream()), + changed_variables(0, problem_.handle_ptr->get_stream()), + next_changed_constraints(0, problem_.handle_ptr->get_stream()), cnst_slack(0, problem_.handle_ptr->get_stream()), vars_bnd(0, problem_.handle_ptr->get_stream()), tmp_act(0, problem_.handle_ptr->get_stream()), @@ -94,44 +98,6 @@ load_balanced_bounds_presolve_t::~load_balanced_bounds_presolve_t() if (upd_bnd_graph_created) { cudaGraphExecDestroy(upd_bnd_exec); } } -template -std::pair sub_warp_meta(rmm::cuda_stream_view stream, - rmm::device_uvector& d_warp_offsets, - rmm::device_uvector& d_warp_id_offsets, - const std::vector& bin_offsets, - i_t w_t_r) -{ - // 1, 2, 4, 8, 16 - auto sub_warp_bin_count = 5; - std::vector warp_counts(sub_warp_bin_count); - - std::vector warp_offsets(warp_counts.size() + 1); - std::vector warp_id_offsets(warp_counts.size() + 1); - - for (size_t i = 0; i < warp_id_offsets.size(); ++i) { - warp_id_offsets[i] = bin_offsets[i + std::log2(w_t_r) + 1]; - } - warp_id_offsets[0] = bin_offsets[0]; - - i_t non_empty_bin_count = 0; - for (size_t i = 0; i < warp_counts.size(); ++i) { - warp_counts[i] = - raft::ceildiv((warp_id_offsets[i + 1] - warp_id_offsets[i]) * (1 << i), raft::WarpSize); - if (warp_counts[i] != 0) { non_empty_bin_count++; } - } - - warp_offsets[0] = 0; - for (size_t i = 1; i < warp_offsets.size(); ++i) { - warp_offsets[i] = warp_offsets[i - 1] + warp_counts[i - 1]; - } - expand_device_copy(d_warp_offsets, warp_offsets, stream); - expand_device_copy(d_warp_id_offsets, warp_id_offsets, stream); - - // If there is only 1 bin active, then there is no need to add logic to determine which warps work - // on which bin - return std::make_pair(non_empty_bin_count == 1, warp_offsets.back()); -} - template void load_balanced_bounds_presolve_t::copy_input_bounds( const load_balanced_problem_t& problem) @@ -165,55 +131,6 @@ void load_balanced_bounds_presolve_t::update_device_bounds( raft::copy(vars_bnd.data(), host_bounds.data(), host_bounds.size(), handle_ptr->get_stream()); } -template -bool build_graph(managed_stream_pool& streams, - const raft::handle_t* handle_ptr, - cudaGraph_t& graph, - cudaGraphExec_t& graph_exec, - DryRunFunc d_func, - CaptureGraphFunc g_func) -{ - bool graph_created = false; - cudaGraphCreate(&graph, 0); - cudaEvent_t fork_stream_event; - cudaEventCreate(&fork_stream_event); - - cudaStreamBeginCapture(handle_ptr->get_stream(), cudaStreamCaptureModeThreadLocal); - cudaEventRecord(fork_stream_event, handle_ptr->get_stream()); - - // dry-run - managed pool tracks how many streams were issued - d_func(); - streams.wait_issued_on_event(fork_stream_event); - streams.reset_issued(); - - g_func(); - auto activity_done = streams.create_events_on_issued(); - streams.reset_issued(); - for (auto& e : activity_done) { - cudaStreamWaitEvent(handle_ptr->get_stream(), e); - } - - cudaStreamEndCapture(handle_ptr->get_stream(), &graph); - RAFT_CHECK_CUDA(handle_ptr->get_stream()); - - if (graph_exec != nullptr) { - cudaGraphExecDestroy(graph_exec); - cudaGraphInstantiate(&graph_exec, graph, NULL, NULL, 0); - RAFT_CHECK_CUDA(handle_ptr->get_stream()); - } else { - cudaGraphInstantiate(&graph_exec, graph, NULL, NULL, 0); - RAFT_CHECK_CUDA(handle_ptr->get_stream()); - } - - cudaGraphDestroy(graph); - graph_created = true; - - handle_ptr->get_stream().synchronize(); - RAFT_CHECK_CUDA(handle_ptr->get_stream()); - - return graph_created; -} - template void load_balanced_bounds_presolve_t::setup( const load_balanced_problem_t& problem) @@ -223,6 +140,10 @@ void load_balanced_bounds_presolve_t::setup( auto stream = handle_ptr->get_stream(); stream.synchronize(); host_bounds.resize(2 * pb->n_variables); + var_bounds_changed.resize(pb->n_variables, stream); + changed_constraints.resize(pb->n_constraints, stream); + changed_variables.resize(pb->n_variables, stream); + next_changed_constraints.resize(pb->n_constraints, stream); cnst_slack.resize(2 * pb->n_constraints, stream); vars_bnd.resize(2 * pb->n_variables, stream); calc_slack_graph_created = false; @@ -300,15 +221,19 @@ load_balanced_bounds_presolve_t::get_activity_view( const load_balanced_problem_t& pb) { load_balanced_bounds_presolve_t::activity_view_t v; - v.cnst_reorg_ids = make_span(pb.cnst_reorg_ids); - v.coeff = make_span(pb.coefficients); - v.vars = make_span(pb.variables); - v.offsets = make_span(pb.offsets); - v.cnst_bnd = make_span_2(pb.cnst_bounds_data); - v.vars_bnd = make_span_2(vars_bnd); - v.cnst_slack = make_span_2(cnst_slack); - v.nnz = pb.nnz; - v.tolerances = pb.tolerances; + v.cnst_reorg_ids = make_span(pb.cnst_reorg_ids); + v.coeff = make_span(pb.coefficients); + v.vars = make_span(pb.variables); + v.offsets = make_span(pb.offsets); + v.cnst_bnd = make_span_2(pb.cnst_bounds_data); + v.vars_bnd = make_span_2(vars_bnd); + v.cnst_slack = make_span_2(cnst_slack); + v.var_bounds_changed = make_span(var_bounds_changed); + v.changed_constraints = make_span(changed_constraints); + v.changed_variables = make_span(changed_variables); + v.next_changed_constraints = make_span(next_changed_constraints); + v.nnz = pb.nnz; + v.tolerances = pb.tolerances; return v; } @@ -318,16 +243,20 @@ load_balanced_bounds_presolve_t::get_bounds_update_view( const load_balanced_problem_t& pb) { load_balanced_bounds_presolve_t::bounds_update_view_t v; - v.vars_reorg_ids = make_span(pb.vars_reorg_ids); - v.coeff = make_span(pb.reverse_coefficients); - v.cnst = make_span(pb.reverse_constraints); - v.offsets = make_span(pb.reverse_offsets); - v.vars_types = make_span(pb.vars_types); - v.vars_bnd = make_span_2(vars_bnd); - v.cnst_slack = make_span_2(cnst_slack); - v.bounds_changed = bounds_changed.data(); - v.nnz = pb.nnz; - v.tolerances = pb.tolerances; + v.vars_reorg_ids = make_span(pb.vars_reorg_ids); + v.coeff = make_span(pb.reverse_coefficients); + v.cnst = make_span(pb.reverse_constraints); + v.offsets = make_span(pb.reverse_offsets); + v.vars_types = make_span(pb.vars_types); + v.vars_bnd = make_span_2(vars_bnd); + v.cnst_slack = make_span_2(cnst_slack); + v.bounds_changed = bounds_changed.data(); + v.var_bounds_changed = make_span(var_bounds_changed); + v.changed_constraints = make_span(changed_constraints); + v.changed_variables = make_span(changed_variables); + v.next_changed_constraints = make_span(next_changed_constraints); + v.nnz = pb.nnz; + v.tolerances = pb.tolerances; return v; } @@ -436,6 +365,38 @@ bool load_balanced_bounds_presolve_t::update_bounds_from_slack( return (zero < h_bounds_changed); } +template +void load_balanced_bounds_presolve_t::init_changed_constraints( + const raft::handle_t* handle_ptr) +{ + thrust::fill( + handle_ptr->get_thrust_policy(), var_bounds_changed.begin(), var_bounds_changed.end(), 0); + thrust::fill( + handle_ptr->get_thrust_policy(), changed_variables.begin(), changed_variables.end(), 1); + thrust::fill( + handle_ptr->get_thrust_policy(), changed_constraints.begin(), changed_constraints.end(), 1); + thrust::fill(handle_ptr->get_thrust_policy(), + next_changed_constraints.begin(), + next_changed_constraints.end(), + 0); +} + +template +void load_balanced_bounds_presolve_t::prepare_for_next_iteration( + const raft::handle_t* handle_ptr) +{ + std::swap(changed_constraints, next_changed_constraints); + handle_ptr->sync_stream(); + thrust::fill(handle_ptr->get_thrust_policy(), + next_changed_constraints.begin(), + next_changed_constraints.end(), + 0); + thrust::fill( + handle_ptr->get_thrust_policy(), changed_variables.begin(), changed_variables.end(), 0); + thrust::fill( + handle_ptr->get_thrust_policy(), var_bounds_changed.begin(), var_bounds_changed.end(), 0); +} + template termination_criterion_t load_balanced_bounds_presolve_t::bound_update_loop( const raft::handle_t* handle_ptr, timer_t timer) @@ -443,6 +404,7 @@ termination_criterion_t load_balanced_bounds_presolve_t::bound_update_ termination_criterion_t criteria = termination_criterion_t::ITERATION_LIMIT; i_t iter; + init_changed_constraints(handle_ptr); for (iter = 0; iter < settings.iteration_limit; ++iter) { calculate_constraint_slack_iter(handle_ptr); if (!update_bounds_from_slack(handle_ptr)) { @@ -453,6 +415,7 @@ termination_criterion_t load_balanced_bounds_presolve_t::bound_update_ } break; } + prepare_for_next_iteration(handle_ptr); if (timer.check_time_limit()) { criteria = termination_criterion_t::TIME_LIMIT; CUOPT_LOG_DEBUG("Exiting bounds prop because of time limit at iter %d", iter); @@ -600,6 +563,20 @@ void load_balanced_bounds_presolve_t::set_updated_bounds( problem->variable_bounds.data(), vars_bnd.data(), vars_bnd.size(), handle_ptr->get_stream()); } +template +void load_balanced_bounds_presolve_t::set_updated_bounds(rmm::device_uvector& lb, + rmm::device_uvector& ub) +{ + auto& handle_ptr = pb->handle_ptr; + auto out = thrust::make_zip_iterator(thrust::make_tuple(lb.begin(), ub.begin())); + auto bnd_span = make_span_2(vars_bnd); + thrust::transform(handle_ptr->get_thrust_policy(), + bnd_span.begin(), + bnd_span.end(), + out, + [] __device__(auto bnd) { return thrust::make_tuple(bnd.x, bnd.y); }); +} + #if MIP_INSTANTIATE_FLOAT template class load_balanced_bounds_presolve_t; #endif diff --git a/cpp/src/mip/presolve/load_balanced_bounds_presolve.cuh b/cpp/src/mip/presolve/load_balanced_bounds_presolve.cuh index e433821bff..38d654ad13 100644 --- a/cpp/src/mip/presolve/load_balanced_bounds_presolve.cuh +++ b/cpp/src/mip/presolve/load_balanced_bounds_presolve.cuh @@ -108,7 +108,7 @@ class managed_stream_pool { { for (int i = 0; i < end_unsycned + 1; ++i) { streams_[i].synchronize(); - RAFT_CHECK_CUDA(streams_[i].value()); + RAFT_CHECK_CUDA(streams_[i]); } end_unsycned = -1; next_stream = 0; @@ -159,6 +159,8 @@ class load_balanced_bounds_presolve_t { void calculate_activity_graph(bool erase_inf_cnst, bool dry_run = false); void calculate_bounds_update_graph(bool dry_run = false); + void init_changed_constraints(const raft::handle_t* handle_ptr); + void prepare_for_next_iteration(const raft::handle_t* handle_ptr); void calculate_constraint_slack(const raft::handle_t* handle_ptr); void calculate_constraint_slack_iter(const raft::handle_t* handle_ptr); @@ -179,6 +181,7 @@ class load_balanced_bounds_presolve_t { void set_bounds(const std::vector>& var_probe_vals, const raft::handle_t* handle_ptr); void set_updated_bounds(load_balanced_problem_t* problem); + void set_updated_bounds(rmm::device_uvector& lb, rmm::device_uvector& ub); struct activity_view_t { raft::device_span cnst_reorg_ids; @@ -188,6 +191,10 @@ class load_balanced_bounds_presolve_t { raft::device_span cnst_bnd; // new indexing raft::device_span vars_bnd; // old indexing raft::device_span cnst_slack; // old indexing + raft::device_span var_bounds_changed; + raft::device_span changed_constraints; + raft::device_span changed_variables; + raft::device_span next_changed_constraints; i_t nnz; typename mip_solver_settings_t::tolerances_t tolerances; }; @@ -200,6 +207,10 @@ class load_balanced_bounds_presolve_t { raft::device_span vars_types; // new indexing raft::device_span vars_bnd; // old indexing raft::device_span cnst_slack; // old indexing + raft::device_span var_bounds_changed; + raft::device_span changed_constraints; + raft::device_span changed_variables; + raft::device_span next_changed_constraints; i_t* bounds_changed; i_t nnz; typename mip_solver_settings_t::tolerances_t tolerances; @@ -217,6 +228,11 @@ class load_balanced_bounds_presolve_t { rmm::device_scalar bounds_changed; + rmm::device_uvector var_bounds_changed; + rmm::device_uvector changed_constraints; + rmm::device_uvector changed_variables; + rmm::device_uvector next_changed_constraints; + rmm::device_uvector cnst_slack; rmm::device_uvector vars_bnd; rmm::device_uvector tmp_act; diff --git a/cpp/src/mip/presolve/load_balanced_bounds_presolve_helpers.cuh b/cpp/src/mip/presolve/load_balanced_bounds_presolve_helpers.cuh index e8ce0b52e0..82b4397704 100644 --- a/cpp/src/mip/presolve/load_balanced_bounds_presolve_helpers.cuh +++ b/cpp/src/mip/presolve/load_balanced_bounds_presolve_helpers.cuh @@ -12,6 +12,7 @@ #pragma once +#include "load_balanced_bounds_next_constraint_kernels.cuh" #include "load_balanced_bounds_presolve_kernels.cuh" #include "load_balanced_partition_helpers.cuh" @@ -25,6 +26,55 @@ namespace cuopt::linear_programming::detail { +template +bool build_graph(managed_stream_pool& streams, + const raft::handle_t* handle_ptr, + cudaGraph_t& graph, + cudaGraphExec_t& graph_exec, + DryRunFunc d_func, + CaptureGraphFunc g_func) +{ + bool graph_created = false; + cudaGraphCreate(&graph, 0); + cudaEvent_t fork_stream_event; + cudaEventCreate(&fork_stream_event); + + // dry-run - managed pool tracks how many streams were issued + d_func(); + cudaStreamBeginCapture(handle_ptr->get_stream(), cudaStreamCaptureModeGlobal); + cudaEventRecord(fork_stream_event, handle_ptr->get_stream()); + + streams.wait_issued_on_event(fork_stream_event); + streams.reset_issued(); + + g_func(); + auto activity_done = streams.create_events_on_issued(); + streams.reset_issued(); + for (auto& e : activity_done) { + cudaStreamWaitEvent(handle_ptr->get_stream(), e); + } + + cudaStreamEndCapture(handle_ptr->get_stream(), &graph); + RAFT_CHECK_CUDA(handle_ptr->get_stream()); + + if (graph_exec != nullptr) { + cudaGraphExecDestroy(graph_exec); + cudaGraphInstantiate(&graph_exec, graph, NULL, NULL, 0); + RAFT_CHECK_CUDA(handle_ptr->get_stream()); + } else { + cudaGraphInstantiate(&graph_exec, graph, NULL, NULL, 0); + RAFT_CHECK_CUDA(handle_ptr->get_stream()); + } + + cudaGraphDestroy(graph); + graph_created = true; + + handle_ptr->get_stream().synchronize(); + RAFT_CHECK_CUDA(handle_ptr->get_stream()); + + return graph_created; +} + template i_t get_id_offset(const std::vector& bin_offsets, i_t degree_cutoff) { @@ -79,6 +129,44 @@ struct heavy_vertex_meta_t : public thrust::unary_function { } }; +template +std::pair sub_warp_meta(rmm::cuda_stream_view stream, + rmm::device_uvector& d_warp_offsets, + rmm::device_uvector& d_warp_id_offsets, + const std::vector& bin_offsets, + i_t w_t_r) +{ + // 1, 2, 4, 8, 16 + auto sub_warp_bin_count = 5; + std::vector warp_counts(sub_warp_bin_count); + + std::vector warp_offsets(warp_counts.size() + 1); + std::vector warp_id_offsets(warp_counts.size() + 1); + + for (size_t i = 0; i < warp_id_offsets.size(); ++i) { + warp_id_offsets[i] = bin_offsets[i + std::log2(w_t_r) + 1]; + } + warp_id_offsets[0] = bin_offsets[0]; + + i_t non_empty_bin_count = 0; + for (size_t i = 0; i < warp_counts.size(); ++i) { + warp_counts[i] = + raft::ceildiv((warp_id_offsets[i + 1] - warp_id_offsets[i]) * (1 << i), raft::WarpSize); + if (warp_counts[i] != 0) { non_empty_bin_count++; } + } + + warp_offsets[0] = 0; + for (size_t i = 1; i < warp_offsets.size(); ++i) { + warp_offsets[i] = warp_offsets[i - 1] + warp_counts[i - 1]; + } + expand_device_copy(d_warp_offsets, warp_offsets, stream); + expand_device_copy(d_warp_id_offsets, warp_id_offsets, stream); + + // If there is only 1 bin active, then there is no need to add logic to determine which warps work + // on which bin + return std::make_pair(non_empty_bin_count == 1, warp_offsets.back()); +} + template i_t create_heavy_item_block_segments(rmm::cuda_stream_view stream, rmm::device_uvector& vertex_id, @@ -135,8 +223,13 @@ i_t create_heavy_item_block_segments(rmm::cuda_stream_view stream, /// CALCULATE ACTIVITY -template -void calc_activity_heavy_cnst(managed_stream_pool& streams, +template +void calc_activity_heavy_cnst(stream_pool_t& streams, activity_view_t view, raft::device_span tmp_cnst_act, const rmm::device_uvector& heavy_cnst_vertex_ids, @@ -175,8 +268,13 @@ void calc_activity_heavy_cnst(managed_stream_pool& streams, } } -template -void calc_activity_per_block(managed_stream_pool& streams, +template +void calc_activity_per_block(stream_pool_t& streams, activity_view_t view, const std::vector& cnst_bin_offsets, i_t degree_beg, @@ -203,8 +301,12 @@ void calc_activity_per_block(managed_stream_pool& streams, } } -template -void calc_activity_per_block(managed_stream_pool& streams, +template +void calc_activity_per_block(stream_pool_t& streams, activity_view_t view, const std::vector& cnst_bin_offsets, i_t heavy_degree_cutoff, @@ -234,8 +336,9 @@ template -void calc_activity_sub_warp(managed_stream_pool& streams, + typename activity_view_t, + typename stream_pool_t> +void calc_activity_sub_warp(stream_pool_t& streams, activity_view_t view, i_t degree_beg, i_t degree_end, @@ -266,8 +369,9 @@ template -void calc_activity_sub_warp(managed_stream_pool& streams, + typename activity_view_t, + typename stream_pool_t> +void calc_activity_sub_warp(stream_pool_t& streams, activity_view_t view, i_t degree, const std::vector& cnst_bin_offsets, @@ -278,8 +382,12 @@ void calc_activity_sub_warp(managed_stream_pool& streams, streams, view, degree, degree, cnst_bin_offsets, erase_inf_cnst, dry_run); } -template -void calc_activity_sub_warp(managed_stream_pool& streams, +template +void calc_activity_sub_warp(stream_pool_t& streams, activity_view_t view, i_t cnst_sub_warp_count, rmm::device_uvector& warp_cnst_offsets, @@ -306,8 +414,12 @@ void calc_activity_sub_warp(managed_stream_pool& streams, } } -template -void calc_activity_sub_warp(managed_stream_pool& streams, +template +void calc_activity_sub_warp(stream_pool_t& streams, activity_view_t view, bool is_cnst_sub_warp_single_bin, i_t cnst_sub_warp_count, @@ -354,8 +466,13 @@ void calc_activity_sub_warp(managed_stream_pool& streams, /// BOUNDS UPDATE -template -void upd_bounds_heavy_vars(managed_stream_pool& streams, +template +void upd_bounds_heavy_vars(stream_pool_t& streams, bounds_update_view_t view, raft::device_span tmp_vars_bnd, const rmm::device_uvector& heavy_vars_vertex_ids, @@ -382,41 +499,24 @@ void upd_bounds_heavy_vars(managed_stream_pool& streams, auto num_heavy_vars = vars_bin_offsets.back() - heavy_vars_beg_id; finalize_upd_bnd_kernel<<>>( heavy_vars_beg_id, make_span(heavy_vars_block_segments), tmp_vars_bnd, view); - } - } -} - -template -void upd_bounds_heavy_vars(managed_stream_pool& streams, - bounds_update_view_t view, - raft::device_span tmp_vars_bnd, - const rmm::device_uvector& heavy_vars_block_segments, - const std::vector& vars_bin_offsets, - i_t heavy_degree_cutoff, - i_t num_blocks_heavy_vars, - bool dry_run = false) -{ - if (num_blocks_heavy_vars != 0) { - auto heavy_vars_stream = streams.get_stream(); - // TODO : Check heavy_vars_block_segments size for profiling - if (!dry_run) { - auto heavy_vars_beg_id = get_id_offset(vars_bin_offsets, heavy_degree_cutoff); - lb_upd_bnd_heavy_kernel + lb_upd_next_constraint_heavy_kernel <<>>( heavy_vars_beg_id, - make_span(heavy_vars_block_segments, 1, heavy_vars_block_segments.size()), + make_span(heavy_vars_vertex_ids), + make_span(heavy_vars_pseudo_block_ids), heavy_degree_cutoff, - view, - tmp_vars_bnd); - auto num_heavy_vars = vars_bin_offsets.back() - heavy_vars_beg_id; - finalize_upd_bnd_kernel<<>>( - heavy_vars_beg_id, make_span(heavy_vars_block_segments), tmp_vars_bnd, view); + view); } } } -template -void upd_bounds_per_block(managed_stream_pool& streams, +template +void upd_bounds_per_block(stream_pool_t& streams, bounds_update_view_t view, const std::vector& vars_bin_offsets, i_t degree_beg, @@ -433,12 +533,18 @@ void upd_bounds_per_block(managed_stream_pool& streams, if (!dry_run) { lb_upd_bnd_block_kernel <<>>(vars_id_beg, view); + lb_upd_next_constraint_block_kernel + <<>>(vars_id_beg, view); } } } -template -void upd_bounds_per_block(managed_stream_pool& streams, +template +void upd_bounds_per_block(stream_pool_t& streams, bounds_update_view_t view, const std::vector& vars_bin_offsets, i_t heavy_degree_cutoff, @@ -462,8 +568,9 @@ template -void upd_bounds_sub_warp(managed_stream_pool& streams, + typename bounds_update_view_t, + typename stream_pool_t> +void upd_bounds_sub_warp(stream_pool_t& streams, bounds_update_view_t view, i_t degree_beg, i_t degree_end, @@ -480,12 +587,18 @@ void upd_bounds_sub_warp(managed_stream_pool& streams, if (!dry_run) { lb_upd_bnd_sub_warp_kernel <<>>(vars_id_beg, vars_id_end, view); + lb_upd_next_constraint_sub_warp_kernel + <<>>(vars_id_beg, vars_id_end, view); } } } -template -void upd_bounds_sub_warp(managed_stream_pool& streams, +template +void upd_bounds_sub_warp(stream_pool_t& streams, bounds_update_view_t view, i_t vars_sub_warp_count, rmm::device_uvector& warp_vars_offsets, @@ -501,6 +614,9 @@ void upd_bounds_sub_warp(managed_stream_pool& streams, lb_upd_bnd_sub_warp_kernel <<>>( view, make_span(warp_vars_offsets), make_span(warp_vars_id_offsets)); + lb_upd_next_constraint_sub_warp_kernel + <<>>( + view, make_span(warp_vars_offsets), make_span(warp_vars_id_offsets)); } } } @@ -509,8 +625,9 @@ template -void upd_bounds_sub_warp(managed_stream_pool& streams, + typename bounds_update_view_t, + typename stream_pool_t> +void upd_bounds_sub_warp(stream_pool_t& streams, bounds_update_view_t view, i_t degree, const std::vector& vars_bin_offsets, @@ -520,8 +637,12 @@ void upd_bounds_sub_warp(managed_stream_pool& streams, streams, view, degree, degree, vars_bin_offsets, dry_run); } -template -void upd_bounds_sub_warp(managed_stream_pool& streams, +template +void upd_bounds_sub_warp(stream_pool_t& streams, bounds_update_view_t view, bool is_vars_sub_warp_single_bin, i_t vars_sub_warp_count, @@ -549,4 +670,178 @@ void upd_bounds_sub_warp(managed_stream_pool& streams, } } } + +/// MARK_NEXT_CONSTRAINT + +template +void upd_next_constraint_heavy_vars(stream_pool_t& streams, + bounds_update_view_t view, + const rmm::device_uvector& heavy_vars_vertex_ids, + const rmm::device_uvector& heavy_vars_pseudo_block_ids, + const rmm::device_uvector& heavy_vars_block_segments, + const std::vector& vars_bin_offsets, + i_t heavy_degree_cutoff, + i_t num_blocks_heavy_vars, + bool dry_run = false) +{ + if (num_blocks_heavy_vars != 0) { + auto heavy_vars_stream = streams.get_stream(); + if (!dry_run) { + auto heavy_vars_beg_id = get_id_offset(vars_bin_offsets, heavy_degree_cutoff); + lb_upd_next_constraint_heavy_kernel + <<>>( + heavy_vars_beg_id, + make_span(heavy_vars_vertex_ids), + make_span(heavy_vars_pseudo_block_ids), + heavy_degree_cutoff, + view); + } + } +} + +template +void upd_next_constraint_per_block(stream_pool_t& streams, + bounds_update_view_t view, + const std::vector& vars_bin_offsets, + i_t degree_beg, + i_t degree_end, + bool dry_run) +{ + static_assert(block_dim <= 1024, "Cannot launch kernel with more than 1024 threads"); + + auto [vars_id_beg, vars_id_end] = get_id_range(vars_bin_offsets, degree_beg, degree_end); + + auto block_count = vars_id_end - vars_id_beg; + if (block_count > 0) { + auto block_stream = streams.get_stream(); + if (!dry_run) { + lb_upd_next_constraint_block_kernel + <<>>(vars_id_beg, view); + } + } +} + +template +void upd_next_constraint_per_block(stream_pool_t& streams, + bounds_update_view_t view, + const std::vector& vars_bin_offsets, + i_t heavy_degree_cutoff, + bool dry_run = false) +{ + if (view.nnz < 10000) { + upd_next_constraint_per_block(streams, view, vars_bin_offsets, 32, 32, dry_run); + upd_next_constraint_per_block(streams, view, vars_bin_offsets, 64, 64, dry_run); + upd_next_constraint_per_block( + streams, view, vars_bin_offsets, 128, 128, dry_run); + upd_next_constraint_per_block( + streams, view, vars_bin_offsets, 256, 256, dry_run); + } else { + //[1024, heavy_degree_cutoff/2] -> 128 block size + upd_next_constraint_per_block( + streams, view, vars_bin_offsets, 1024, heavy_degree_cutoff / 2, dry_run); + //[64, 512] -> 32 block size + upd_next_constraint_per_block(streams, view, vars_bin_offsets, 128, 512, dry_run); + } +} + +template +void upd_next_constraint_sub_warp(stream_pool_t& streams, + bounds_update_view_t view, + i_t degree_beg, + i_t degree_end, + const std::vector& vars_bin_offsets, + bool dry_run) +{ + constexpr i_t block_dim = 32; + auto vars_per_block = block_dim / threads_per_variable; + auto [vars_id_beg, vars_id_end] = get_id_range(vars_bin_offsets, degree_beg, degree_end); + + auto block_count = raft::ceildiv(vars_id_end - vars_id_beg, vars_per_block); + if (block_count != 0) { + auto sub_warp_stream = streams.get_stream(); + if (!dry_run) { + lb_upd_next_constraint_sub_warp_kernel + <<>>(vars_id_beg, vars_id_end, view); + } + } +} + +template +void upd_next_constraint_sub_warp(stream_pool_t& streams, + bounds_update_view_t view, + i_t vars_sub_warp_count, + rmm::device_uvector& warp_vars_offsets, + rmm::device_uvector& warp_vars_id_offsets, + bool dry_run) +{ + constexpr i_t block_dim = 256; + + auto block_count = raft::ceildiv(vars_sub_warp_count * 32, block_dim); + if (block_count != 0) { + auto sub_warp_stream = streams.get_stream(); + if (!dry_run) { + lb_upd_next_constraint_sub_warp_kernel + <<>>( + view, make_span(warp_vars_offsets), make_span(warp_vars_id_offsets)); + } + } +} + +template +void upd_next_constraint_sub_warp(stream_pool_t& streams, + bounds_update_view_t view, + i_t degree, + const std::vector& vars_bin_offsets, + bool dry_run) +{ + upd_next_constraint_sub_warp( + streams, view, degree, degree, vars_bin_offsets, dry_run); +} + +template +void upd_next_constraint_sub_warp(stream_pool_t& streams, + bounds_update_view_t view, + bool is_vars_sub_warp_single_bin, + i_t vars_sub_warp_count, + rmm::device_uvector& warp_vars_offsets, + rmm::device_uvector& warp_vars_id_offsets, + const std::vector& vars_bin_offsets, + bool dry_run = false) +{ + if (view.nnz < 10000) { + upd_next_constraint_sub_warp(streams, view, 16, vars_bin_offsets, dry_run); + upd_next_constraint_sub_warp(streams, view, 8, vars_bin_offsets, dry_run); + upd_next_constraint_sub_warp(streams, view, 4, vars_bin_offsets, dry_run); + upd_next_constraint_sub_warp(streams, view, 2, vars_bin_offsets, dry_run); + upd_next_constraint_sub_warp(streams, view, 1, vars_bin_offsets, dry_run); + } else { + if (is_vars_sub_warp_single_bin) { + upd_next_constraint_sub_warp(streams, view, 64, vars_bin_offsets, dry_run); + upd_next_constraint_sub_warp(streams, view, 32, vars_bin_offsets, dry_run); + upd_next_constraint_sub_warp(streams, view, 16, vars_bin_offsets, dry_run); + upd_next_constraint_sub_warp(streams, view, 8, vars_bin_offsets, dry_run); + upd_next_constraint_sub_warp(streams, view, 1, 4, vars_bin_offsets, dry_run); + } else { + upd_next_constraint_sub_warp( + streams, view, vars_sub_warp_count, warp_vars_offsets, warp_vars_id_offsets, dry_run); + } + } +} + } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/load_balanced_bounds_presolve_kernels.cuh b/cpp/src/mip/presolve/load_balanced_bounds_presolve_kernels.cuh index f51444f895..c408e019c3 100644 --- a/cpp/src/mip/presolve/load_balanced_bounds_presolve_kernels.cuh +++ b/cpp/src/mip/presolve/load_balanced_bounds_presolve_kernels.cuh @@ -13,32 +13,18 @@ #pragma once #include +#include +#include "load_balanced_bounds_common_kernels.cuh" namespace cuopt::linear_programming::detail { -template -__device__ __forceinline__ void detect_range_sub_warp(i_t* id_warp_beg, - i_t* id_range_end, - i_t* threads_per_item, - raft::device_span warp_offsets, - raft::device_span bin_offsets) -{ - i_t warp_id = (blockDim.x * blockIdx.x + threadIdx.x) / raft::WarpSize; - i_t lane_id = threadIdx.x & 31; - bool pred = false; - if (lane_id < warp_offsets.size()) { pred = (warp_id >= warp_offsets[lane_id]); } - unsigned int m = __ballot_sync(0xffffffff, pred); - i_t seg = 31 - __clz(m); - i_t it_per_warp = (1 << (5 - seg)); // item per warp = raft::WarpSize/(2^seg) - if (5 - seg < 0) { - *threads_per_item = 0; - return; - } - i_t beg = bin_offsets[seg] + (warp_id - warp_offsets[seg]) * it_per_warp; - i_t end = bin_offsets[seg + 1]; - *id_warp_beg = beg; - *id_range_end = end; - *threads_per_item = (1 << seg); +template +inline __device__ thrust::pair skip_cnst_calc(upd_view_t upd_0, + upd_view_t upd_1, + i_t cnst_idx) +{ + return thrust::make_pair((upd_0.changed_constraints[cnst_idx] == i_t{0}), + (upd_1.changed_constraints[cnst_idx] == i_t{0})); } template +template +__device__ f_t2 calc_act(activity_view_t view, upd_view_t upd, i_t tid, i_t beg, i_t end) +{ + auto act = f_t2{0., 0.}; + for (i_t i = tid + beg; i < end; i += MAX_EDGE_PER_CNST) { + auto coeff = view.coeff[i]; + auto var = view.vars[i]; + + atomicExch(&upd.changed_variables[var], 1); + + auto bounds = upd.vars_bnd[var]; + auto min_contrib = bounds.x; + auto max_contrib = bounds.y; + if (coeff < 0.0) { + min_contrib = bounds.y; + max_contrib = bounds.x; + } + act.x += coeff * min_contrib; + act.y += coeff * max_contrib; + } + return act; +} + +template +__device__ thrust::pair calc_act( + activity_view_t view, upd_view_t upd_0, upd_view_t upd_1, i_t tid, i_t beg, i_t end) +{ + auto act_0 = f_t2{0., 0.}; + auto act_1 = f_t2{0., 0.}; + for (i_t i = tid + beg; i < end; i += MAX_EDGE_PER_CNST) { + auto coeff = view.coeff[i]; + auto var = view.vars[i]; + + atomicExch(&upd_0.changed_variables[var], 1); + atomicExch(&upd_1.changed_variables[var], 1); + + auto bounds_0 = upd_0.vars_bnd[var]; + auto bounds_1 = upd_1.vars_bnd[var]; + auto min_contrib_0 = bounds_0.x; + auto max_contrib_0 = bounds_0.y; + auto min_contrib_1 = bounds_1.x; + auto max_contrib_1 = bounds_1.y; + if (coeff < 0.0) { + min_contrib_0 = bounds_0.y; + max_contrib_0 = bounds_0.x; + min_contrib_1 = bounds_1.y; + max_contrib_1 = bounds_1.x; + } + act_0.x += coeff * min_contrib_0; + act_0.y += coeff * max_contrib_0; + act_1.x += coeff * min_contrib_1; + act_1.y += coeff * max_contrib_1; + } + return thrust::make_pair(act_0, act_1); +} + +template __global__ void lb_calc_act_heavy_kernel(i_t id_range_beg, raft::device_span ids, raft::device_span pseudo_block_ids, i_t work_per_block, activity_view_t view, - raft::device_span tmp_cnst_act) + upd_view_t upd_0, + upd_view_t upd_1) { - auto idx = ids[blockIdx.x] + id_range_beg; + auto idx = ids[blockIdx.x] + id_range_beg; + auto cnst_idx = view.cnst_reorg_ids[idx]; + auto skip_calc = skip_cnst_calc(upd_0, upd_1, cnst_idx); + if (thrust::get<0>(skip_calc) && thrust::get<1>(skip_calc)) { return; } auto pseudo_block_id = pseudo_block_ids[blockIdx.x]; i_t item_off_beg = view.offsets[idx] + work_per_block * pseudo_block_id; i_t item_off_end = min(item_off_beg + work_per_block, view.offsets[idx + 1]); @@ -82,36 +145,58 @@ __global__ void lb_calc_act_heavy_kernel(i_t id_range_beg, typedef cub::BlockReduce BlockReduce; __shared__ typename BlockReduce::TempStorage temp_storage; - auto act = calc_act(view, threadIdx.x, item_off_beg, item_off_end); - - act.x = BlockReduce(temp_storage).Sum(act.x); - __syncthreads(); - act.y = BlockReduce(temp_storage).Sum(act.y); - - // don't subtract constraint bounds yet - // to be done in post processing in finalize_calc_act_kernel - if (threadIdx.x == 0) { tmp_cnst_act[blockIdx.x] = act; } + if (thrust::get<1>(skip_calc)) { + auto act = calc_act(view, upd_0, threadIdx.x, item_off_beg, item_off_end); + + act.x = BlockReduce(temp_storage).Sum(act.x); + __syncthreads(); + act.y = BlockReduce(temp_storage).Sum(act.y); + + // don't subtract constraint bounds yet + // to be done in post processing in finalize_calc_act_kernel + if (threadIdx.x == 0) { upd_0.tmp_cnst_slack[blockIdx.x] = act; } + } else if (thrust::get<0>(skip_calc)) { + auto act = calc_act(view, upd_1, threadIdx.x, item_off_beg, item_off_end); + + act.x = BlockReduce(temp_storage).Sum(act.x); + __syncthreads(); + act.y = BlockReduce(temp_storage).Sum(act.y); + + // don't subtract constraint bounds yet + // to be done in post processing in finalize_calc_act_kernel + if (threadIdx.x == 0) { upd_1.tmp_cnst_slack[blockIdx.x] = act; } + } else { + auto act = + calc_act(view, upd_0, upd_1, threadIdx.x, item_off_beg, item_off_end); + + thrust::get<0>(act).x = BlockReduce(temp_storage).Sum(thrust::get<0>(act).x); + __syncthreads(); + thrust::get<0>(act).y = BlockReduce(temp_storage).Sum(thrust::get<0>(act).y); + __syncthreads(); + thrust::get<1>(act).x = BlockReduce(temp_storage).Sum(thrust::get<1>(act).x); + __syncthreads(); + thrust::get<1>(act).y = BlockReduce(temp_storage).Sum(thrust::get<1>(act).y); + if (threadIdx.x == 0) { + upd_0.tmp_cnst_slack[blockIdx.x] = thrust::get<0>(act); + upd_1.tmp_cnst_slack[blockIdx.x] = thrust::get<1>(act); + } + } } template __global__ void lb_calc_act_heavy_kernel(i_t id_range_beg, - raft::device_span item_block_segments, + raft::device_span ids, + raft::device_span pseudo_block_ids, i_t work_per_block, activity_view_t view, raft::device_span tmp_cnst_act) { - __shared__ i_t id_map; - __shared__ i_t pseudo_block_id; - if (threadIdx.x == 0) { - id_map = thrust::upper_bound( - thrust::seq, item_block_segments.begin(), item_block_segments.end(), blockIdx.x) - - item_block_segments.begin(); - pseudo_block_id = blockIdx.x - item_block_segments[id_map - 1]; - } - __syncthreads(); - auto idx = id_range_beg + id_map; - i_t item_off_beg = view.offsets[idx] + work_per_block * pseudo_block_id; - i_t item_off_end = min(item_off_beg + work_per_block, view.offsets[idx + 1]); + auto idx = ids[blockIdx.x] + id_range_beg; + auto cnst_idx = view.cnst_reorg_ids[idx]; + if (view.changed_constraints[cnst_idx] == 0) { return; } + auto pseudo_block_id = pseudo_block_ids[blockIdx.x]; + i_t item_off_beg = view.offsets[idx] + work_per_block * pseudo_block_id; + i_t item_off_end = min(item_off_beg + work_per_block, view.offsets[idx + 1]); typedef cub::BlockReduce BlockReduce; __shared__ typename BlockReduce::TempStorage temp_storage; @@ -140,6 +225,125 @@ inline __device__ void write_cnst_slack( view.cnst_slack[cnst_idx] = cnst_prop; } +template +inline __device__ void write_cnst_slack(activity_view_t view, + i_t cnst_idx, + f_t2 cnst_lb_ub, + f_t2 act) +{ + auto cnst_prop = f_t2{cnst_lb_ub.y - act.x, cnst_lb_ub.x - act.y}; + view.cnst_slack[cnst_idx] = cnst_prop; +} + +template +inline __device__ void write_cnst_slack( + activity_view_t view, i_t cnst_idx, f_t2 cnst_lb_ub, f_t2 act, f_t eps) +{ + auto cnst_prop = f_t2{cnst_lb_ub.y - act.x, cnst_lb_ub.x - act.y}; + if ((0 > cnst_prop.x + eps) || (eps < cnst_prop.y)) { + cnst_prop.x = std::numeric_limits::quiet_NaN(); + } + view.cnst_slack[cnst_idx] = cnst_prop; +} + +template +__global__ void finalize_calc_act_kernel(i_t heavy_cnst_beg_id, + raft::device_span item_offsets, + activity_view_t view, + upd_view_t upd_0, + upd_view_t upd_1) +{ + using warp_reduce = cub::WarpReduce; + __shared__ typename warp_reduce::TempStorage temp_storage; + i_t idx = heavy_cnst_beg_id + blockIdx.x; + i_t cnst_idx = view.cnst_reorg_ids[idx]; + auto skip_calc = skip_cnst_calc(upd_0, upd_1, cnst_idx); + if (thrust::get<0>(skip_calc) && thrust::get<1>(skip_calc)) { return; } + auto cnst_lb_ub = view.cnst_bnd[idx]; + [[maybe_unused]] f_t eps = {}; + if constexpr (erase_inf_cnst) { + eps = get_cstr_tolerance(cnst_lb_ub.x, + cnst_lb_ub.y, + view.tolerances.absolute_tolerance, + view.tolerances.relative_tolerance); + } + + // assumes cnst_bnd[i].x has ub and cnst_bnd[i].y has lb + i_t item_off_beg = item_offsets[blockIdx.x]; + i_t item_off_end = item_offsets[blockIdx.x + 1]; + + if (thrust::get<1>(skip_calc)) { + f_t2 cnst_prop = f_t2{0., 0.}; + // assumes tmp_act[i].x has min activity and tmp_act[i].y has max activity + for (i_t i = threadIdx.x + item_off_beg; i < item_off_end; i += blockDim.x) { + auto act = upd_0.tmp_cnst_slack[i]; + cnst_prop.x += act.x; + cnst_prop.y += act.y; + } + cnst_prop.x = warp_reduce(temp_storage).Sum(cnst_prop.x); + __syncwarp(); + cnst_prop.y = warp_reduce(temp_storage).Sum(cnst_prop.y); + if (threadIdx.x == 0) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, cnst_prop, eps); + } else { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, cnst_prop); + } + } + } else if (thrust::get<0>(skip_calc)) { + f_t2 cnst_prop = f_t2{0., 0.}; + // assumes tmp_act[i].x has min activity and tmp_act[i].y has max activity + for (i_t i = threadIdx.x + item_off_beg; i < item_off_end; i += blockDim.x) { + auto act = upd_1.tmp_cnst_slack[i]; + cnst_prop.x += act.x; + cnst_prop.y += act.y; + } + cnst_prop.x = warp_reduce(temp_storage).Sum(cnst_prop.x); + __syncwarp(); + cnst_prop.y = warp_reduce(temp_storage).Sum(cnst_prop.y); + if (threadIdx.x == 0) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, cnst_prop, eps); + } else { + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, cnst_prop); + } + } + } else { + f_t2 cnst_prop_0 = f_t2{0., 0.}; + f_t2 cnst_prop_1 = f_t2{0., 0.}; + // assumes tmp_act[i].x has min activity and tmp_act[i].y has max activity + for (i_t i = threadIdx.x + item_off_beg; i < item_off_end; i += blockDim.x) { + auto act = upd_0.tmp_cnst_slack[i]; + cnst_prop_0.x += act.x; + cnst_prop_0.y += act.y; + act = upd_1.tmp_cnst_slack[i]; + cnst_prop_1.x += act.x; + cnst_prop_1.y += act.y; + } + cnst_prop_0.x = warp_reduce(temp_storage).Sum(cnst_prop_0.x); + __syncwarp(); + cnst_prop_0.y = warp_reduce(temp_storage).Sum(cnst_prop_0.y); + __syncwarp(); + cnst_prop_1.x = warp_reduce(temp_storage).Sum(cnst_prop_1.x); + __syncwarp(); + cnst_prop_1.y = warp_reduce(temp_storage).Sum(cnst_prop_1.y); + if (threadIdx.x == 0) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, cnst_prop_0, eps); + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, cnst_prop_1, eps); + } else { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, cnst_prop_0); + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, cnst_prop_1); + } + } + } +} + template __global__ void finalize_calc_act_kernel(i_t heavy_cnst_beg_id, raft::device_span item_offsets, @@ -148,8 +352,9 @@ __global__ void finalize_calc_act_kernel(i_t heavy_cnst_beg_id, { using warp_reduce = cub::WarpReduce; __shared__ typename warp_reduce::TempStorage temp_storage; - i_t idx = heavy_cnst_beg_id + blockIdx.x; - i_t cnst_idx = view.cnst_reorg_ids[idx]; + i_t idx = heavy_cnst_beg_id + blockIdx.x; + i_t cnst_idx = view.cnst_reorg_ids[idx]; + if (view.changed_constraints[cnst_idx] == 0) { return; } auto cnst_lb_ub = view.cnst_bnd[idx]; [[maybe_unused]] f_t eps = {}; if constexpr (erase_inf_cnst) { @@ -177,6 +382,45 @@ __global__ void finalize_calc_act_kernel(i_t heavy_cnst_beg_id, } } +// template +//__global__ void finalize_calc_act_kernel(i_t heavy_cnst_beg_id, +// raft::device_span item_offsets, +// activity_view_t view, +// upd_view_t upd) +//{ +// using warp_reduce = cub::WarpReduce; +// __shared__ typename warp_reduce::TempStorage temp_storage; +// i_t idx = heavy_cnst_beg_id + blockIdx.x; +// i_t cnst_idx = view.cnst_reorg_ids[idx]; +// if (view.changed_constraints[cnst_idx] == 0) { return; } +// auto cnst_lb_ub = view.cnst_bnd[idx]; +// [[maybe_unused]] f_t eps = {}; +// if constexpr (erase_inf_cnst) { +// eps = get_cstr_tolerance(cnst_lb_ub.x, +// cnst_lb_ub.y, +// view.tolerances.absolute_tolerance, +// view.tolerances.relative_tolerance); +// } +// +// // assumes cnst_bnd[i].x has ub and cnst_bnd[i].y has lb +// i_t item_off_beg = item_offsets[blockIdx.x]; +// i_t item_off_end = item_offsets[blockIdx.x + 1]; +// f_t2 cnst_prop = f_t2{0., 0.}; +// // assumes tmp_act[i].x has min activity and tmp_act[i].y has max activity +// for (i_t i = threadIdx.x + item_off_beg; i < item_off_end; i += blockDim.x) { +// auto act = upd.tmp_act[i]; +// cnst_prop.x += act.x; +// cnst_prop.y += act.y; +// } +// cnst_prop.x = warp_reduce(temp_storage).Sum(cnst_prop.x); +// __syncwarp(); +// cnst_prop.y = warp_reduce(temp_storage).Sum(cnst_prop.y); +// if (threadIdx.x == 0) { +// write_cnst_slack(upd, cnst_idx, cnst_lb_ub, cnst_prop, eps); +// } +// } + template (view, cnst_idx, cnst_lb_ub, act, eps); } } +template +__global__ void lb_calc_act_block_kernel(i_t id_range_beg, activity_view_t view, upd_view_t upd) + +{ + i_t idx = id_range_beg + blockIdx.x; + i_t cnst_idx = view.cnst_reorg_ids[idx]; + if (view.changed_constraints[cnst_idx] == 0) { return; } + auto cnst_lb_ub = view.cnst_bnd[idx]; + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + [[maybe_unused]] f_t eps = {}; + if constexpr (erase_inf_cnst) { + eps = get_cstr_tolerance(cnst_lb_ub.x, + cnst_lb_ub.y, + view.tolerances.absolute_tolerance, + view.tolerances.relative_tolerance); + } + + typedef cub::BlockReduce BlockReduce; + __shared__ typename BlockReduce::TempStorage temp_storage; + + auto act = calc_act(view, upd, threadIdx.x, item_off_beg, item_off_end); + + act.x = BlockReduce(temp_storage).Sum(act.x); + __syncthreads(); + act.y = BlockReduce(temp_storage).Sum(act.y); + + if (threadIdx.x == 0) { write_cnst_slack(upd, cnst_idx, cnst_lb_ub, act, eps); } +} + +template +__global__ void lb_calc_act_block_kernel(i_t id_range_beg, + activity_view_t view, + upd_view_t upd_0, + upd_view_t upd_1) + +{ + i_t idx = id_range_beg + blockIdx.x; + i_t cnst_idx = view.cnst_reorg_ids[idx]; + auto skip_calc = skip_cnst_calc(upd_0, upd_1, cnst_idx); + if (thrust::get<0>(skip_calc) && thrust::get<1>(skip_calc)) { return; } + auto cnst_lb_ub = view.cnst_bnd[idx]; + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + [[maybe_unused]] f_t eps = {}; + if constexpr (erase_inf_cnst) { + eps = get_cstr_tolerance(cnst_lb_ub.x, + cnst_lb_ub.y, + view.tolerances.absolute_tolerance, + view.tolerances.relative_tolerance); + } + + typedef cub::BlockReduce BlockReduce; + __shared__ typename BlockReduce::TempStorage temp_storage; + + if (thrust::get<1>(skip_calc)) { + auto act = calc_act(view, upd_0, threadIdx.x, item_off_beg, item_off_end); + + act.x = BlockReduce(temp_storage).Sum(act.x); + __syncthreads(); + act.y = BlockReduce(temp_storage).Sum(act.y); + + if (threadIdx.x == 0) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, act, eps); + } else { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, act); + } + } + } else if (thrust::get<0>(skip_calc)) { + auto act = calc_act(view, upd_1, threadIdx.x, item_off_beg, item_off_end); + + act.x = BlockReduce(temp_storage).Sum(act.x); + __syncthreads(); + act.y = BlockReduce(temp_storage).Sum(act.y); + + if (threadIdx.x == 0) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, act, eps); + } else { + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, act); + } + } + } else { + auto act = + calc_act(view, upd_0, upd_1, threadIdx.x, item_off_beg, item_off_end); + + thrust::get<0>(act).x = BlockReduce(temp_storage).Sum(thrust::get<0>(act).x); + __syncthreads(); + thrust::get<0>(act).y = BlockReduce(temp_storage).Sum(thrust::get<0>(act).y); + __syncthreads(); + thrust::get<1>(act).x = BlockReduce(temp_storage).Sum(thrust::get<1>(act).x); + __syncthreads(); + thrust::get<1>(act).y = BlockReduce(temp_storage).Sum(thrust::get<1>(act).y); + + if (threadIdx.x == 0) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, thrust::get<0>(act), eps); + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, thrust::get<1>(act), eps); + } else { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, thrust::get<0>(act)); + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, thrust::get<1>(act)); + } + } + } +} + template -__global__ void lb_calc_act_sub_warp_kernel_test(i_t id_range_beg, - i_t id_range_end, - activity_view_t view) + typename activity_view_t, + typename upd_view_t> +__global__ void lb_calc_act_sub_warp_kernel( + i_t id_range_beg, i_t id_range_end, activity_view_t view, upd_view_t upd_0, upd_view_t upd_1) { constexpr i_t ids_per_block = BDIM / MAX_EDGE_PER_CNST; i_t id_beg = blockIdx.x * ids_per_block + id_range_beg; i_t idx = id_beg + (threadIdx.x / MAX_EDGE_PER_CNST); + // bool continue_calc = (idx < id_range_end); + auto skip_calc = thrust::make_pair(!(idx < id_range_end), !(idx < id_range_end)); i_t cnst_idx; - f_t eps; + [[maybe_unused]] f_t eps = {}; f_t2 cnst_lb_ub; if (idx < id_range_end) { cnst_idx = view.cnst_reorg_ids[idx]; + skip_calc = skip_cnst_calc(upd_0, upd_1, cnst_idx); cnst_lb_ub = view.cnst_bnd[idx]; if constexpr (erase_inf_cnst) { eps = get_cstr_tolerance(cnst_lb_ub.x, @@ -245,20 +612,48 @@ __global__ void lb_calc_act_sub_warp_kernel_test(i_t id_range_beg, using warp_reduce = cub::WarpReduce; __shared__ typename warp_reduce::TempStorage temp_storage; - auto act = f_t2{0., 0.}; + auto act = thrust::make_pair(f_t2{0., 0.}, f_t2{0., 0.}); - if (idx < id_range_end) { - i_t item_off_beg = view.offsets[idx]; - i_t item_off_end = view.offsets[idx + 1]; - act = calc_act(view, p_tid, item_off_beg, item_off_end); + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + if (thrust::get<1>(skip_calc)) { + thrust::get<0>(act) = + calc_act(view, upd_0, p_tid, item_off_beg, item_off_end); + } else if (thrust::get<0>(skip_calc)) { + thrust::get<1>(act) = + calc_act(view, upd_1, p_tid, item_off_beg, item_off_end); + } else { + act = calc_act( + view, upd_0, upd_1, p_tid, item_off_beg, item_off_end); } - - act.x = warp_reduce(temp_storage).Sum(act.x); + thrust::get<0>(act).x = warp_reduce(temp_storage).Sum(thrust::get<0>(act).x); __syncwarp(); - act.y = warp_reduce(temp_storage).Sum(act.y); + thrust::get<0>(act).y = warp_reduce(temp_storage).Sum(thrust::get<0>(act).y); + __syncwarp(); + thrust::get<1>(act).x = warp_reduce(temp_storage).Sum(thrust::get<1>(act).x); + __syncwarp(); + thrust::get<1>(act).y = warp_reduce(temp_storage).Sum(thrust::get<1>(act).y); - if (head_flag && (idx < id_range_end)) { - write_cnst_slack(view, cnst_idx, cnst_lb_ub, act, eps); + if (head_flag && thrust::get<1>(skip_calc)) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, thrust::get<0>(act), eps); + } else { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, thrust::get<0>(act)); + } + } else if (head_flag && thrust::get<0>(skip_calc)) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, thrust::get<1>(act), eps); + } else { + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, thrust::get<1>(act)); + } + } else if (head_flag) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, thrust::get<0>(act), eps); + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, thrust::get<1>(act), eps); + } else { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, thrust::get<0>(act)); + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, thrust::get<1>(act)); + } } } @@ -276,11 +671,13 @@ __global__ void lb_calc_act_sub_warp_kernel(i_t id_range_beg, constexpr i_t ids_per_block = BDIM / MAX_EDGE_PER_CNST; i_t id_beg = blockIdx.x * ids_per_block + id_range_beg; i_t idx = id_beg + (threadIdx.x / MAX_EDGE_PER_CNST); + bool continue_calc = (idx < id_range_end); i_t cnst_idx; [[maybe_unused]] f_t eps = {}; f_t2 cnst_lb_ub; - if (idx < id_range_end) { - cnst_idx = view.cnst_reorg_ids[idx]; + if (continue_calc) { + cnst_idx = view.cnst_reorg_ids[idx]; + if (view.changed_constraints[cnst_idx] == 0) { continue_calc = false; } cnst_lb_ub = view.cnst_bnd[idx]; if constexpr (erase_inf_cnst) { eps = get_cstr_tolerance(cnst_lb_ub.x, @@ -298,7 +695,7 @@ __global__ void lb_calc_act_sub_warp_kernel(i_t id_range_beg, auto act = f_t2{0., 0.}; - if (idx < id_range_end) { + if (continue_calc) { i_t item_off_beg = view.offsets[idx]; i_t item_off_end = view.offsets[idx + 1]; act = calc_act(view, p_tid, item_off_beg, item_off_end); @@ -308,7 +705,7 @@ __global__ void lb_calc_act_sub_warp_kernel(i_t id_range_beg, __syncwarp(); act.y = warp_reduce(temp_storage).Sum(act.y); - if (head_flag && (idx < id_range_end)) { + if (head_flag && continue_calc) { write_cnst_slack(view, cnst_idx, cnst_lb_ub, act, eps); } } @@ -322,14 +719,16 @@ template __device__ void calc_act_sub_warp(i_t id_warp_beg, i_t id_range_end, activity_view_t view) { - i_t lane_id = (threadIdx.x & 31); - i_t idx = id_warp_beg + (lane_id / MAX_EDGE_PER_CNST); + i_t lane_id = (threadIdx.x & 31); + i_t idx = id_warp_beg + (lane_id / MAX_EDGE_PER_CNST); + bool continue_calc = (idx < id_range_end); i_t cnst_idx; [[maybe_unused]] f_t eps = {}; f_t2 cnst_lb_ub; - if (idx < id_range_end) { + if (continue_calc) { cnst_idx = view.cnst_reorg_ids[idx]; cnst_lb_ub = view.cnst_bnd[idx]; + if (view.changed_constraints[cnst_idx] == 0) { continue_calc = false; } if constexpr (erase_inf_cnst) { eps = get_cstr_tolerance(cnst_lb_ub.x, cnst_lb_ub.y, @@ -346,7 +745,7 @@ __device__ void calc_act_sub_warp(i_t id_warp_beg, i_t id_range_end, activity_vi auto act = f_t2{0., 0.}; - if (idx < id_range_end) { + if (continue_calc) { i_t item_off_beg = view.offsets[idx]; i_t item_off_end = view.offsets[idx + 1]; act = calc_act(view, p_tid, item_off_beg, item_off_end); @@ -356,7 +755,7 @@ __device__ void calc_act_sub_warp(i_t id_warp_beg, i_t id_range_end, activity_vi __syncwarp(); act.y = warp_reduce(temp_storage).Sum(act.y); - if (head_flag && (idx < id_range_end)) { + if (head_flag && continue_calc) { write_cnst_slack(view, cnst_idx, cnst_lb_ub, act, eps); } } @@ -388,66 +787,236 @@ __global__ void lb_calc_act_sub_warp_kernel(activity_view_t view, } } -/// BOUNDS UPDATE - -template -__device__ f_t2 update_bounds(bounds_update_view_t view, i_t tid, i_t beg, i_t end, f_t2 init) + i_t BDIM, + i_t MAX_EDGE_PER_CNST, + typename activity_view_t, + typename upd_view_t> +__device__ void calc_act_sub_warp( + i_t id_warp_beg, i_t id_range_end, activity_view_t view, upd_view_t upd_0, upd_view_t upd_1) { - f_t2 bounds = init; + i_t lane_id = (threadIdx.x & 31); + i_t idx = id_warp_beg + (lane_id / MAX_EDGE_PER_CNST); + auto skip_calc = thrust::make_pair(!(idx < id_range_end), !(idx < id_range_end)); + i_t cnst_idx; + [[maybe_unused]] f_t eps = {}; + f_t2 cnst_lb_ub; + if (idx < id_range_end) { + cnst_idx = view.cnst_reorg_ids[idx]; + cnst_lb_ub = view.cnst_bnd[idx]; + skip_calc = skip_cnst_calc(upd_0, upd_1, cnst_idx); + if constexpr (erase_inf_cnst) { + eps = get_cstr_tolerance(cnst_lb_ub.x, + cnst_lb_ub.y, + view.tolerances.absolute_tolerance, + view.tolerances.relative_tolerance); + } + } + i_t p_tid = lane_id & (MAX_EDGE_PER_CNST - 1); - const auto old_lb = bounds.x; - const auto old_ub = bounds.y; + i_t head_flag = (p_tid == 0); - for (i_t i = tid + beg; i < end; i += MAX_EDGE_PER_VAR) { - auto a = view.coeff[i]; - auto cnst_idx = view.cnst[i]; + using warp_reduce = cub::WarpReduce; + __shared__ typename warp_reduce::TempStorage temp_storage; - // cnst_slack[cnst_idx].x now has cnst_ub - min_a - // cnst_slack[cnst_idx].y now has cnst_lb - max_a - auto cnstr_data = view.cnst_slack[cnst_idx]; - auto cnstr_ub_minus_min_a = cnstr_data.x; - auto cnstr_lb_minus_max_a = cnstr_data.y; - // don't propagate over constraints that are infeasible - if (isnan(cnstr_data.x)) { continue; } + auto act = thrust::make_pair(f_t2{0., 0.}, f_t2{0., 0.}); - f_t min_contrib = old_lb; - f_t max_contrib = old_ub; - if (a < 0.0) { - min_contrib = old_ub; - max_contrib = old_lb; - } + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; - auto delta_min_act = (cnstr_ub_minus_min_a + (a * min_contrib)) / a; - auto delta_max_act = (cnstr_lb_minus_max_a + (a * max_contrib)) / a; + if (thrust::get<1>(skip_calc)) { + thrust::get<0>(act) = + calc_act(view, upd_0, p_tid, item_off_beg, item_off_end); + } else if (thrust::get<0>(skip_calc)) { + thrust::get<1>(act) = + calc_act(view, upd_1, p_tid, item_off_beg, item_off_end); + } else { + act = calc_act( + view, upd_0, upd_1, p_tid, item_off_beg, item_off_end); + } + thrust::get<0>(act).x = warp_reduce(temp_storage).Sum(thrust::get<0>(act).x); + __syncwarp(); + thrust::get<0>(act).y = warp_reduce(temp_storage).Sum(thrust::get<0>(act).y); + __syncwarp(); + thrust::get<1>(act).x = warp_reduce(temp_storage).Sum(thrust::get<1>(act).x); + __syncwarp(); + thrust::get<1>(act).y = warp_reduce(temp_storage).Sum(thrust::get<1>(act).y); - f_t lb_contrib = delta_max_act; - f_t ub_contrib = delta_min_act; - if (a < 0.0) { - lb_contrib = delta_min_act; - ub_contrib = delta_max_act; + if (head_flag && thrust::get<1>(skip_calc)) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, thrust::get<0>(act), eps); + } else { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, thrust::get<0>(act)); + } + } else if (head_flag && thrust::get<0>(skip_calc)) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, thrust::get<1>(act), eps); + } else { + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, thrust::get<1>(act)); + } + } else if (head_flag) { + if constexpr (erase_inf_cnst) { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, thrust::get<0>(act), eps); + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, thrust::get<1>(act), eps); + } else { + write_cnst_slack(upd_0, cnst_idx, cnst_lb_ub, thrust::get<0>(act)); + write_cnst_slack(upd_1, cnst_idx, cnst_lb_ub, thrust::get<1>(act)); } - bounds.x = max(bounds.x, lb_contrib); - bounds.y = min(bounds.y, ub_contrib); } - - return bounds; } -template -inline __device__ void write_updated_bounds( - f_t2* ptr, bool is_int, bounds_update_view_t view, f_t2 bounds, f_t2 old_bounds) +template +__global__ void lb_calc_act_sub_warp_kernel(activity_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + raft::device_span warp_cnst_offsets, + raft::device_span warp_cnst_id_offsets) { - auto threshold = 1e3 * view.tolerances.absolute_tolerance; - if (is_int) { - bounds.x = ceil(bounds.x - view.tolerances.integrality_tolerance); + i_t id_warp_beg, id_range_end, threads_per_constraints; + detect_range_sub_warp( + &id_warp_beg, &id_range_end, &threads_per_constraints, warp_cnst_offsets, warp_cnst_id_offsets); + + if (threads_per_constraints == 1) { + calc_act_sub_warp( + id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_constraints == 2) { + calc_act_sub_warp( + id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_constraints == 4) { + calc_act_sub_warp( + id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_constraints == 8) { + calc_act_sub_warp( + id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_constraints == 16) { + calc_act_sub_warp( + id_warp_beg, id_range_end, view, upd_0, upd_1); + } +} + +/// BOUNDS UPDATE + +template +inline __device__ void update_bounds_per_cnst(f_t2& bounds, + f_t2 old_bounds, + f_t coeff, + i_t cnst_idx, + raft::device_span cnst_slack, + raft::device_span changed_constraints) +{ + // cnst_slack[cnst_idx].x now has cnst_ub - min_a + // cnst_slack[cnst_idx].y now has cnst_lb - max_a + auto cnstr_data = cnst_slack[cnst_idx]; + bool unchanged = (changed_constraints[cnst_idx] == 0); + auto cnstr_ub_minus_min_a = cnstr_data.x; + auto cnstr_lb_minus_max_a = cnstr_data.y; + // don't propagate over constraints that are infeasible + if (unchanged || isnan(cnstr_data.x)) { return; } + + f_t min_contrib = old_bounds.x; + f_t max_contrib = old_bounds.y; + if (coeff < 0.0) { + min_contrib = old_bounds.y; + max_contrib = old_bounds.x; + } + + auto delta_min_act = (cnstr_ub_minus_min_a + (coeff * min_contrib)) / coeff; + auto delta_max_act = (cnstr_lb_minus_max_a + (coeff * max_contrib)) / coeff; + + f_t lb_contrib = delta_max_act; + f_t ub_contrib = delta_min_act; + if (coeff < 0.0) { + lb_contrib = delta_min_act; + ub_contrib = delta_max_act; + } + bounds.x = max(bounds.x, lb_contrib); + bounds.y = min(bounds.y, ub_contrib); +} + +template +__device__ thrust::pair update_bounds(bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + i_t tid, + i_t beg, + i_t end, + thrust::pair old_bounds) +{ + auto bounds = old_bounds; + + for (i_t i = tid + beg; i < end; i += MAX_EDGE_PER_VAR) { + auto a = view.coeff[i]; + auto cnst_idx = view.cnst[i]; + + update_bounds_per_cnst(thrust::get<0>(bounds), + thrust::get<0>(old_bounds), + a, + cnst_idx, + upd_0.cnst_slack, + upd_0.changed_constraints); + update_bounds_per_cnst(thrust::get<1>(bounds), + thrust::get<1>(old_bounds), + a, + cnst_idx, + upd_1.cnst_slack, + upd_1.changed_constraints); + } + + return bounds; +} + +template +__device__ f_t2 update_bounds(bounds_update_view_t view, + raft::device_span cnst_slack, + raft::device_span changed_constraints, + i_t tid, + i_t beg, + i_t end, + f_t2 old_bounds) +{ + f_t2 bounds = old_bounds; + + for (i_t i = tid + beg; i < end; i += MAX_EDGE_PER_VAR) { + auto a = view.coeff[i]; + auto cnst_idx = view.cnst[i]; + + update_bounds_per_cnst(bounds, old_bounds, a, cnst_idx, cnst_slack, changed_constraints); + } + + return bounds; +} + +template +inline __device__ bool write_updated_bounds( + f_t2* ptr, bool is_int, bounds_update_view_t view, f_t2 bounds, f_t2 old_bounds) +{ + bool changed = false; + auto threshold = 1e3 * view.tolerances.absolute_tolerance; + if (is_int) { + bounds.x = ceil(bounds.x - view.tolerances.integrality_tolerance); bounds.y = floor(bounds.y + view.tolerances.integrality_tolerance); } auto lb_updated = (fabs(bounds.x - old_bounds.x) > threshold); auto ub_updated = (fabs(bounds.y - old_bounds.y) > threshold); + if ((bounds.x != old_bounds.x) || (bounds.y != old_bounds.y)) { changed = true; } if (lb_updated) { old_bounds.x = bounds.x; } if (ub_updated) { old_bounds.y = bounds.y; } @@ -455,6 +1024,30 @@ inline __device__ void write_updated_bounds( *ptr = old_bounds; if (lb_updated || ub_updated) { atomicAdd(view.bounds_changed, 1); } + return changed; +} + +template +inline __device__ bool write_updated_bounds( + f_t2* ptr, bool is_int, bounds_update_view_t view, upd_view_t upd, f_t2 bounds, f_t2 old_bounds) +{ + bool changed = false; + auto threshold = 1e3 * view.tolerances.absolute_tolerance; + if (is_int) { + bounds.x = ceil(bounds.x - view.tolerances.integrality_tolerance); + bounds.y = floor(bounds.y + view.tolerances.integrality_tolerance); + } + auto lb_updated = (fabs(bounds.x - old_bounds.x) > threshold); + auto ub_updated = (fabs(bounds.y - old_bounds.y) > threshold); + if ((bounds.x != old_bounds.x) || (bounds.y != old_bounds.y)) { changed = true; } + + if (lb_updated) { old_bounds.x = bounds.x; } + if (ub_updated) { old_bounds.y = bounds.y; } + + *ptr = old_bounds; + + if (lb_updated || ub_updated) { atomicAdd(upd.bounds_changed, 1); } + return changed; } template @@ -463,13 +1056,17 @@ __global__ void lb_upd_bnd_heavy_kernel(i_t id_range_beg, raft::device_span pseudo_block_ids, i_t work_per_block, bounds_update_view_t view, - raft::device_span tmp_bnd) + raft::device_span tmp_vars_bnd) { auto idx = ids[blockIdx.x] + id_range_beg; auto pseudo_block_id = pseudo_block_ids[blockIdx.x]; auto var_idx = view.vars_reorg_ids[idx]; // x is lb, y is ub - auto old_bounds = view.vars_bnd[var_idx]; + auto old_bounds = view.vars_bnd[var_idx]; + if (view.changed_variables[var_idx] == 0) { + tmp_vars_bnd[blockIdx.x] = old_bounds; + return; + } bool is_int = (view.vars_types[idx] == var_t::INTEGER); i_t item_off_beg = view.offsets[idx] + work_per_block * pseudo_block_id; i_t item_off_end = min(item_off_beg + work_per_block, view.offsets[idx + 1]); @@ -480,31 +1077,39 @@ __global__ void lb_upd_bnd_heavy_kernel(i_t id_range_beg, // if it is a set variable then don't propagate the bound // consider continuous vars as set if their bounds cross or equal if (old_bounds.x + view.tolerances.integrality_tolerance >= old_bounds.y) { - tmp_bnd[blockIdx.x] = old_bounds; + tmp_vars_bnd[blockIdx.x] = old_bounds; return; } - auto bounds = - update_bounds(view, threadIdx.x, item_off_beg, item_off_end, old_bounds); + auto bounds = update_bounds(view, + view.cnst_slack, + view.changed_constraints, + threadIdx.x, + item_off_beg, + item_off_end, + old_bounds); bounds.x = BlockReduce(temp_storage).Reduce(bounds.x, cub::Max()); __syncthreads(); bounds.y = BlockReduce(temp_storage).Reduce(bounds.y, cub::Min()); if (threadIdx.x == 0) { - write_updated_bounds(&tmp_bnd[blockIdx.x], is_int, view, bounds, old_bounds); + bool changed = + write_updated_bounds(&tmp_vars_bnd[blockIdx.x], is_int, view, bounds, old_bounds); + atomicExch(&view.var_bounds_changed[var_idx], 1); } } template __global__ void finalize_upd_bnd_kernel(i_t heavy_vars_beg_id, raft::device_span item_offsets, - raft::device_span tmp_bnd, + raft::device_span tmp_vars_bnd, bounds_update_view_t view) { using warp_reduce = cub::WarpReduce; __shared__ typename warp_reduce::TempStorage temp_storage; i_t idx = heavy_vars_beg_id + blockIdx.x; i_t var_idx = view.vars_reorg_ids[idx]; + if (view.changed_variables[var_idx] == 0) { return; } // assumes cnst_bnd[i].x has ub and cnst_bnd[i].y has lb i_t item_off_beg = item_offsets[blockIdx.x]; @@ -512,7 +1117,7 @@ __global__ void finalize_upd_bnd_kernel(i_t heavy_vars_beg_id, f_t2 bounds = f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; // assumes tmp_act[i].x has min activity and tmp_act[i].y has max activity for (i_t i = threadIdx.x + item_off_beg; i < item_off_end; i += blockDim.x) { - auto bnd = tmp_bnd[i]; + auto bnd = tmp_vars_bnd[i]; bounds.x = max(bounds.x, bnd.x); bounds.y = min(bounds.y, bnd.y); } @@ -527,6 +1132,7 @@ __global__ void lb_upd_bnd_block_kernel(i_t id_range_beg, bounds_update_view_t v { i_t idx = id_range_beg + blockIdx.x; i_t var_idx = view.vars_reorg_ids[idx]; + if (view.changed_variables[var_idx] == 0) { return; } // x is lb, y is ub auto old_bounds = view.vars_bnd[var_idx]; bool is_int = (view.vars_types[idx] == var_t::INTEGER); @@ -539,15 +1145,21 @@ __global__ void lb_upd_bnd_block_kernel(i_t id_range_beg, bounds_update_view_t v // if it is a set variable then don't propagate the bound // consider continuous vars as set if their bounds cross or equal if (old_bounds.x + view.tolerances.integrality_tolerance >= old_bounds.y) { return; } - auto bounds = - update_bounds(view, threadIdx.x, item_off_beg, item_off_end, old_bounds); + auto bounds = update_bounds(view, + view.cnst_slack, + view.changed_constraints, + threadIdx.x, + item_off_beg, + item_off_end, + old_bounds); bounds.x = BlockReduce(temp_storage).Reduce(bounds.x, cub::Max()); __syncthreads(); bounds.y = BlockReduce(temp_storage).Reduce(bounds.y, cub::Min()); if (threadIdx.x == 0) { - write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, old_bounds); + bool changed = write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, old_bounds); + if (changed) { view.var_bounds_changed[var_idx] = changed; } } } @@ -565,10 +1177,12 @@ __global__ void lb_upd_bnd_sub_warp_kernel(i_t id_range_beg, i_t id_range_end, a i_t var_idx; auto old_bounds = f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; - auto bounds = old_bounds; - bool is_int = false; - if (idx < id_range_end) { - var_idx = view.vars_reorg_ids[idx]; + auto bounds = old_bounds; + bool is_int = false; + bool continue_calc = (idx < id_range_end); + if (continue_calc) { + var_idx = view.vars_reorg_ids[idx]; + if (view.changed_variables[var_idx] == 0) { continue_calc = false; } old_bounds = view.vars_bnd[var_idx]; is_int = (view.vars_types[idx] == var_t::INTEGER); } @@ -579,7 +1193,7 @@ __global__ void lb_upd_bnd_sub_warp_kernel(i_t id_range_beg, i_t id_range_end, a using warp_reduce = cub::WarpReduce; __shared__ typename warp_reduce::TempStorage temp_storage; - if (idx < id_range_end) { + if (continue_calc) { // if it is a set variable then don't propagate the bound // consider continuous vars as set if their bounds cross or equal if (old_bounds.x + view.tolerances.integrality_tolerance >= old_bounds.y) { head_flag = 0; } @@ -587,16 +1201,22 @@ __global__ void lb_upd_bnd_sub_warp_kernel(i_t id_range_beg, i_t id_range_end, a i_t item_off_beg = view.offsets[idx]; i_t item_off_end = view.offsets[idx + 1]; - bounds = update_bounds( - view, p_tid, item_off_beg, item_off_end, old_bounds); + bounds = update_bounds(view, + view.cnst_slack, + view.changed_constraints, + p_tid, + item_off_beg, + item_off_end, + old_bounds); } bounds.x = warp_reduce(temp_storage).Reduce(bounds.x, cub::Max()); __syncwarp(); bounds.y = warp_reduce(temp_storage).Reduce(bounds.y, cub::Min()); - if (head_flag && (idx < id_range_end)) { - write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, old_bounds); + if (head_flag && continue_calc) { + bool changed = write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, old_bounds); + if (changed) { view.var_bounds_changed[var_idx] = changed; } } } @@ -613,10 +1233,12 @@ __device__ void upd_bnd_sub_warp(i_t id_warp_beg, i_t id_range_end, bounds_updat i_t var_idx; auto old_bounds = f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; - auto bounds = old_bounds; - bool is_int = false; - if (idx < id_range_end) { - var_idx = view.vars_reorg_ids[idx]; + auto bounds = old_bounds; + bool is_int = false; + bool continue_calc = (idx < id_range_end); + if (continue_calc) { + var_idx = view.vars_reorg_ids[idx]; + if (view.changed_variables[var_idx] == 0) { continue_calc = false; } old_bounds = view.vars_bnd[var_idx]; is_int = (view.vars_types[idx] == var_t::INTEGER); } @@ -629,7 +1251,7 @@ __device__ void upd_bnd_sub_warp(i_t id_warp_beg, i_t id_range_end, bounds_updat using warp_reduce = cub::WarpReduce; __shared__ typename warp_reduce::TempStorage temp_storage; - if (idx < id_range_end) { + if (continue_calc) { // if it is a set variable then don't propagate the bound // consider continuous vars as set if their bounds cross or equal if (old_bounds.x + view.tolerances.integrality_tolerance >= old_bounds.y) { head_flag = 0; } @@ -637,19 +1259,595 @@ __device__ void upd_bnd_sub_warp(i_t id_warp_beg, i_t id_range_end, bounds_updat i_t item_off_beg = view.offsets[idx]; i_t item_off_end = view.offsets[idx + 1]; - bounds = update_bounds( - view, p_tid, item_off_beg, item_off_end, old_bounds); + bounds = update_bounds(view, + view.cnst_slack, + view.changed_constraints, + p_tid, + item_off_beg, + item_off_end, + old_bounds); } bounds.x = warp_reduce(temp_storage).Reduce(bounds.x, cub::Max()); __syncwarp(); bounds.y = warp_reduce(temp_storage).Reduce(bounds.y, cub::Min()); - if (head_flag && (idx < id_range_end)) { - write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, old_bounds); + if (head_flag && continue_calc) { + bool changed = write_updated_bounds(&view.vars_bnd[var_idx], is_int, view, bounds, old_bounds); + if (changed) { view.var_bounds_changed[var_idx] = changed; } + } +} + +#if 1 +template +inline __device__ bool skip_update(f_t2 bnd, f_t int_tol) +{ + return (bnd.x + int_tol >= bnd.y); +} + +template +inline __device__ thrust::pair skip_update( + thrust::pair bnd, upd_view_t upd_0, upd_view_t upd_1, i_t var_idx, f_t int_tol) +{ + return thrust::make_pair((thrust::get<0>(bnd).x + int_tol >= thrust::get<0>(bnd).y) || + (upd_0.changed_variables[var_idx] == 0), + + (thrust::get<1>(bnd).x + int_tol >= thrust::get<1>(bnd).y) || + (upd_1.changed_variables[var_idx] == 0)); +} + +template +inline __device__ thrust::pair skip_update(upd_view_t upd_0, + upd_view_t upd_1, + i_t var_idx) +{ + return thrust::make_pair((upd_0.var_bounds_changed[var_idx] == 0), + (upd_1.var_bounds_changed[var_idx] == 0)); +} + +template +__global__ void lb_upd_bnd_heavy_kernel(i_t id_range_beg, + raft::device_span ids, + raft::device_span pseudo_block_ids, + i_t work_per_block, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1) +{ + auto idx = ids[blockIdx.x] + id_range_beg; + auto pseudo_block_id = pseudo_block_ids[blockIdx.x]; + auto var_idx = view.vars_reorg_ids[idx]; + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + auto old_bounds = thrust::make_pair(upd_0.vars_bnd[var_idx], upd_1.vars_bnd[var_idx]); + // auto skip_calc = thrust::make_pair( + // skip_update(thrust::get<0>(old_bounds), view.tolerances.integrality_tolerance) || + // (upd_0.changed_variables[var_idx] == 0), skip_update(thrust::get<1>(old_bounds), + // view.tolerances.integrality_tolerance) || (upd_1.changed_variables[var_idx] == 0)); + auto skip_calc = + skip_update(old_bounds, upd_0, upd_1, var_idx, view.tolerances.integrality_tolerance); + + typedef cub::BlockReduce BlockReduce; + __shared__ typename BlockReduce::TempStorage temp_storage; + + if (thrust::get<0>(skip_calc) && thrust::get<1>(skip_calc)) { + if (threadIdx.x == 0) { + upd_0.tmp_vars_bnd[blockIdx.x] = thrust::get<0>(old_bounds); + upd_1.tmp_vars_bnd[blockIdx.x] = thrust::get<1>(old_bounds); + } + return; + } else if (thrust::get<0>(skip_calc)) { + if (threadIdx.x == 0) { upd_0.tmp_vars_bnd[blockIdx.x] = thrust::get<0>(old_bounds); } + bool is_int = (view.vars_types[idx] == var_t::INTEGER); + i_t item_off_beg = view.offsets[idx] + work_per_block * pseudo_block_id; + i_t item_off_end = min(item_off_beg + work_per_block, view.offsets[idx + 1]); + auto bounds = update_bounds(view, + upd_1.cnst_slack, + upd_1.changed_constraints, + threadIdx.x, + item_off_beg, + item_off_end, + thrust::get<1>(old_bounds)); + bounds.x = BlockReduce(temp_storage).Reduce(bounds.x, cub::Max()); + __syncthreads(); + bounds.y = BlockReduce(temp_storage).Reduce(bounds.y, cub::Min()); + if (threadIdx.x == 0) { + bool changed = write_updated_bounds( + &upd_1.tmp_vars_bnd[blockIdx.x], is_int, view, upd_1, bounds, thrust::get<1>(old_bounds)); + atomicExch(&upd_1.var_bounds_changed[var_idx], 1); + } + } else if (thrust::get<1>(skip_calc)) { + if (threadIdx.x == 0) { upd_1.tmp_vars_bnd[blockIdx.x] = thrust::get<1>(old_bounds); } + bool is_int = (view.vars_types[idx] == var_t::INTEGER); + i_t item_off_beg = view.offsets[idx] + work_per_block * pseudo_block_id; + i_t item_off_end = min(item_off_beg + work_per_block, view.offsets[idx + 1]); + auto bounds = update_bounds(view, + upd_0.cnst_slack, + upd_0.changed_constraints, + threadIdx.x, + item_off_beg, + item_off_end, + thrust::get<0>(old_bounds)); + bounds.x = BlockReduce(temp_storage).Reduce(bounds.x, cub::Max()); + __syncthreads(); + bounds.y = BlockReduce(temp_storage).Reduce(bounds.y, cub::Min()); + if (threadIdx.x == 0) { + bool changed = write_updated_bounds( + &upd_0.tmp_vars_bnd[blockIdx.x], is_int, view, upd_0, bounds, thrust::get<0>(old_bounds)); + atomicExch(&upd_0.var_bounds_changed[var_idx], 1); + } + } else { + bool is_int = (view.vars_types[idx] == var_t::INTEGER); + i_t item_off_beg = view.offsets[idx] + work_per_block * pseudo_block_id; + i_t item_off_end = min(item_off_beg + work_per_block, view.offsets[idx + 1]); + auto bounds = update_bounds( + view, upd_0, upd_1, threadIdx.x, item_off_beg, item_off_end, old_bounds); + thrust::get<0>(bounds).x = + BlockReduce(temp_storage).Reduce(thrust::get<0>(bounds).x, cub::Max()); + __syncthreads(); + thrust::get<0>(bounds).y = + BlockReduce(temp_storage).Reduce(thrust::get<0>(bounds).y, cub::Min()); + __syncthreads(); + thrust::get<1>(bounds).x = + BlockReduce(temp_storage).Reduce(thrust::get<1>(bounds).x, cub::Max()); + __syncthreads(); + thrust::get<1>(bounds).y = + BlockReduce(temp_storage).Reduce(thrust::get<1>(bounds).y, cub::Min()); + if (threadIdx.x == 0) { + bool changed = write_updated_bounds(&upd_0.tmp_vars_bnd[blockIdx.x], + is_int, + view, + upd_0, + thrust::get<0>(bounds), + thrust::get<0>(old_bounds)); + atomicExch(&upd_0.var_bounds_changed[var_idx], 1); + changed = write_updated_bounds(&upd_1.tmp_vars_bnd[blockIdx.x], + is_int, + view, + upd_1, + thrust::get<1>(bounds), + thrust::get<1>(old_bounds)); + atomicExch(&upd_1.var_bounds_changed[var_idx], 1); + } + } +} + +template +__global__ void finalize_upd_bnd_kernel(i_t heavy_vars_beg_id, + raft::device_span item_offsets, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1) +{ + i_t idx = heavy_vars_beg_id + blockIdx.x; + i_t var_idx = view.vars_reorg_ids[idx]; + + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + // auto skip_calc = thrust::make_pair( + // (upd_0.var_bounds_changed[var_idx] == 0) || (upd_0.changed_variables[var_idx] == 0), + // (upd_1.var_bounds_changed[var_idx] == 0) || (upd_1.changed_variables[var_idx] == 0)); + auto skip_calc = skip_update(upd_0, upd_1, var_idx); + + using warp_reduce = cub::WarpReduce; + __shared__ typename warp_reduce::TempStorage temp_storage; + if (thrust::get<0>(skip_calc) && thrust::get<1>(skip_calc)) { + return; + } else if (thrust::get<0>(skip_calc)) { + // assumes cnst_bnd[i].x has ub and cnst_bnd[i].y has lb + i_t item_off_beg = item_offsets[blockIdx.x]; + i_t item_off_end = item_offsets[blockIdx.x + 1]; + f_t2 bounds = f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + // assumes tmp_act[i].x has min activity and tmp_act[i].y has max activity + for (i_t i = threadIdx.x + item_off_beg; i < item_off_end; i += blockDim.x) { + auto bnd = upd_1.tmp_vars_bnd[i]; + bounds.x = max(bounds.x, bnd.x); + bounds.y = min(bounds.y, bnd.y); + } + bounds.x = warp_reduce(temp_storage).Reduce(bounds.x, cub::Max()); + __syncwarp(); + bounds.y = warp_reduce(temp_storage).Reduce(bounds.y, cub::Min()); + if (threadIdx.x == 0) { upd_1.vars_bnd[var_idx] = bounds; } + } else if (thrust::get<1>(skip_calc)) { + // assumes cnst_bnd[i].x has ub and cnst_bnd[i].y has lb + i_t item_off_beg = item_offsets[blockIdx.x]; + i_t item_off_end = item_offsets[blockIdx.x + 1]; + f_t2 bounds = f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + // assumes tmp_act[i].x has min activity and tmp_act[i].y has max activity + for (i_t i = threadIdx.x + item_off_beg; i < item_off_end; i += blockDim.x) { + auto bnd = upd_0.tmp_vars_bnd[i]; + bounds.x = max(bounds.x, bnd.x); + bounds.y = min(bounds.y, bnd.y); + } + bounds.x = warp_reduce(temp_storage).Reduce(bounds.x, cub::Max()); + __syncwarp(); + bounds.y = warp_reduce(temp_storage).Reduce(bounds.y, cub::Min()); + if (threadIdx.x == 0) { upd_0.vars_bnd[var_idx] = bounds; } + } else { + // assumes cnst_bnd[i].x has ub and cnst_bnd[i].y has lb + i_t item_off_beg = item_offsets[blockIdx.x]; + i_t item_off_end = item_offsets[blockIdx.x + 1]; + f_t2 bounds_0 = + f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + f_t2 bounds_1 = + f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + // assumes tmp_act[i].x has min activity and tmp_act[i].y has max activity + for (i_t i = threadIdx.x + item_off_beg; i < item_off_end; i += blockDim.x) { + auto bnd_0 = upd_0.tmp_vars_bnd[i]; + bounds_0.x = max(bounds_0.x, bnd_0.x); + bounds_0.y = min(bounds_0.y, bnd_0.y); + auto bnd_1 = upd_1.tmp_vars_bnd[i]; + bounds_1.x = max(bounds_1.x, bnd_1.x); + bounds_1.y = min(bounds_1.y, bnd_1.y); + } + bounds_0.x = warp_reduce(temp_storage).Reduce(bounds_0.x, cub::Max()); + __syncwarp(); + bounds_0.y = warp_reduce(temp_storage).Reduce(bounds_0.y, cub::Min()); + __syncwarp(); + bounds_1.x = warp_reduce(temp_storage).Reduce(bounds_1.x, cub::Max()); + __syncwarp(); + bounds_1.y = warp_reduce(temp_storage).Reduce(bounds_1.y, cub::Min()); + if (threadIdx.x == 0) { + upd_0.vars_bnd[var_idx] = bounds_0; + upd_1.vars_bnd[var_idx] = bounds_1; + } + } +} + +template +__device__ void upd_bnd_sub_warp( + i_t id_warp_beg, i_t id_range_end, bounds_update_view_t view, upd_view_t upd_0, upd_view_t upd_1) +{ + i_t lane_id = (threadIdx.x & 31); + i_t idx = id_warp_beg + (lane_id / MAX_EDGE_PER_VAR); + i_t var_idx; + auto old_bounds = thrust::make_pair( + f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}); + auto bounds = old_bounds; + bool is_int = false; + + auto skip_calc = thrust::make_pair(!(idx < id_range_end), !(idx < id_range_end)); + + if (idx < id_range_end) { + var_idx = view.vars_reorg_ids[idx]; + thrust::get<0>(old_bounds) = upd_0.vars_bnd[var_idx]; + thrust::get<1>(old_bounds) = upd_1.vars_bnd[var_idx]; + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + // skip unchanged variables + skip_calc = + skip_update(old_bounds, upd_0, upd_1, var_idx, view.tolerances.integrality_tolerance); + is_int = (view.vars_types[idx] == var_t::INTEGER); + } + // Equivalent to + // i_t p_tid = threadIdx.x % MAX_EDGE_PER_VAR; + i_t p_tid = lane_id & (MAX_EDGE_PER_VAR - 1); + + i_t head_flag = (p_tid == 0); + + using warp_reduce = cub::WarpReduce; + __shared__ typename warp_reduce::TempStorage temp_storage; + + if (thrust::get<0>(skip_calc) && thrust::get<1>(skip_calc)) { + return; + } else if (thrust::get<0>(skip_calc)) { + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + + thrust::get<1>(bounds) = + update_bounds(view, + upd_1.cnst_slack, + upd_1.changed_constraints, + p_tid, + item_off_beg, + item_off_end, + thrust::get<1>(old_bounds)); + } else if (thrust::get<1>(skip_calc)) { + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + + thrust::get<0>(bounds) = + update_bounds(view, + upd_0.cnst_slack, + upd_0.changed_constraints, + p_tid, + item_off_beg, + item_off_end, + thrust::get<0>(old_bounds)); + } else { + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + + bounds = update_bounds( + view, upd_0, upd_1, p_tid, item_off_beg, item_off_end, old_bounds); + } + + thrust::get<0>(bounds).x = warp_reduce(temp_storage).Reduce(thrust::get<0>(bounds).x, cub::Max()); + __syncwarp(); + thrust::get<0>(bounds).y = warp_reduce(temp_storage).Reduce(thrust::get<0>(bounds).y, cub::Min()); + __syncwarp(); + thrust::get<1>(bounds).x = warp_reduce(temp_storage).Reduce(thrust::get<1>(bounds).x, cub::Max()); + __syncwarp(); + thrust::get<1>(bounds).y = warp_reduce(temp_storage).Reduce(thrust::get<1>(bounds).y, cub::Min()); + + if (head_flag && !thrust::get<0>(skip_calc)) { + bool changed = write_updated_bounds(&upd_0.vars_bnd[var_idx], + is_int, + view, + upd_0, + thrust::get<0>(bounds), + thrust::get<0>(old_bounds)); + if (changed) { upd_0.var_bounds_changed[var_idx] = changed; } + } + if (head_flag && !thrust::get<1>(skip_calc)) { + bool changed = write_updated_bounds(&upd_1.vars_bnd[var_idx], + is_int, + view, + upd_1, + thrust::get<1>(bounds), + thrust::get<1>(old_bounds)); + if (changed) { upd_1.var_bounds_changed[var_idx] = changed; } + } +} + +template +__global__ void lb_upd_bnd_sub_warp_kernel(bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1, + raft::device_span warp_vars_offsets, + raft::device_span warp_vars_id_offsets) +{ + i_t id_warp_beg, id_range_end, threads_per_variable; + detect_range_sub_warp( + &id_warp_beg, &id_range_end, &threads_per_variable, warp_vars_offsets, warp_vars_id_offsets); + + if (threads_per_variable == 1) { + upd_bnd_sub_warp(id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_variable == 2) { + upd_bnd_sub_warp(id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_variable == 4) { + upd_bnd_sub_warp(id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_variable == 8) { + upd_bnd_sub_warp(id_warp_beg, id_range_end, view, upd_0, upd_1); + } else if (threads_per_variable == 16) { + upd_bnd_sub_warp(id_warp_beg, id_range_end, view, upd_0, upd_1); } } +template +__global__ void lb_upd_bnd_sub_warp_kernel( + i_t id_range_beg, i_t id_range_end, bounds_update_view_t view, upd_view_t upd_0, upd_view_t upd_1) +{ + constexpr i_t ids_per_block = BDIM / MAX_EDGE_PER_VAR; + i_t id_beg = blockIdx.x * ids_per_block + id_range_beg; + i_t idx = id_beg + (threadIdx.x / MAX_EDGE_PER_VAR); + i_t var_idx; + auto old_bounds = thrust::make_pair( + f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + f_t2{-std::numeric_limits::infinity(), std::numeric_limits::infinity()}); + auto bounds = old_bounds; + bool is_int = false; + + auto skip_calc = thrust::make_pair(!(idx < id_range_end), !(idx < id_range_end)); + + if (idx < id_range_end) { + var_idx = view.vars_reorg_ids[idx]; + thrust::get<0>(old_bounds) = upd_0.vars_bnd[var_idx]; + thrust::get<1>(old_bounds) = upd_1.vars_bnd[var_idx]; + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + // skip unchanged variables + skip_calc = + skip_update(old_bounds, upd_0, upd_1, var_idx, view.tolerances.integrality_tolerance); + is_int = (view.vars_types[idx] == var_t::INTEGER); + } + i_t p_tid = threadIdx.x % MAX_EDGE_PER_VAR; + + i_t head_flag = (p_tid == 0); + + using warp_reduce = cub::WarpReduce; + __shared__ typename warp_reduce::TempStorage temp_storage; + + if (thrust::get<0>(skip_calc) && thrust::get<1>(skip_calc)) { + return; + } else if (thrust::get<0>(skip_calc)) { + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + + thrust::get<1>(bounds) = + update_bounds(view, + upd_1.cnst_slack, + upd_1.changed_constraints, + p_tid, + item_off_beg, + item_off_end, + thrust::get<1>(old_bounds)); + } else if (thrust::get<1>(skip_calc)) { + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + + thrust::get<0>(bounds) = + update_bounds(view, + upd_0.cnst_slack, + upd_0.changed_constraints, + p_tid, + item_off_beg, + item_off_end, + thrust::get<0>(old_bounds)); + } else { + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + + bounds = update_bounds( + view, upd_0, upd_1, p_tid, item_off_beg, item_off_end, old_bounds); + } + + thrust::get<0>(bounds).x = warp_reduce(temp_storage).Reduce(thrust::get<0>(bounds).x, cub::Max()); + __syncwarp(); + thrust::get<0>(bounds).y = warp_reduce(temp_storage).Reduce(thrust::get<0>(bounds).y, cub::Min()); + __syncwarp(); + thrust::get<1>(bounds).x = warp_reduce(temp_storage).Reduce(thrust::get<1>(bounds).x, cub::Max()); + __syncwarp(); + thrust::get<1>(bounds).y = warp_reduce(temp_storage).Reduce(thrust::get<1>(bounds).y, cub::Min()); + + if (head_flag && !thrust::get<0>(skip_calc)) { + bool changed = write_updated_bounds(&upd_0.vars_bnd[var_idx], + is_int, + view, + upd_0, + thrust::get<0>(bounds), + thrust::get<0>(old_bounds)); + if (changed) { upd_0.var_bounds_changed[var_idx] = changed; } + } + if (head_flag && !thrust::get<1>(skip_calc)) { + bool changed = write_updated_bounds(&upd_1.vars_bnd[var_idx], + is_int, + view, + upd_1, + thrust::get<1>(bounds), + thrust::get<1>(old_bounds)); + if (changed) { upd_1.var_bounds_changed[var_idx] = changed; } + } +} + +template +__global__ void lb_upd_bnd_block_kernel(i_t id_range_beg, + bounds_update_view_t view, + upd_view_t upd_0, + upd_view_t upd_1) +{ + i_t idx = id_range_beg + blockIdx.x; + i_t var_idx = view.vars_reorg_ids[idx]; + auto old_bounds = thrust::make_pair(upd_0.vars_bnd[var_idx], upd_1.vars_bnd[var_idx]); + + // if it is a set variable then don't propagate the bound + // consider continuous vars as set if their bounds cross or equal + // auto skip_calc = thrust::make_pair( + // skip_update(thrust::get<0>(old_bounds), view.tolerances.integrality_tolerance) || + // (upd_0.changed_variables[var_idx] == 0), skip_update(thrust::get<1>(old_bounds), + // view.tolerances.integrality_tolerance) || (upd_1.changed_variables[var_idx] == 0)); + auto skip_calc = + skip_update(old_bounds, upd_0, upd_1, var_idx, view.tolerances.integrality_tolerance); + + typedef cub::BlockReduce BlockReduce; + __shared__ typename BlockReduce::TempStorage temp_storage; + + if (thrust::get<0>(skip_calc) && thrust::get<1>(skip_calc)) { + return; + } else if (thrust::get<0>(skip_calc)) { + bool is_int = (view.vars_types[idx] == var_t::INTEGER); + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + auto bounds = update_bounds(view, + upd_1.cnst_slack, + upd_1.changed_constraints, + threadIdx.x, + item_off_beg, + item_off_end, + thrust::get<1>(old_bounds)); + + bounds.x = BlockReduce(temp_storage).Reduce(bounds.x, cub::Max()); + __syncthreads(); + bounds.y = BlockReduce(temp_storage).Reduce(bounds.y, cub::Min()); + + if (threadIdx.x == 0) { + bool changed = write_updated_bounds( + &upd_1.vars_bnd[var_idx], is_int, view, upd_1, bounds, thrust::get<1>(old_bounds)); + if (changed) { upd_1.var_bounds_changed[var_idx] = changed; } + } + } else if (thrust::get<1>(skip_calc)) { + bool is_int = (view.vars_types[idx] == var_t::INTEGER); + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + auto bounds = update_bounds(view, + upd_0.cnst_slack, + upd_0.changed_constraints, + threadIdx.x, + item_off_beg, + item_off_end, + thrust::get<0>(old_bounds)); + + bounds.x = BlockReduce(temp_storage).Reduce(bounds.x, cub::Max()); + __syncthreads(); + bounds.y = BlockReduce(temp_storage).Reduce(bounds.y, cub::Min()); + + if (threadIdx.x == 0) { + bool changed = write_updated_bounds( + &upd_0.vars_bnd[var_idx], is_int, view, upd_0, bounds, thrust::get<0>(old_bounds)); + if (changed) { upd_0.var_bounds_changed[var_idx] = changed; } + } + } else { + bool is_int = (view.vars_types[idx] == var_t::INTEGER); + i_t item_off_beg = view.offsets[idx]; + i_t item_off_end = view.offsets[idx + 1]; + auto bounds = update_bounds( + view, upd_0, upd_1, threadIdx.x, item_off_beg, item_off_end, old_bounds); + + thrust::get<0>(bounds).x = + BlockReduce(temp_storage).Reduce(thrust::get<0>(bounds).x, cub::Max()); + __syncthreads(); + thrust::get<0>(bounds).y = + BlockReduce(temp_storage).Reduce(thrust::get<0>(bounds).y, cub::Min()); + __syncthreads(); + thrust::get<1>(bounds).x = + BlockReduce(temp_storage).Reduce(thrust::get<1>(bounds).x, cub::Max()); + __syncthreads(); + thrust::get<1>(bounds).y = + BlockReduce(temp_storage).Reduce(thrust::get<1>(bounds).y, cub::Min()); + + if (threadIdx.x == 0) { + bool changed = write_updated_bounds(&upd_0.vars_bnd[var_idx], + is_int, + view, + upd_0, + thrust::get<0>(bounds), + thrust::get<0>(old_bounds)); + if (changed) { upd_0.var_bounds_changed[var_idx] = changed; } + changed = write_updated_bounds(&upd_1.vars_bnd[var_idx], + is_int, + view, + upd_1, + thrust::get<1>(bounds), + thrust::get<1>(old_bounds)); + if (changed) { upd_1.var_bounds_changed[var_idx] = changed; } + } + } +} + +#endif + template __global__ void lb_upd_bnd_sub_warp_kernel(bounds_update_view_t view, raft::device_span warp_vars_offsets, diff --git a/cpp/src/mip/presolve/probing_cache.cu b/cpp/src/mip/presolve/probing_cache.cu index cb7bfdae2c..bf0cb0c51f 100644 --- a/cpp/src/mip/presolve/probing_cache.cu +++ b/cpp/src/mip/presolve/probing_cache.cu @@ -18,11 +18,13 @@ #include "probing_cache.cuh" #include -#include +#include #include #include #include +#include +#include #include #include @@ -135,6 +137,49 @@ bool probing_cache_t::contains(i_t var_id) return probing_cache.count(var_id) > 0; } +template +void inline insert_current_probing_to_cache(i_t var_idx, + const val_interval_t& probe_val, + bound_presolve_t& bound_presolve, + const std::vector& original_lb, + const std::vector& original_ub, + const std::vector& modified_bounds, + const std::vector& h_integer_indices, + std::atomic& n_implied_singletons) +{ + f_t int_tol = bound_presolve.context.settings.get_integrality_tolerance(); + + cache_entry_t cache_item; + cache_item.val_interval = probe_val; + for (auto impacted_var_idx : h_integer_indices) { + if (original_lb[impacted_var_idx] != modified_bounds[2 * impacted_var_idx] || + original_ub[impacted_var_idx] != modified_bounds[2 * impacted_var_idx + 1]) { + if (integer_equal(modified_bounds[2 * impacted_var_idx], + modified_bounds[2 * impacted_var_idx + 1], + int_tol)) { + ++n_implied_singletons; + } + cuopt_assert(modified_bounds[2 * impacted_var_idx] >= original_lb[impacted_var_idx], + "Lower bound must be greater than or equal to original lower bound"); + cuopt_assert(modified_bounds[2 * impacted_var_idx + 1] <= original_ub[impacted_var_idx], + "Upper bound must be less than or equal to original upper bound"); + cached_bound_t new_bound{modified_bounds[2 * impacted_var_idx], + modified_bounds[2 * impacted_var_idx + 1]}; + cache_item.var_to_cached_bound_map.insert({impacted_var_idx, new_bound}); + } + } + { + std::lock_guard lock(bound_presolve.probing_cache.probing_cache_mutex); + if (!bound_presolve.probing_cache.probing_cache.count(var_idx) > 0) { + std::array, 2> entries_per_var; + entries_per_var[0] = cache_item; + bound_presolve.probing_cache.probing_cache.insert({var_idx, entries_per_var}); + } else { + bound_presolve.probing_cache.probing_cache[var_idx][1] = cache_item; + } + } +} + template void inline insert_current_probing_to_cache(i_t var_idx, const val_interval_t& probe_val, @@ -344,7 +389,7 @@ template void compute_cache_for_var(i_t var_idx, bound_presolve_t& bound_presolve, problem_t& problem, - multi_probe_t& multi_probe_presolve, + lb_multi_probe_t& multi_probe_presolve, const std::vector& h_var_lower_bounds, const std::vector& h_var_upper_bounds, const std::vector& h_integer_indices, @@ -355,8 +400,7 @@ void compute_cache_for_var(i_t var_idx, RAFT_CUDA_TRY(cudaSetDevice(device_id)); // test if we need per thread handle raft::handle_t handle{}; - std::vector h_improved_lower_bounds(h_var_lower_bounds.size()); - std::vector h_improved_upper_bounds(h_var_upper_bounds.size()); + std::vector h_improved_bounds(2 * h_var_lower_bounds.size()); std::pair, val_interval_t> probe_vals; f_t lb = h_var_lower_bounds[var_idx]; f_t ub = h_var_upper_bounds[var_idx]; @@ -407,8 +451,7 @@ void compute_cache_for_var(i_t var_idx, } } } - auto bounds_presolve_result = - multi_probe_presolve.solve_for_interval(problem, var_interval_vals, &handle); + auto bounds_presolve_result = multi_probe_presolve.solve_for_interval(var_interval_vals, &handle); if (bounds_presolve_result != termination_criterion_t::NO_UPDATE) { CUOPT_LOG_TRACE("Adding cached bounds for var %d", var_idx); } @@ -418,23 +461,16 @@ void compute_cache_for_var(i_t var_idx, // save the impacted bounds if (bounds_presolve_result != termination_criterion_t::NO_UPDATE) { const auto& probe_val = i == 0 ? probe_vals.first : probe_vals.second; - auto& d_lb = i == 0 ? multi_probe_presolve.upd_0.lb : multi_probe_presolve.upd_1.lb; - auto& d_ub = i == 0 ? multi_probe_presolve.upd_0.ub : multi_probe_presolve.upd_1.ub; - raft::copy(h_improved_lower_bounds.data(), - d_lb.data(), - h_improved_lower_bounds.size(), - handle.get_stream()); - raft::copy(h_improved_upper_bounds.data(), - d_ub.data(), - h_improved_upper_bounds.size(), - handle.get_stream()); + auto& d_bnds = + i == 0 ? multi_probe_presolve.upd_0.vars_bnd : multi_probe_presolve.upd_1.vars_bnd; + raft::copy( + h_improved_bounds.data(), d_bnds.data(), h_improved_bounds.size(), handle.get_stream()); insert_current_probing_to_cache(var_idx, probe_val, bound_presolve, h_var_lower_bounds, h_var_upper_bounds, - h_improved_lower_bounds, - h_improved_upper_bounds, + h_improved_bounds, h_integer_indices, n_of_implied_singletons); } @@ -462,14 +498,16 @@ void compute_probing_cache(bound_presolve_t& bound_presolve, const size_t max_threads = 10; omp_set_num_threads(max_threads); + load_balanced_problem_t lb_problem(problem); + // Create a vector of multi_probe_t objects - std::vector> multi_probe_presolve_pool; + std::vector>> multi_probe_presolve_pool; // Initialize multi_probe_presolve_pool for (size_t i = 0; i < max_threads; i++) { - multi_probe_presolve_pool.emplace_back(bound_presolve.context); - multi_probe_presolve_pool[i].resize(problem); - multi_probe_presolve_pool[i].compute_stats = false; + multi_probe_presolve_pool.emplace_back( + std::make_unique>(lb_problem, bound_presolve.context)); + multi_probe_presolve_pool[i]->compute_stats = false; } // Atomic variables for tracking progress @@ -486,12 +524,12 @@ void compute_probing_cache(bound_presolve_t& bound_presolve, int thread_idx = omp_get_thread_num(); CUOPT_LOG_TRACE("Computing probing cache for var %d on thread %d", var_idx, thread_idx); - auto& multi_probe_presolve = multi_probe_presolve_pool[thread_idx]; + auto multi_probe_presolve = multi_probe_presolve_pool[thread_idx].get(); compute_cache_for_var(var_idx, bound_presolve, problem, - multi_probe_presolve, + *multi_probe_presolve, h_var_lower_bounds, h_var_upper_bounds, h_integer_indices, diff --git a/cpp/tests/linear_programming/mip_unit_tests/multi_probe_test.cu b/cpp/tests/linear_programming/mip_unit_tests/multi_probe_test.cu index 014543fe72..5dad8cb76a 100644 --- a/cpp/tests/linear_programming/mip_unit_tests/multi_probe_test.cu +++ b/cpp/tests/linear_programming/mip_unit_tests/multi_probe_test.cu @@ -22,7 +22,10 @@ #include #include #include +#include +#include #include +#include #include #include #include @@ -98,6 +101,32 @@ convert_probe_tuple(std::tuple, std::vector, std::vecto return std::make_pair(std::move(probe_first), std::move(probe_second)); } +std::tuple, std::vector, std::vector, std::vector> +bounds_probe_results(detail::load_balanced_bounds_presolve_t& bnd_prb_0, + detail::load_balanced_bounds_presolve_t& bnd_prb_1, + detail::problem_t& problem, + const std::pair>, + std::vector>>& probe) +{ + auto& probe_first = std::get<0>(probe); + auto& probe_second = std::get<1>(probe); + rmm::device_uvector b_lb_0(problem.n_variables, problem.handle_ptr->get_stream()); + rmm::device_uvector b_ub_0(problem.n_variables, problem.handle_ptr->get_stream()); + rmm::device_uvector b_lb_1(problem.n_variables, problem.handle_ptr->get_stream()); + rmm::device_uvector b_ub_1(problem.n_variables, problem.handle_ptr->get_stream()); + bnd_prb_0.solve(probe_first); + bnd_prb_1.solve(probe_second); + bnd_prb_0.set_updated_bounds(b_lb_0, b_ub_0); + bnd_prb_1.set_updated_bounds(b_lb_1, b_ub_1); + + auto h_lb_0 = host_copy(b_lb_0); + auto h_ub_0 = host_copy(b_ub_0); + auto h_lb_1 = host_copy(b_lb_1); + auto h_ub_1 = host_copy(b_ub_1); + return std::make_tuple( + std::move(h_lb_0), std::move(h_ub_0), std::move(h_lb_1), std::move(h_ub_1)); +} + std::tuple, std::vector, std::vector, std::vector> bounds_probe_results(detail::bound_presolve_t& bnd_prb_0, detail::bound_presolve_t& bnd_prb_1, @@ -126,11 +155,11 @@ bounds_probe_results(detail::bound_presolve_t& bnd_prb_0, std::tuple, std::vector, std::vector, std::vector> multi_probe_results( - detail::multi_probe_t& prb, + detail::lb_multi_probe_t& prb, detail::problem_t& problem, const std::tuple, std::vector, std::vector>& probe_tuple) { - prb.solve(problem, probe_tuple); + prb.solve(probe_tuple); rmm::device_uvector m_lb_0(problem.n_variables, problem.handle_ptr->get_stream()); rmm::device_uvector m_ub_0(problem.n_variables, problem.handle_ptr->get_stream()); rmm::device_uvector m_lb_1(problem.n_variables, problem.handle_ptr->get_stream()); @@ -169,9 +198,11 @@ void test_multi_probe(std::string path) problem.reverse_constraints, true); detail::mip_solver_t solver(problem, default_settings, scaling, cuopt::timer_t(0)); + detail::load_balanced_problem_t lb_problem(problem); detail::bound_presolve_t bnd_prb_0(solver.context); detail::bound_presolve_t bnd_prb_1(solver.context); - detail::multi_probe_t multi_probe_prs(solver.context); + + detail::lb_multi_probe_t multi_probe_prs(lb_problem, solver.context); auto probe_tuple = select_k_random(problem, 100); auto bounds_probe_vals = convert_probe_tuple(probe_tuple); @@ -181,28 +212,11 @@ void test_multi_probe(std::string path) auto [m_lb_0, m_ub_0, m_lb_1, m_ub_1] = multi_probe_results(multi_probe_prs, problem, probe_tuple); - auto bnd_min_act_0 = host_copy(bnd_prb_0.upd.min_activity); - auto bnd_max_act_0 = host_copy(bnd_prb_0.upd.max_activity); - auto bnd_min_act_1 = host_copy(bnd_prb_1.upd.min_activity); - auto bnd_max_act_1 = host_copy(bnd_prb_1.upd.max_activity); - - auto mlp_min_act_0 = host_copy(multi_probe_prs.upd_0.min_activity); - auto mlp_max_act_0 = host_copy(multi_probe_prs.upd_0.max_activity); - auto mlp_min_act_1 = host_copy(multi_probe_prs.upd_1.min_activity); - auto mlp_max_act_1 = host_copy(multi_probe_prs.upd_1.max_activity); - - for (int i = 0; i < (int)bnd_min_act_0.size(); ++i) { - EXPECT_DOUBLE_EQ(bnd_min_act_0[i], mlp_min_act_0[i]); - EXPECT_DOUBLE_EQ(bnd_max_act_0[i], mlp_max_act_0[i]); - EXPECT_DOUBLE_EQ(bnd_min_act_1[i], mlp_min_act_1[i]); - EXPECT_DOUBLE_EQ(bnd_max_act_1[i], mlp_max_act_1[i]); - } - for (int i = 0; i < (int)bnd_lb_0.size(); ++i) { - EXPECT_DOUBLE_EQ(bnd_lb_0[i], m_lb_0[i]); - EXPECT_DOUBLE_EQ(bnd_ub_0[i], m_ub_0[i]); - EXPECT_DOUBLE_EQ(bnd_lb_1[i], m_lb_1[i]); - EXPECT_DOUBLE_EQ(bnd_ub_1[i], m_ub_1[i]); + EXPECT_DOUBLE_EQ(bnd_lb_0[i], m_lb_0[i]) << "index 0 " << i << std::endl; + EXPECT_DOUBLE_EQ(bnd_ub_0[i], m_ub_0[i]) << "index 0 " << i << std::endl; + EXPECT_DOUBLE_EQ(bnd_lb_1[i], m_lb_1[i]) << "index 1 " << i << std::endl; + EXPECT_DOUBLE_EQ(bnd_ub_1[i], m_ub_1[i]) << "index 1 " << i << std::endl; } }