diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index c46f09258..6ce9a4f4d 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1886,6 +1886,7 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( root_crossover_settings, original_lp_.lower, original_lp_.upper, + exploration_stats_.start_time, basic_list, nonbasic_list, crossover_vstatus_); @@ -2285,6 +2286,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut mutex_original_lp_.lock(); remove_cuts(original_lp_, settings_, + exploration_stats_.start_time, Arow_, new_slacks_, original_rows, diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 84ee986c4..a13d5cedc 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -106,6 +106,8 @@ class branch_and_bound_t { void set_concurrent_lp_root_solve(bool enable) { enable_concurrent_lp_root_solve_ = enable; } + bool stop_for_time_limit(mip_solution_t& solution); + // Repair a low-quality solution from the heuristics. bool repair_solution(const std::vector& leaf_edge_norms, const std::vector& potential_solution, diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index 095f46238..ee7e2f780 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -314,6 +314,9 @@ void strong_branching(const user_problem_t& original_problem, pc.strong_branch_up.assign(fractional.size(), 0); pc.num_strong_branches_completed = 0; + const f_t elapsed_time = toc(start_time); + if (elapsed_time > settings.time_limit) { return; } + if (settings.mip_batch_pdlp_strong_branching) { settings.log.printf("Batch PDLP strong branching enabled\n"); @@ -334,10 +337,16 @@ void strong_branching(const user_problem_t& original_problem, fraction_values.push_back(original_root_soln_x[j]); } - const auto mps_model = simplex_problem_to_mps_data_model(original_problem); + const auto mps_model = simplex_problem_to_mps_data_model(original_problem); + const f_t batch_elapsed_time = toc(start_time); + const f_t batch_remaining_time = + std::max(static_cast(0.0), settings.time_limit - batch_elapsed_time); + if (batch_remaining_time <= 0.0) { return; } + pdlp_solver_settings_t pdlp_settings; + pdlp_settings.time_limit = batch_remaining_time; const raft::handle_t batch_pdlp_handle; const auto solutions = - batch_pdlp_solve(&batch_pdlp_handle, mps_model, fractional, fraction_values); + batch_pdlp_solve(&batch_pdlp_handle, mps_model, fractional, fraction_values, pdlp_settings); f_t batch_pdlp_strong_branching_time = toc(start_batch); // Find max iteration on how many are done accross the batch diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index cd3dbf000..ad1d34f7a 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -2501,20 +2501,21 @@ i_t add_cuts(const simplex_solver_settings_t& settings, } template -void remove_cuts(lp_problem_t& lp, - const simplex_solver_settings_t& settings, - csr_matrix_t& Arow, - std::vector& new_slacks, - i_t original_rows, - std::vector& var_types, - std::vector& vstatus, - std::vector& edge_norms, - std::vector& x, - std::vector& y, - std::vector& z, - std::vector& basic_list, - std::vector& nonbasic_list, - basis_update_mpf_t& basis_update) +i_t remove_cuts(lp_problem_t& lp, + const simplex_solver_settings_t& settings, + f_t start_time, + csr_matrix_t& Arow, + std::vector& new_slacks, + i_t original_rows, + std::vector& var_types, + std::vector& vstatus, + std::vector& edge_norms, + std::vector& x, + std::vector& y, + std::vector& z, + std::vector& basic_list, + std::vector& nonbasic_list, + basis_update_mpf_t& basis_update) { std::vector cuts_to_remove; cuts_to_remove.reserve(lp.num_rows - original_rows); @@ -2644,9 +2645,13 @@ void remove_cuts(lp_problem_t& lp, lp.A.col_start[lp.A.n]); basis_update.resize(lp.num_rows); - basis_update.refactor_basis( - lp.A, settings, lp.lower, lp.upper, basic_list, nonbasic_list, vstatus); + i_t refactor_status = basis_update.refactor_basis( + lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + if (refactor_status == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; } + if (refactor_status == TIME_LIMIT_RETURN) { return TIME_LIMIT_RETURN; } } + + return 0; } template @@ -2793,20 +2798,21 @@ template int add_cuts(const simplex_solver_settings_t& settings, std::vector& vstatus, std::vector& edge_norms); -template void remove_cuts(lp_problem_t& lp, - const simplex_solver_settings_t& settings, - csr_matrix_t& Arow, - std::vector& new_slacks, - int original_rows, - std::vector& var_types, - std::vector& vstatus, - std::vector& edge_norms, - std::vector& x, - std::vector& y, - std::vector& z, - std::vector& basic_list, - std::vector& nonbasic_list, - basis_update_mpf_t& basis_update); +template int remove_cuts(lp_problem_t& lp, + const simplex_solver_settings_t& settings, + double start_time, + csr_matrix_t& Arow, + std::vector& new_slacks, + int original_rows, + std::vector& var_types, + std::vector& vstatus, + std::vector& edge_norms, + std::vector& x, + std::vector& y, + std::vector& z, + std::vector& basic_list, + std::vector& nonbasic_list, + basis_update_mpf_t& basis_update); template void read_saved_solution_for_cut_verification( const lp_problem_t& lp, diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index a4a36d75b..4f55e96e4 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -461,19 +461,20 @@ i_t add_cuts(const simplex_solver_settings_t& settings, std::vector& edge_norms); template -void remove_cuts(lp_problem_t& lp, - const simplex_solver_settings_t& settings, - csr_matrix_t& Arow, - std::vector& new_slacks, - i_t original_rows, - std::vector& var_types, - std::vector& vstatus, - std::vector& edge_norms, - std::vector& x, - std::vector& y, - std::vector& z, - std::vector& basic_list, - std::vector& nonbasic_list, - basis_update_mpf_t& basis_update); +i_t remove_cuts(lp_problem_t& lp, + const simplex_solver_settings_t& settings, + f_t start_time, + csr_matrix_t& Arow, + std::vector& new_slacks, + i_t original_rows, + std::vector& var_types, + std::vector& vstatus, + std::vector& edge_norms, + std::vector& x, + std::vector& y, + std::vector& z, + std::vector& basic_list, + std::vector& nonbasic_list, + basis_update_mpf_t& basis_update); } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/basis_solves.cpp b/cpp/src/dual_simplex/basis_solves.cpp index e7bc4ae66..c5fee4e10 100644 --- a/cpp/src/dual_simplex/basis_solves.cpp +++ b/cpp/src/dual_simplex/basis_solves.cpp @@ -160,6 +160,7 @@ template i_t factorize_basis(const csc_matrix_t& A, const simplex_solver_settings_t& settings, const std::vector& basic_list, + f_t start_time, csc_matrix_t& L, csc_matrix_t& U, std::vector& p, @@ -383,15 +384,16 @@ i_t factorize_basis(const csc_matrix_t& A, settings, settings.threshold_partial_pivoting_tol, identity, + start_time, S_col_perm, SL, SU, S_perm_inv, work_estimate); if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { - settings.log.printf("Concurrent halt\n"); return CONCURRENT_HALT_RETURN; } + if (Srank < 0) { return Srank; } if (Srank != Sdim) { // Get the rank deficient columns deficient.clear(); @@ -618,7 +620,14 @@ i_t factorize_basis(const csc_matrix_t& A, q.resize(m); work_estimate += m; f_t fact_start = tic(); - rank = right_looking_lu(A, settings, medium_tol, basic_list, q, L, U, pinv, work_estimate); + rank = + right_looking_lu(A, settings, medium_tol, basic_list, start_time, q, L, U, pinv, work_estimate); + if (rank < 0) { + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + return CONCURRENT_HALT_RETURN; + } + return rank; + } inverse_permutation(pinv, p); work_estimate += 3 * pinv.size(); @@ -638,7 +647,6 @@ i_t factorize_basis(const csc_matrix_t& A, work_estimate += 3 * (m - rank); } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { - settings.log.printf("Concurrent halt\n"); return CONCURRENT_HALT_RETURN; } if (verbose) { @@ -938,6 +946,7 @@ template void get_basis_from_vstatus(int m, template int factorize_basis(const csc_matrix_t& A, const simplex_solver_settings_t& settings, const std::vector& basis_list, + double start_time, csc_matrix_t& L, csc_matrix_t& U, std::vector& p, diff --git a/cpp/src/dual_simplex/basis_solves.hpp b/cpp/src/dual_simplex/basis_solves.hpp index fbeba565f..76a9c396b 100644 --- a/cpp/src/dual_simplex/basis_solves.hpp +++ b/cpp/src/dual_simplex/basis_solves.hpp @@ -30,6 +30,7 @@ template i_t factorize_basis(const csc_matrix_t& A, const simplex_solver_settings_t& settings, const std::vector& basis_list, + f_t start_time, csc_matrix_t& L, csc_matrix_t& U, std::vector& p, diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 6af315d2f..9c56ada50 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -2343,6 +2343,7 @@ int basis_update_mpf_t::refactor_basis( const simplex_solver_settings_t& settings, const std::vector& lower, const std::vector& upper, + f_t start_time, std::vector& basic_list, std::vector& nonbasic_list, std::vector& vstatus) @@ -2360,6 +2361,7 @@ int basis_update_mpf_t::refactor_basis( i_t status = factorize_basis(A, settings, basic_list, + start_time, L0_, U0_, row_permutation_, @@ -2369,6 +2371,7 @@ int basis_update_mpf_t::refactor_basis( slacks_needed, work_estimate_); if (status == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; } + if (status == TIME_LIMIT_RETURN) { return TIME_LIMIT_RETURN; } if (status == -1) { settings.log.debug("Initial factorization failed\n"); basis_repair(A, @@ -2403,6 +2406,7 @@ int basis_update_mpf_t::refactor_basis( status = factorize_basis(A, settings, basic_list, + start_time, L0_, U0_, row_permutation_, @@ -2412,6 +2416,7 @@ int basis_update_mpf_t::refactor_basis( slacks_needed, work_estimate_); if (status == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; } + if (status == TIME_LIMIT_RETURN) { return TIME_LIMIT_RETURN; } if (status == -1) { #ifdef CHECK_L_FACTOR if (L0_.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); } diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp index bb2b42d8a..f313e15c1 100644 --- a/cpp/src/dual_simplex/basis_updates.hpp +++ b/cpp/src/dual_simplex/basis_updates.hpp @@ -382,6 +382,7 @@ class basis_update_mpf_t { const simplex_solver_settings_t& settings, const std::vector& lower, const std::vector& upper, + f_t start_time, std::vector& basic_list, std::vector& nonbasic_list, std::vector& vstatus); diff --git a/cpp/src/dual_simplex/crossover.cpp b/cpp/src/dual_simplex/crossover.cpp index 15f013a5b..988c9c50a 100644 --- a/cpp/src/dual_simplex/crossover.cpp +++ b/cpp/src/dual_simplex/crossover.cpp @@ -25,7 +25,7 @@ namespace { crossover_status_t return_to_status(int status) { - if (status == -1) { + if (status == TIME_LIMIT_RETURN) { return crossover_status_t::TIME_LIMIT; } else if (status == CONCURRENT_HALT_RETURN) { return crossover_status_t::CONCURRENT_LIMIT; @@ -506,10 +506,22 @@ i_t dual_push(const lp_problem_t& lp, std::vector deficient; std::vector slacks_needed; f_t work_estimate = 0; - i_t rank = factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed, work_estimate); + i_t rank = factorize_basis(lp.A, + settings, + basic_list, + start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + work_estimate); if (rank == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; + } else if (rank < 0) { + return rank; } else if (rank != m) { settings.log.printf("Failed to factorize basis. rank %d m %d\n", rank, m); basis_repair(lp.A, @@ -523,13 +535,22 @@ i_t dual_push(const lp_problem_t& lp, superbasic_list, vstatus, work_estimate); - rank = factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed, work_estimate); + rank = factorize_basis(lp.A, + settings, + basic_list, + start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + work_estimate); if (rank == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; - } else if (rank == -1) { - settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m); - return -1; + } else if (rank < 0) { + return rank; } else { settings.log.printf("Basis repaired\n"); } @@ -562,7 +583,7 @@ i_t dual_push(const lp_problem_t& lp, } if (toc(start_time) > settings.time_limit) { settings.log.printf("Crossover time exceeded\n"); - return -1; + return TIME_LIMIT_RETURN; } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Concurrent halt\n"); @@ -863,10 +884,22 @@ i_t primal_push(const lp_problem_t& lp, std::vector deficient; std::vector slacks_needed; f_t work_estimate = 0; - i_t rank = factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed, work_estimate); + i_t rank = factorize_basis(lp.A, + settings, + basic_list, + start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + work_estimate); if (rank == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; + } else if (rank < 0) { + return rank; } else if (rank != m) { settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m); basis_repair(lp.A, @@ -883,13 +916,22 @@ i_t primal_push(const lp_problem_t& lp, // We need to be careful. As basis_repair may have changed the superbasic list find_primal_superbasic_variables( lp, settings, solution, solution, vstatus, nonbasic_list, superbasic_list); - rank = factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed, work_estimate); + rank = factorize_basis(lp.A, + settings, + basic_list, + start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + work_estimate); if (rank == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; - } else if (rank == -1) { - settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m); - return -1; + } else if (rank < 0) { + return rank; } else { settings.log.debug("Basis repaired\n"); } @@ -919,7 +961,7 @@ i_t primal_push(const lp_problem_t& lp, if (toc(start_time) > settings.time_limit) { settings.log.printf("Crossover time limit exceeded\n"); - return -1; + return TIME_LIMIT_RETURN; } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Concurrent halt\n"); @@ -1228,9 +1270,19 @@ crossover_status_t crossover(const lp_problem_t& lp, std::vector deficient; std::vector slacks_needed; - rank = factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed, work_estimate); - if (rank == CONCURRENT_HALT_RETURN) { return crossover_status_t::CONCURRENT_LIMIT; } + rank = factorize_basis(lp.A, + settings, + basic_list, + start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + work_estimate); + if (rank < 0) { return return_to_status(rank); } if (rank != m) { settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m); basis_repair(lp.A, @@ -1244,13 +1296,23 @@ crossover_status_t crossover(const lp_problem_t& lp, superbasic_list, vstatus, work_estimate); - rank = factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed, work_estimate); + rank = factorize_basis(lp.A, + settings, + basic_list, + start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + work_estimate); if (rank == CONCURRENT_HALT_RETURN) { return crossover_status_t::CONCURRENT_LIMIT; - } else if (rank == -1) { + } else if (rank < 0) { settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m); - return crossover_status_t::NUMERICAL_ISSUES; + return return_to_status(rank); } else { settings.log.debug("Basis repaired\n"); } @@ -1400,10 +1462,20 @@ crossover_status_t crossover(const lp_problem_t& lp, nonbasic_list.clear(); superbasic_list.clear(); get_basis_from_vstatus(m, vstatus, basic_list, nonbasic_list, superbasic_list); - rank = factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed, work_estimate); - if (rank == CONCURRENT_HALT_RETURN) { - return crossover_status_t::CONCURRENT_LIMIT; + rank = factorize_basis(lp.A, + settings, + basic_list, + start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + work_estimate); + if (rank < 0) { + return return_to_status(rank); } else if (rank != m) { settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m); basis_repair(lp.A, @@ -1417,13 +1489,21 @@ crossover_status_t crossover(const lp_problem_t& lp, superbasic_list, vstatus, work_estimate); - rank = factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed, work_estimate); - if (rank == CONCURRENT_HALT_RETURN) { - return crossover_status_t::CONCURRENT_LIMIT; - } else if (rank == -1) { + rank = factorize_basis(lp.A, + settings, + basic_list, + start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + work_estimate); + if (rank < 0) { settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m); - return crossover_status_t::NUMERICAL_ISSUES; + return return_to_status(rank); } else { settings.log.debug("Basis repaired\n"); } diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index d136c500f..426d9a753 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -2535,10 +2535,11 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, assert(superbasic_list.size() == 0); assert(nonbasic_list.size() == n - m); - if (ft.refactor_basis(lp.A, settings, lp.lower, lp.upper, basic_list, nonbasic_list, vstatus) > - 0) { - return dual::status_t::NUMERICAL; - } + i_t refactor_status = ft.refactor_basis( + lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + if (refactor_status == CONCURRENT_HALT_RETURN) { return dual::status_t::CONCURRENT_LIMIT; } + if (refactor_status == TIME_LIMIT_RETURN) { return dual::status_t::TIME_LIMIT; } + if (refactor_status > 0) { return dual::status_t::NUMERICAL; } if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; } } @@ -3406,15 +3407,24 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, PHASE2_NVTX_RANGE("DualSimplex::refactorization"); num_refactors++; bool should_recompute_x = false; - if (ft.refactor_basis( - lp.A, settings, lp.lower, lp.upper, basic_list, nonbasic_list, vstatus) > 0) { + i_t refactor_status = ft.refactor_basis( + lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + if (refactor_status == CONCURRENT_HALT_RETURN) { return dual::status_t::CONCURRENT_LIMIT; } + if (refactor_status == TIME_LIMIT_RETURN) { return dual::status_t::TIME_LIMIT; } + if (refactor_status > 0) { should_recompute_x = true; settings.log.printf("Failed to factorize basis. Iteration %d\n", iter); if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; } - i_t count = 0; - i_t deficient_size; - while ((deficient_size = ft.refactor_basis( - lp.A, settings, lp.lower, lp.upper, basic_list, nonbasic_list, vstatus)) > 0) { + i_t count = 0; + i_t deficient_size = 0; + while (true) { + deficient_size = ft.refactor_basis( + lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + if (deficient_size == CONCURRENT_HALT_RETURN) { + return dual::status_t::CONCURRENT_LIMIT; + } + if (deficient_size == TIME_LIMIT_RETURN) { return dual::status_t::TIME_LIMIT; } + if (deficient_size <= 0) { break; } settings.log.printf("Failed to repair basis. Iteration %d. %d deficient columns.\n", iter, static_cast(deficient_size)); @@ -3425,6 +3435,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, count++; if (count > 10) { return dual::status_t::NUMERICAL; } } + if (deficient_size < 0) { return dual::status_t::NUMERICAL; } settings.log.printf("Successfully repaired basis. Iteration %d\n", iter); } diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 8d4a533e9..b9ee41951 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -445,6 +445,7 @@ i_t find_dependent_rows(lp_problem_t& problem, i_t pivots = right_looking_lu_row_permutation_only(C, settings, 1e-13, tic(), q, pinv); if (pivots == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; } + if (pivots == TIME_LIMIT_RETURN) { return TIME_LIMIT_RETURN; } if (pivots < m) { settings.log.printf("Found %d dependent rows\n", m - pivots); const i_t num_dependent = m - pivots; @@ -1101,6 +1102,7 @@ i_t presolve(const lp_problem_t& original, f_t dependent_row_start = tic(); const i_t independent_rows = find_dependent_rows(problem, settings, dependent_rows, infeasible); if (independent_rows == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; } + if (independent_rows == TIME_LIMIT_RETURN) { return TIME_LIMIT_RETURN; } if (infeasible != kOk) { settings.log.printf("Found problem infeasible in presolve\n"); return -1; @@ -1161,7 +1163,9 @@ void crush_primal_solution(const user_problem_t& user_problem, const std::vector& new_slacks, std::vector& solution) { - solution.resize(problem.num_cols, 0.0); + // Re-crush can be called with a reused output vector; make sure all entries, + // including previously added slacks, are reset before writing new values. + solution.assign(problem.num_cols, 0.0); for (i_t j = 0; j < user_problem.num_cols; j++) { solution[j] = user_solution[j]; } @@ -1200,7 +1204,8 @@ void crush_primal_solution_with_slack(const user_problem_t& user_probl const std::vector& new_slacks, std::vector& solution) { - solution.resize(problem.num_cols, 0.0); + // Re-crush can be called with a reused output vector; clear stale entries first. + solution.assign(problem.num_cols, 0.0); for (i_t j = 0; j < user_problem.num_cols; j++) { solution[j] = user_solution[j]; } diff --git a/cpp/src/dual_simplex/primal.cpp b/cpp/src/dual_simplex/primal.cpp index b1b6be0b7..d4c6743dc 100644 --- a/cpp/src/dual_simplex/primal.cpp +++ b/cpp/src/dual_simplex/primal.cpp @@ -295,10 +295,25 @@ primal::status_t primal_phase2(i_t phase, std::vector q(m); std::vector deficient; std::vector slacks_needed; - i_t rank = factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed, work_estimate); + i_t rank = factorize_basis(lp.A, + settings, + basic_list, + start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + work_estimate); if (rank == CONCURRENT_HALT_RETURN) { return primal::status_t::CONCURRENT_LIMIT; + } else if (rank == TIME_LIMIT_RETURN) { + return primal::status_t::TIME_LIMIT; + } else if (rank < 0) { + return toc(start_time) > settings.time_limit ? primal::status_t::TIME_LIMIT + : primal::status_t::NUMERICAL; } else if (rank != m) { settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m); basis_repair(lp.A, @@ -312,13 +327,26 @@ primal::status_t primal_phase2(i_t phase, superbasic_list, vstatus, work_estimate); - rank = factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed, work_estimate); + rank = factorize_basis(lp.A, + settings, + basic_list, + start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + work_estimate); if (rank == CONCURRENT_HALT_RETURN) { return primal::status_t::CONCURRENT_LIMIT; - } else if (rank == -1) { + } else if (rank == TIME_LIMIT_RETURN) { + return primal::status_t::TIME_LIMIT; + } else if (rank < 0) { settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m); - return primal::status_t::NUMERICAL; + return toc(start_time) > settings.time_limit ? primal::status_t::TIME_LIMIT + : primal::status_t::NUMERICAL; } else { settings.log.debug("Basis repaired\n"); } diff --git a/cpp/src/dual_simplex/right_looking_lu.cpp b/cpp/src/dual_simplex/right_looking_lu.cpp index 8b62612d0..657ebc476 100644 --- a/cpp/src/dual_simplex/right_looking_lu.cpp +++ b/cpp/src/dual_simplex/right_looking_lu.cpp @@ -625,6 +625,7 @@ i_t right_looking_lu(const csc_matrix_t& A, const simplex_solver_settings_t& settings, f_t tol, const std::vector& column_list, + f_t start_time, std::vector& q, csc_matrix_t& L, csc_matrix_t& U, @@ -693,7 +694,10 @@ i_t right_looking_lu(const csc_matrix_t& A, i_t pivots = 0; for (i_t k = 0; k < n; ++k) { - if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return -1; } + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + return CONCURRENT_HALT_RETURN; + } + if (toc(start_time) > settings.time_limit) { return TIME_LIMIT_RETURN; } // Find pivot that satisfies // abs(pivot) >= abstol, // abs(pivot) >= threshold_tol * max abs[pivot column] @@ -1196,11 +1200,7 @@ i_t right_looking_lu_row_permutation_only(const csc_matrix_t& A, toc(factorization_start_time)); last_print = tic(); } - if (toc(factorization_start_time) > settings.time_limit) { - settings.log.printf("Right-looking LU factorization time exceeded\n"); - return -1; - } - + if (toc(start_time) > settings.time_limit) { return TIME_LIMIT_RETURN; } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Concurrent halt\n"); return CONCURRENT_HALT_RETURN; @@ -1235,6 +1235,7 @@ template int right_looking_lu(const csc_matrix_t& A, const simplex_solver_settings_t& settings, double tol, const std::vector& column_list, + double start_time, std::vector& q, csc_matrix_t& L, csc_matrix_t& U, diff --git a/cpp/src/dual_simplex/right_looking_lu.hpp b/cpp/src/dual_simplex/right_looking_lu.hpp index 39182eb3f..5f0bf570b 100644 --- a/cpp/src/dual_simplex/right_looking_lu.hpp +++ b/cpp/src/dual_simplex/right_looking_lu.hpp @@ -19,6 +19,7 @@ i_t right_looking_lu(const csc_matrix_t& A, const simplex_solver_settings_t& settings, f_t tol, const std::vector& column_list, + f_t start_time, std::vector& q, csc_matrix_t& L, csc_matrix_t& U, diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index c49ce2be6..d300d6011 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -154,6 +154,7 @@ lp_status_t solve_linear_program_with_advanced_basis( ok = presolve(original_lp, settings, presolved_lp, presolve_info); } if (ok == CONCURRENT_HALT_RETURN) { return lp_status_t::CONCURRENT_LIMIT; } + if (ok == TIME_LIMIT_RETURN) { return lp_status_t::TIME_LIMIT; } if (ok == -1) { return lp_status_t::INFEASIBLE; } constexpr bool write_out_matlab = false; @@ -349,6 +350,7 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us lp_problem_t presolved_lp(user_problem.handle_ptr, 1, 1, 1); const i_t ok = presolve(original_lp, barrier_settings, presolved_lp, presolve_info); if (ok == CONCURRENT_HALT_RETURN) { return lp_status_t::CONCURRENT_LIMIT; } + if (ok == TIME_LIMIT_RETURN) { return lp_status_t::TIME_LIMIT; } if (ok == -1) { return lp_status_t::INFEASIBLE; } // Apply columns scaling to the presolve LP diff --git a/cpp/src/dual_simplex/types.hpp b/cpp/src/dual_simplex/types.hpp index 9de33ed3b..ea46a1f67 100644 --- a/cpp/src/dual_simplex/types.hpp +++ b/cpp/src/dual_simplex/types.hpp @@ -21,5 +21,7 @@ constexpr float64_t inf = std::numeric_limits::infinity(); // We return this constant to signal that a concurrent halt has occurred #define CONCURRENT_HALT_RETURN -2 +// We return this constant to signal that a time limit has occurred +#define TIME_LIMIT_RETURN -3 } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index baa2b6d2f..7fd8533f8 100644 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu @@ -67,8 +67,12 @@ void rins_t::node_callback(const std::vector& solution, f_t objec // opportunistic early test w/ atomic to avoid having to take the lock if (!rins_thread->cpu_thread_done) return; std::lock_guard lock(rins_mutex); - if (rins_thread->cpu_thread_done && dm.population.current_size() > 0 && - dm.population.is_feasible()) { + bool population_ready = false; + if (rins_thread->cpu_thread_done) { + std::lock_guard pop_lock(dm.population.write_mutex); + population_ready = dm.population.current_size() > 0 && dm.population.is_feasible(); + } + if (population_ready) { lp_optimal_solution = solution; rins_thread->start_cpu_solver(); } @@ -99,8 +103,6 @@ void rins_t::run_rins() { if (total_calls == 0) RAFT_CUDA_TRY(cudaSetDevice(context.handle_ptr->get_device())); - if (!dm.population.is_feasible()) return; - cuopt_assert(lp_optimal_solution.size() == problem_copy->n_variables, "Assignment size mismatch"); cuopt_assert(problem_copy->handle_ptr == &rins_handle, "Handle mismatch"); // Do not make assertions based on problem_ptr. The original problem may have been modified within @@ -111,13 +113,14 @@ void rins_t::run_rins() // "Problem size mismatch"); // cuopt_assert(problem_copy->n_binary_vars == problem_ptr->n_binary_vars, "Problem size // mismatch"); - cuopt_assert(dm.population.current_size() > 0, "No solutions in population"); solution_t best_sol(*problem_copy); rins_handle.sync_stream(); // copy the best from the population into a solution_t in the RINS stream { std::lock_guard lock(dm.population.write_mutex); + if (!dm.population.is_feasible()) return; + cuopt_assert(dm.population.current_size() > 0, "No solutions in population"); auto& best_feasible_ref = dm.population.best_feasible(); cuopt_assert(best_feasible_ref.assignment.size() == best_sol.assignment.size(), "Assignment size mismatch"); diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 425afee46..3e2251171 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -316,9 +316,11 @@ solution_t mip_solver_t::run_solver() if (!is_feasible.value(sol.handle_ptr->get_stream())) { CUOPT_LOG_ERROR( "Solution is not feasible due to variable bounds, returning infeasible solution!"); + context.stats.total_solve_time = timer_.elapsed_time(); context.problem_ptr->post_process_solution(sol); return sol; } + context.stats.total_solve_time = timer_.elapsed_time(); context.problem_ptr->post_process_solution(sol); dm.rins.stop_rins(); return sol; diff --git a/cpp/tests/linear_programming/unit_tests/presolve_test.cu b/cpp/tests/linear_programming/unit_tests/presolve_test.cu index 69fbb17bb..97dfe0957 100644 --- a/cpp/tests/linear_programming/unit_tests/presolve_test.cu +++ b/cpp/tests/linear_programming/unit_tests/presolve_test.cu @@ -166,7 +166,6 @@ TEST(pslp_presolve, postsolve_accuracy_afiro) TEST(pslp_presolve, postsolve_dual_accuracy_afiro) { const raft::handle_t handle_{}; - constexpr double tolerance = 1e-5; auto path = make_path_absolute("linear_programming/afiro_original.mps"); auto mps_data_model = cuopt::mps_parser::parse_mps(path, true); @@ -351,7 +350,6 @@ TEST(pslp_presolve, postsolve_reduced_costs) TEST(pslp_presolve, postsolve_multiple_problems) { const raft::handle_t handle_{}; - constexpr double tolerance = 1e-4; std::vector> instances{ {"afiro_original", -464.75314},