From aa4d55c5c9ef7b8ef78d88fea36e381c4efb6f58 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Mon, 16 Mar 2026 16:59:00 -0700 Subject: [PATCH] Use a single dual simplex pivot to initialize pseudocost --- cpp/src/branch_and_bound/branch_and_bound.cpp | 5 + cpp/src/branch_and_bound/pseudo_costs.cpp | 286 +++++++++++++++- cpp/src/branch_and_bound/pseudo_costs.hpp | 7 +- cpp/src/dual_simplex/phase2.cpp | 306 ++++++++++-------- cpp/src/dual_simplex/phase2.hpp | 23 ++ cpp/src/math_optimization/solver_settings.cu | 2 +- 6 files changed, 481 insertions(+), 148 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index ba2b63983e..dac6c82a24 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2454,10 +2454,15 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut exploration_stats_.start_time, var_types_, root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, fractional, root_objective_, root_vstatus_, edge_norms_, + basic_list, + nonbasic_list, + basis_update, pc_); } diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index ee7e2f7803..d86c06740e 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -22,6 +22,256 @@ namespace cuopt::linear_programming::dual_simplex { namespace { + +template +struct objective_estimate_t { + f_t down_obj; + f_t up_obj; +}; + +template +f_t compute_step_length(const simplex_solver_settings_t& settings, + const std::vector& vstatus, + const std::vector& z, + const std::vector& delta_z, + const std::vector& delta_z_indices) +{ + f_t step_length = inf; + f_t pivot_tol = settings.pivot_tol; + const i_t nz = delta_z_indices.size(); + for (i_t h = 0; h < nz; h++) { + const i_t j = delta_z_indices[h]; + if (vstatus[j] == variable_status_t::NONBASIC_LOWER && delta_z[j] < -pivot_tol) { + const f_t ratio = -z[j] / delta_z[j]; + if (ratio < step_length) { + step_length = ratio; + } + } else if (vstatus[j] == variable_status_t::NONBASIC_UPPER && delta_z[j] > pivot_tol) { + const f_t ratio = -z[j] / delta_z[j]; + if (ratio < step_length) { + step_length = ratio; + } + } + } + return step_length; +} + + +template +objective_estimate_t single_pivot_objective_estimate( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const csc_matrix_t& A_transpose, + const std::vector& vstatus, + i_t variable_j, + i_t basic_j, + f_t root_obj, + const std::vector& x, + const std::vector& y, + const std::vector& z, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& nonbasic_mark, + basis_update_mpf_t& basis_factors, + std::vector& workspace, + std::vector& delta_z, + f_t& work_estimate) +{ + // Compute the objective estimate for the down and up branches of variable j + + // Down branch + i_t direction = -1; + sparse_vector_t e_k(lp.num_rows, 0); + e_k.i.push_back(basic_j); + e_k.x.push_back(-f_t(direction)); + + sparse_vector_t delta_y(lp.num_rows, 0); + basis_factors.b_transpose_solve(e_k, delta_y); + + // Compute delta_z_N = -N^T * delta_y + i_t delta_y_nz0 = 0; + const i_t nz_delta_y = delta_y.i.size(); + for (i_t k = 0; k < nz_delta_y; k++) { + if (std::abs(delta_y.x[k]) > 1e-12) { delta_y_nz0++; } + } + work_estimate += nz_delta_y; + const f_t delta_y_nz_percentage = delta_y_nz0 / static_cast(lp.num_rows) * 100.0; + const bool use_transpose = delta_y_nz_percentage <= 30.0; + std::vector delta_z_indices; + // delta_z starts out all zero + if (use_transpose) { + compute_delta_z(A_transpose, + delta_y, + variable_j, + -1, + nonbasic_mark, + workspace, + delta_z_indices, + delta_z, + work_estimate); + } else { + std::vector delta_y_dense(lp.num_rows, 0); + delta_y.to_dense(delta_y_dense); + compute_reduced_cost_update(lp, + basic_list, + nonbasic_list, + delta_y_dense, + variable_j, + -1, + workspace, + delta_z_indices, + delta_z, + work_estimate); + } + + + + // Verify dual feasibility +#ifdef CHECK_DUAL_FEASIBILITY + { + std::vector dual_residual = z; + for (i_t j = 0; j < lp.num_cols; j++) { + dual_residual[j] -= lp.objective[j]; + } + matrix_transpose_vector_multiply(lp.A, 1.0, y, 1.0, dual_residual); + f_t dual_residual_norm = vector_norm_inf(dual_residual); + settings.log.printf("Dual residual norm: %e\n", dual_residual_norm); + } +#endif + + // Compute the step-length + f_t step_length = compute_step_length(settings, vstatus, z, delta_z, delta_z_indices); + + // Handle the leaving variable case + + f_t delta_obj = step_length * (x[variable_j] - std::floor(x[variable_j])); +#ifdef CHECK_DELTA_OBJ + f_t delta_obj_check = 0.0; + for (i_t k = 0; k < delta_y.i.size(); k++) { + delta_obj_check += lp.rhs[delta_y.i[k]] * delta_y.x[k]; + } + for (i_t h = 0; h < delta_z_indices.size(); h++) { + const i_t j = delta_z_indices[h]; + if (vstatus[j] == variable_status_t::NONBASIC_LOWER) { + delta_obj_check += lp.lower[j] * delta_z[j]; + } else if (vstatus[j] == variable_status_t::NONBASIC_UPPER) { + delta_obj_check += lp.upper[j] * delta_z[j]; + } + } + delta_obj_check += std::floor(x[variable_j]) * delta_z[variable_j]; + delta_obj_check *= step_length; + if (std::abs(delta_obj_check - delta_obj) > 1e-6) { + settings.log.printf("Delta obj check %e. Delta obj %e. Step length %e.\n", delta_obj_check, delta_obj, step_length); + } +#endif + + f_t down_obj = root_obj + delta_obj; + settings.log.printf("Down branch %d. Step length: %e. Objective estimate: %e. Delta obj: %e. Root obj: %e.\n", variable_j, step_length, down_obj, delta_obj, root_obj); + + // Up branch + direction = 1; + // Negate delta_z + for (i_t j : delta_z_indices) + { + delta_z[j] *= -1.0; + } + + // Compute the step-length + step_length = compute_step_length(settings, vstatus, z, delta_z, delta_z_indices); + + delta_obj = step_length * (std::ceil(x[variable_j]) - x[variable_j]); + f_t up_obj = root_obj + delta_obj; + settings.log.printf("Up branch %d. Step length: %e. Objective estimate: %e. Delta obj: %e. Root obj: %e.\n", variable_j, step_length, up_obj, delta_obj, root_obj); + + delta_z_indices.push_back(variable_j); + + // Clear delta_z + for (i_t j : delta_z_indices) + { + delta_z[j] = 0.0; + workspace[j] = 0; + } + +#ifdef CHECK_DELTA_Z + for (i_t j = 0; j < lp.num_cols; j++) { + if (delta_z[j] != 0.0) { + settings.log.printf("Delta z %d: %e\n", j, delta_z[j]); + } + } + for (i_t j = 0; j < lp.num_cols; j++) { + if (workspace[j] != 0) { + settings.log.printf("Workspace %d: %d\n", j, workspace[j]); + } + } +#endif + + + return objective_estimate_t{.down_obj = down_obj, .up_obj = up_obj}; +} + + +template +void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& vstatus, + const std::vector& x, + const std::vector& y, + const std::vector& z, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& fractional, + f_t root_obj, + basis_update_mpf_t& basis_factors, + pseudo_costs_t& pc) +{ + i_t m = lp.num_rows; + i_t n = lp.num_cols; + + csc_matrix_t A_transpose(m, n, 0); + lp.A.transpose(A_transpose); + + std::vector delta_z(n, 0); + std::vector workspace(n, 0); + + f_t work_estimate = 0; + + std::vector basic_map(n, -1); + for (i_t i = 0; i < m; i++) { + basic_map[basic_list[i]] = i; + } + + std::vector nonbasic_mark(n, -1); + for (i_t i = 0; i < n - m; i++) { + nonbasic_mark[nonbasic_list[i]] = i; + } + + for (i_t k = 0; k < fractional.size(); k++) { + const i_t j = fractional[k]; + objective_estimate_t estimate = single_pivot_objective_estimate(lp, + settings, + A_transpose, + vstatus, + j, + basic_map[j], + root_obj, + x, + y, + z, + basic_list, + nonbasic_list, + nonbasic_mark, + basis_factors, + workspace, + delta_z, + work_estimate); + pc.strong_branch_down[k] = estimate.down_obj; + pc.strong_branch_up[k] = estimate.up_obj; + } +} + + + + template void strong_branch_helper(i_t start, i_t end, @@ -302,11 +552,16 @@ void strong_branching(const user_problem_t& original_problem, const simplex_solver_settings_t& settings, f_t start_time, const std::vector& var_types, - const std::vector root_soln, + const std::vector& root_x, + const std::vector& root_y, + const std::vector& root_z, const std::vector& fractional, f_t root_obj, const std::vector& root_vstatus, const std::vector& edge_norms, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& basis_factors, pseudo_costs_t& pc) { pc.resize(original_lp.num_cols); @@ -317,7 +572,7 @@ void strong_branching(const user_problem_t& original_problem, const f_t elapsed_time = toc(start_time); if (elapsed_time > settings.time_limit) { return; } - if (settings.mip_batch_pdlp_strong_branching) { + if (settings.mip_batch_pdlp_strong_branching == 1) { settings.log.printf("Batch PDLP strong branching enabled\n"); f_t start_batch = tic(); @@ -328,7 +583,7 @@ void strong_branching(const user_problem_t& original_problem, // Convert the root_soln to the original problem space std::vector original_root_soln_x; - uncrush_primal_solution(original_problem, original_lp, root_soln, original_root_soln_x); + uncrush_primal_solution(original_problem, original_lp, root_x, original_root_soln_x); std::vector fraction_values; @@ -394,6 +649,20 @@ void strong_branching(const user_problem_t& original_problem, pc.strong_branch_down[k] = obj_down - root_obj; pc.strong_branch_up[k] = obj_up - root_obj; } + } else if (settings.mip_batch_pdlp_strong_branching == 2) { + initialize_pseudo_costs_with_estimate(original_lp, + settings, + root_vstatus, + root_x, + root_y, + root_z, + basic_list, + nonbasic_list, + fractional, + root_obj, + basis_factors, + pc); + } else { settings.log.printf("Strong branching using %d threads and %ld fractional variables\n", settings.num_threads, @@ -429,7 +698,7 @@ void strong_branching(const user_problem_t& original_problem, var_types, fractional, root_obj, - root_soln, + root_x, root_vstatus, edge_norms, pc); @@ -438,7 +707,7 @@ void strong_branching(const user_problem_t& original_problem, settings.log.printf("Strong branching completed in %.2fs\n", toc(strong_branching_start_time)); } - pc.update_pseudo_costs_from_strong_branching(fractional, root_soln); + pc.update_pseudo_costs_from_strong_branching(fractional, root_x); } template @@ -781,11 +1050,16 @@ template void strong_branching(const user_problem_t& o const simplex_solver_settings_t& settings, double start_time, const std::vector& var_types, - const std::vector root_soln, + const std::vector& root_x, + const std::vector& root_y, + const std::vector& root_z, const std::vector& fractional, double root_obj, const std::vector& root_vstatus, const std::vector& edge_norms, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& basis_factors, pseudo_costs_t& pc); #endif diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index 6b6c6917b6..e34dc5b6a4 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -522,11 +522,16 @@ void strong_branching(const user_problem_t& original_problem, const simplex_solver_settings_t& settings, f_t start_time, const std::vector& var_types, - const std::vector root_soln, + const std::vector& root_x, + const std::vector& root_y, + const std::vector& root_z, const std::vector& fractional, f_t root_obj, const std::vector& root_vstatus, const std::vector& edge_norms, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& basis_factors, pseudo_costs_t& pc); } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 43429ba2de..eca9dbd251 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -82,6 +82,130 @@ class nvtx_range_guard { bool active_; }; + +template +void compute_reduced_cost_update(const lp_problem_t& lp, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& delta_y, + i_t leaving_index, + i_t direction, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + f_t& work_estimate) +{ + const i_t m = lp.num_rows; + const i_t n = lp.num_cols; + + size_t nnzs_processed = 0; + // delta_zB = sigma*ei + for (i_t k = 0; k < m; k++) { + const i_t j = basic_list[k]; + delta_z[j] = 0; + } + work_estimate += 2 * m; + delta_z[leaving_index] = direction; + // delta_zN = -N'*delta_y + const i_t num_nonbasic = n - m; + for (i_t k = 0; k < num_nonbasic; k++) { + const i_t j = nonbasic_list[k]; + // z_j <- -A(:, j)'*delta_y + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + f_t dot = 0.0; + for (i_t p = col_start; p < col_end; ++p) { + dot += lp.A.x[p] * delta_y[lp.A.i[p]]; + } + nnzs_processed += col_end - col_start; + + delta_z[j] = -dot; + if (dot != 0.0) { + delta_z_indices.push_back(j); // Note delta_z_indices has n elements reserved + delta_z_mark[j] = 1; + } + } + work_estimate += 3 * num_nonbasic; + work_estimate += 3 * nnzs_processed; + work_estimate += 2 * delta_z_indices.size(); +} + +template +void compute_delta_z(const csc_matrix_t& A_transpose, + const sparse_vector_t& delta_y, + i_t leaving_index, + i_t direction, + const std::vector& nonbasic_mark, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + f_t& work_estimate) +{ + // delta_zN = - N'*delta_y + const i_t nz_delta_y = delta_y.i.size(); + size_t nnz_processed = 0; + size_t nonbasic_marked = 0; + for (i_t k = 0; k < nz_delta_y; k++) { + const i_t i = delta_y.i[k]; + const f_t delta_y_i = delta_y.x[k]; + if (std::abs(delta_y_i) < 1e-12) { continue; } + const i_t row_start = A_transpose.col_start[i]; + const i_t row_end = A_transpose.col_start[i + 1]; + nnz_processed += row_end - row_start; + for (i_t p = row_start; p < row_end; ++p) { + const i_t j = A_transpose.i[p]; + if (nonbasic_mark[j] >= 0) { + delta_z[j] -= delta_y_i * A_transpose.x[p]; + nonbasic_marked++; + if (!delta_z_mark[j]) { + delta_z_mark[j] = 1; + delta_z_indices.push_back(j); + } + } + } + } + work_estimate += 4 * nz_delta_y; + work_estimate += 2 * nnz_processed; + work_estimate += 3 * nonbasic_marked; + work_estimate += 2 * delta_z_indices.size(); + + // delta_zB = sigma*ei + delta_z[leaving_index] = direction; + +#ifdef CHECK_CHANGE_IN_REDUCED_COST + const i_t m = A_transpose.n; + const i_t n = A_transpose.m; + std::vector delta_y_dense(m); + delta_y.to_dense(delta_y_dense); + std::vector delta_z_check(n); + std::vector delta_z_mark_check(n, 0); + std::vector delta_z_indices_check; + phase2::compute_reduced_cost_update(lp, + basic_list, + nonbasic_list, + delta_y_dense, + leaving_index, + direction, + delta_z_mark_check, + delta_z_indices_check, + delta_z_check, + work_estimate); + f_t error_check = 0.0; + for (i_t k = 0; k < n; ++k) { + const f_t diff = std::abs(delta_z[k] - delta_z_check[k]); + if (diff > 1e-6) { + printf("delta_z error %d transpose %e no transpose %e diff %e\n", + k, + delta_z[k], + delta_z_check[k], + diff); + } + error_check = std::max(error_check, diff); + } + if (error_check > 1e-6) { printf("delta_z error %e\n", error_check); } +#endif +} + namespace phase2 { // Computes vectors farkas_y, farkas_zl, farkas_zu that satisfy @@ -322,128 +446,7 @@ void initial_perturbation(const lp_problem_t& lp, n); } -template -void compute_reduced_cost_update(const lp_problem_t& lp, - const std::vector& basic_list, - const std::vector& nonbasic_list, - const std::vector& delta_y, - i_t leaving_index, - i_t direction, - std::vector& delta_z_mark, - std::vector& delta_z_indices, - std::vector& delta_z, - f_t& work_estimate) -{ - const i_t m = lp.num_rows; - const i_t n = lp.num_cols; - - size_t nnzs_processed = 0; - // delta_zB = sigma*ei - for (i_t k = 0; k < m; k++) { - const i_t j = basic_list[k]; - delta_z[j] = 0; - } - work_estimate += 2 * m; - delta_z[leaving_index] = direction; - // delta_zN = -N'*delta_y - const i_t num_nonbasic = n - m; - for (i_t k = 0; k < num_nonbasic; k++) { - const i_t j = nonbasic_list[k]; - // z_j <- -A(:, j)'*delta_y - const i_t col_start = lp.A.col_start[j]; - const i_t col_end = lp.A.col_start[j + 1]; - f_t dot = 0.0; - for (i_t p = col_start; p < col_end; ++p) { - dot += lp.A.x[p] * delta_y[lp.A.i[p]]; - } - nnzs_processed += col_end - col_start; - - delta_z[j] = -dot; - if (dot != 0.0) { - delta_z_indices.push_back(j); // Note delta_z_indices has n elements reserved - delta_z_mark[j] = 1; - } - } - work_estimate += 3 * num_nonbasic; - work_estimate += 3 * nnzs_processed; - work_estimate += 2 * delta_z_indices.size(); -} - -template -void compute_delta_z(const csc_matrix_t& A_transpose, - const sparse_vector_t& delta_y, - i_t leaving_index, - i_t direction, - std::vector& nonbasic_mark, - std::vector& delta_z_mark, - std::vector& delta_z_indices, - std::vector& delta_z, - f_t& work_estimate) -{ - // delta_zN = - N'*delta_y - const i_t nz_delta_y = delta_y.i.size(); - size_t nnz_processed = 0; - size_t nonbasic_marked = 0; - for (i_t k = 0; k < nz_delta_y; k++) { - const i_t i = delta_y.i[k]; - const f_t delta_y_i = delta_y.x[k]; - if (std::abs(delta_y_i) < 1e-12) { continue; } - const i_t row_start = A_transpose.col_start[i]; - const i_t row_end = A_transpose.col_start[i + 1]; - nnz_processed += row_end - row_start; - for (i_t p = row_start; p < row_end; ++p) { - const i_t j = A_transpose.i[p]; - if (nonbasic_mark[j] >= 0) { - delta_z[j] -= delta_y_i * A_transpose.x[p]; - nonbasic_marked++; - if (!delta_z_mark[j]) { - delta_z_mark[j] = 1; - delta_z_indices.push_back(j); - } - } - } - } - work_estimate += 4 * nz_delta_y; - work_estimate += 2 * nnz_processed; - work_estimate += 3 * nonbasic_marked; - work_estimate += 2 * delta_z_indices.size(); - // delta_zB = sigma*ei - delta_z[leaving_index] = direction; - -#ifdef CHECK_CHANGE_IN_REDUCED_COST - const i_t m = A_transpose.n; - const i_t n = A_transpose.m; - std::vector delta_y_dense(m); - delta_y.to_dense(delta_y_dense); - std::vector delta_z_check(n); - std::vector delta_z_mark_check(n, 0); - std::vector delta_z_indices_check; - phase2::compute_reduced_cost_update(lp, - basic_list, - nonbasic_list, - delta_y_dense, - leaving_index, - direction, - delta_z_mark_check, - delta_z_indices_check, - delta_z_check, - work_estimate); - f_t error_check = 0.0; - for (i_t k = 0; k < n; ++k) { - const f_t diff = std::abs(delta_z[k] - delta_z_check[k]); - if (diff > 1e-6) { - printf("delta_z error %d transpose %e no transpose %e diff %e\n", - k, - delta_z[k], - delta_z_check[k], - diff); - } - error_check = std::max(error_check, diff); - } - if (error_check > 1e-6) { printf("delta_z error %e\n", error_check); } -#endif -} template void compute_reduced_costs(const std::vector& objective, @@ -2932,30 +2935,30 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, PHASE2_NVTX_RANGE("DualSimplex::delta_z"); if (use_transpose) { sparse_delta_z++; - phase2::compute_delta_z(A_transpose, - delta_y_sparse, - leaving_index, - direction, - nonbasic_mark, - delta_z_mark, - delta_z_indices, - delta_z, - phase2_work_estimate); + compute_delta_z(A_transpose, + delta_y_sparse, + leaving_index, + direction, + nonbasic_mark, + delta_z_mark, + delta_z_indices, + delta_z, + phase2_work_estimate); } else { dense_delta_z++; // delta_zB = sigma*ei delta_y_sparse.to_dense(delta_y); phase2_work_estimate += delta_y.size(); - phase2::compute_reduced_cost_update(lp, - basic_list, - nonbasic_list, - delta_y, - leaving_index, - direction, - delta_z_mark, - delta_z_indices, - delta_z, - phase2_work_estimate); + compute_reduced_cost_update(lp, + basic_list, + nonbasic_list, + delta_y, + leaving_index, + direction, + delta_z_mark, + delta_z_indices, + delta_z, + phase2_work_estimate); } } timers.delta_z_time += timers.stop_timer(); @@ -3634,6 +3637,29 @@ template dual::status_t dual_phase2_with_advanced_basis( int& iter, std::vector& steepest_edge_norms, work_limit_context_t* work_unit_context); + +template +void compute_reduced_cost_update(const lp_problem_t& lp, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& delta_y, + int leaving_index, + int direction, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + double& work_estimate); + +template +void compute_delta_z(const csc_matrix_t& A_transpose, + const sparse_vector_t& delta_y, + int leaving_index, + int direction, + const std::vector& nonbasic_mark, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + double& work_estimate); #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp index 7f01eb3cf7..5db797449c 100644 --- a/cpp/src/dual_simplex/phase2.hpp +++ b/cpp/src/dual_simplex/phase2.hpp @@ -81,4 +81,27 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, std::vector& delta_y_steepest_edge, work_limit_context_t* work_unit_context = nullptr); +template +void compute_reduced_cost_update(const lp_problem_t& lp, + const std::vector& basic_list, + const std::vector& nonbasic_list, + const std::vector& delta_y, + i_t leaving_index, + i_t direction, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + f_t& work_estimate); + +template +void compute_delta_z(const csc_matrix_t& A_transpose, + const sparse_vector_t& delta_y, + i_t leaving_index, + i_t direction, + const std::vector& nonbasic_mark, + std::vector& delta_z_mark, + std::vector& delta_z_indices, + std::vector& delta_z, + f_t& work_estimate); + } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index a60d508fac..13a8f74145 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -99,7 +99,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_NUM_GPUS, &pdlp_settings.num_gpus, 1, 2, 1}, {CUOPT_NUM_GPUS, &mip_settings.num_gpus, 1, 2, 1}, - {CUOPT_MIP_BATCH_PDLP_STRONG_BRANCHING, &mip_settings.mip_batch_pdlp_strong_branching, 0, 1, 0}, + {CUOPT_MIP_BATCH_PDLP_STRONG_BRANCHING, &mip_settings.mip_batch_pdlp_strong_branching, 0, 2, 0}, {CUOPT_PRESOLVE, reinterpret_cast(&pdlp_settings.presolver), CUOPT_PRESOLVE_DEFAULT, CUOPT_PRESOLVE_PSLP, CUOPT_PRESOLVE_DEFAULT}, {CUOPT_PRESOLVE, reinterpret_cast(&mip_settings.presolver), CUOPT_PRESOLVE_DEFAULT, CUOPT_PRESOLVE_PSLP, CUOPT_PRESOLVE_DEFAULT}, {CUOPT_MIP_DETERMINISM_MODE, &mip_settings.determinism_mode, CUOPT_MODE_OPPORTUNISTIC, CUOPT_MODE_DETERMINISTIC, CUOPT_MODE_OPPORTUNISTIC},