diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 86ed6965c8..5a807de308 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -53,6 +53,7 @@ #define CUOPT_ELIMINATE_DENSE_COLUMNS "eliminate_dense_columns" #define CUOPT_CUDSS_DETERMINISTIC "cudss_deterministic" #define CUOPT_PRESOLVE "presolve" +#define CUOPT_INITIAL_PERTURBATION "initial_perturbation" #define CUOPT_MIP_PROBING "mip_probing" #define CUOPT_DUAL_POSTSOLVE "dual_postsolve" #define CUOPT_MIP_DETERMINISM_MODE "mip_determinism_mode" @@ -198,7 +199,8 @@ #define CUOPT_METHOD_PDLP 1 #define CUOPT_METHOD_DUAL_SIMPLEX 2 #define CUOPT_METHOD_BARRIER 3 -#define CUOPT_METHOD_UNSET 4 +#define CUOPT_METHOD_PRIMAL 4 +#define CUOPT_METHOD_UNSET 5 /* @brief PDLP precision mode constants */ #define CUOPT_PDLP_DEFAULT_PRECISION -1 diff --git a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp index ffcf3fad7a..80197e8a85 100644 --- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp @@ -59,6 +59,7 @@ enum pdlp_solver_mode_t : int { * PDLP: Use the PDLP method. * DualSimplex: Use the dual simplex method. * Barrier: Use the barrier method + * Primal: Use the (experimental) primal simplex method. * Unset: The value was not set. * * @note Default method is Concurrent. @@ -68,6 +69,7 @@ enum method_t : int { PDLP = CUOPT_METHOD_PDLP, DualSimplex = CUOPT_METHOD_DUAL_SIMPLEX, Barrier = CUOPT_METHOD_BARRIER, + Primal = CUOPT_METHOD_PRIMAL, Unset = CUOPT_METHOD_UNSET }; @@ -79,6 +81,7 @@ inline std::string method_to_string(method_t method) case method_t::PDLP: return "PDLP"; case method_t::Barrier: return "Barrier"; case method_t::Concurrent: return "Concurrent"; + case method_t::Primal: return "Primal Simplex"; default: return "Unset"; } } @@ -294,6 +297,7 @@ class pdlp_solver_settings_t { i_t augmented{-1}; i_t dualize{-1}; i_t ordering{-1}; + i_t initial_perturbation{-1}; i_t barrier_dual_initial_point{-1}; i_t postsolve_info{-1}; // Ruiz equilibration for QCQP (barrier) scaling: -1 automatic (row/column diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 4dc6bc67a8..6471071fbd 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -26,8 +26,10 @@ #include #include #include +#include #include #include +#include #include #include @@ -88,6 +90,7 @@ i_t fractional_variables(const simplex_solver_settings_t& settings, { const i_t n = x.size(); assert(x.size() == var_types.size()); + fractional.clear(); for (i_t j = 0; j < n; ++j) { if (is_fractional(x[j], var_types[j], settings.integer_tol)) { fractional.push_back(j); } } @@ -705,10 +708,19 @@ bool branch_and_bound_t::repair_solution(const std::vector& edge_ lp_settings.set_log(false); lp_settings.inside_mip = 2; std::vector leaf_edge_norms = edge_norms; + f_t repair_work_estimate = 0.0; // should probably set the cut off here lp_settings.cut_off - dual_status_t lp_status = simplex::dual_phase2( - 2, 0, lp_start_time, repair_lp, lp_settings, vstatus, lp_solution, iter, leaf_edge_norms); - repaired_solution = lp_solution.x; + dual_status_t lp_status = simplex::dual_phase2(2, + 0, + lp_start_time, + repair_lp, + lp_settings, + vstatus, + lp_solution, + iter, + repair_work_estimate, + leaf_edge_norms); + repaired_solution = lp_solution.x; if (lp_status == dual_status_t::OPTIMAL) { f_t primal_error; @@ -869,6 +881,9 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& exploration_stats_.lexical_reduction_fixings_applied.load(), exploration_stats_.lexical_reduction_pruned_nodes.load()); } + if (integer_pivots_.load() > 0) { + settings_.log.print_format("Number of integer pivots: {}\n", integer_pivots_.load()); + } if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) { solver_status_ = mip_status_t::OPTIMAL; @@ -1604,8 +1619,9 @@ dual_status_t branch_and_bound_t::solve_node_lp( feasible = apply_symmetry_reductions(node_ptr, worker, stats); if (feasible) { - i_t node_iter = 0; - f_t lp_start_time = tic(); + i_t node_iter = 0; + f_t lp_start_time = tic(); + f_t node_work_estimate = 0.0; lp_status = dual_phase2_with_advanced_basis(2, 0, @@ -1619,6 +1635,7 @@ dual_status_t branch_and_bound_t::solve_node_lp( worker->nonbasic_list, worker->leaf_solution, node_iter, + node_work_estimate, worker->leaf_edge_norms); if (lp_status == dual_status_t::NUMERICAL) { @@ -1632,13 +1649,28 @@ dual_status_t branch_and_bound_t::solve_node_lp( worker->basic_list, worker->nonbasic_list, worker->leaf_vstatus, - worker->leaf_edge_norms); + worker->leaf_edge_norms, + node_work_estimate); lp_status = convert_lp_status_to_dual_status(second_status); } stats.total_lp_solve_time += toc(lp_start_time); stats.total_simplex_iters += node_iter; + + if (lp_status == dual_status_t::OPTIMAL) { + std::vector fractional; + i_t num_fractional = + fractional_variables(settings_, worker->leaf_solution.x, var_types_, fractional); + pivot_out_integer_variables(worker->leaf_problem, + worker->basic_list, + worker->nonbasic_list, + worker->leaf_vstatus, + worker->leaf_solution, + worker->basis_factors, + num_fractional, + fractional); + } } } @@ -2778,7 +2810,8 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( basis_update_mpf_t& basis_update, std::vector& basic_list, std::vector& nonbasic_list, - std::vector& edge_norms) + std::vector& edge_norms, + f_t& work_estimate) { lp_status_t root_status; @@ -2794,6 +2827,7 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( nonbasic_list, root_vstatus_, edge_norms_, + work_estimate, nullptr); } @@ -3090,6 +3124,7 @@ auto branch_and_bound_t::do_cut_pass( bool initialize_basis = false; lp_settings.concurrent_halt = NULL; f_t dual_phase2_start_time = tic(); + f_t cut_work_estimate = 0.0; dual_status_t cut_status = dual_phase2_with_advanced_basis(2, 0, initialize_basis, @@ -3102,6 +3137,7 @@ auto branch_and_bound_t::do_cut_pass( nonbasic_list, root_relax_soln_, iter, + cut_work_estimate, edge_norms_); exploration_stats_.total_simplex_iters += iter; f_t dual_phase2_time = toc(dual_phase2_start_time); @@ -3125,7 +3161,8 @@ auto branch_and_bound_t::do_cut_pass( basic_list, nonbasic_list, root_vstatus_, - edge_norms_); + edge_norms_, + cut_work_estimate); if (scratch_status == lp_status_t::OPTIMAL) { // We recovered cut_status = convert_lp_status_to_dual_status(scratch_status); @@ -3141,6 +3178,27 @@ auto branch_and_bound_t::do_cut_pass( } root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + // Refresh fractional info after re-solving with cuts; the pre-cut count is stale. + num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + + pivot_out_integer_variables(original_lp_, + basic_list, + nonbasic_list, + root_vstatus_, + root_relax_soln_, + basis_update, + num_fractional, + fractional); + + dual_degenerate_feasibility_pump(original_lp_, + basic_list, + nonbasic_list, + root_vstatus_, + root_relax_soln_, + basis_update, + num_fractional, + fractional); + if (settings_.benchmark_info_ptr != nullptr) { settings_.benchmark_info_ptr->root_lp_with_cuts = compute_user_objective(original_lp_, root_objective_); @@ -3206,6 +3264,694 @@ auto branch_and_bound_t::do_cut_pass( return {cut_pass_action_t::CONTINUE, mip_status_t::UNSET}; } +template +bool branch_and_bound_t::check_for_dual_degeneracy( + const simplex::lp_solution_t& solution, + const std::vector& nonbasic_list, + std::vector& zero_reduced_costs_vars, + std::vector& zero_reduced_costs_vars_nonbasic_index) +{ + const i_t num_nonbasics = nonbasic_list.size(); + for (i_t k = 0; k < num_nonbasics; k++) { + const i_t j = nonbasic_list[k]; + if (std::abs(solution.z[j]) <= settings_.tight_tol) { + zero_reduced_costs_vars.push_back(j); + zero_reduced_costs_vars_nonbasic_index.push_back(k); + } + } + return !zero_reduced_costs_vars.empty(); +} + +template +void branch_and_bound_t::dual_degenerate_feasibility_pump( + const simplex::lp_problem_t& lp, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& vstatus, + simplex::lp_solution_t& soln, + simplex::basis_update_mpf_t& basis_update, + i_t& num_fractional, + std::vector& fractional) +{ + f_t dual_degenerate_feasibility_pump_start_time = tic(); + std::vector zero_reduced_costs_vars; + std::vector zero_reduced_costs_vars_nonbasic_index; + bool dual_degenerate = check_for_dual_degeneracy( + soln, nonbasic_list, zero_reduced_costs_vars, zero_reduced_costs_vars_nonbasic_index); + if (!dual_degenerate) { return; } + + // Construct a new LP problem + // minimize p^T x + // subject to B x_B + N_z x_z = b - N x_N + // l_B <= x_B <= u_B + // l_z <= x_z <= u_z + // + // where B is the basic matrix, N is the nonbasic matrix, b is the right-hand side, + + const i_t m = lp.num_rows; + const i_t n = lp.num_rows + zero_reduced_costs_vars.size(); + + i_t nnz = 0; + for (i_t j = 0; j < lp.num_cols; j++) { + if (vstatus[j] == variable_status_t::BASIC || std::abs(soln.z[j]) <= settings_.tight_tol) { + nnz += lp.A.col_start[j + 1] - lp.A.col_start[j]; + } + } + simplex::lp_problem_t lp_reduced(lp.handle_ptr, m, n, nnz); + csc_matrix_t& A_reduced = lp_reduced.A; + std::vector original_col_to_reduced_col(lp.num_cols, -1); + i_t nz = 0; + i_t reduced_col = 0; + for (i_t j = 0; j < lp.num_cols; j++) { + if (vstatus[j] == variable_status_t::BASIC || std::abs(soln.z[j]) <= settings_.tight_tol) { + original_col_to_reduced_col[j] = reduced_col; + A_reduced.col_start[reduced_col] = nz; + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + for (i_t p = col_start; p < col_end; p++) { + const i_t i = lp.A.i[p]; + const f_t value = lp.A.x[p]; + A_reduced.i[nz] = i; + A_reduced.x[nz] = value; + nz++; + } + lp_reduced.lower[reduced_col] = lp.lower[j]; + lp_reduced.upper[reduced_col] = lp.upper[j]; + reduced_col++; + } + } + A_reduced.col_start[reduced_col] = nz; + + std::vector b_reduced = lp.rhs; + for (i_t j = 0; j < lp.num_cols; j++) { + if (vstatus[j] == variable_status_t::BASIC || std::abs(soln.z[j]) <= settings_.tight_tol) { + // PASS + } else { + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + for (i_t p = col_start; p < col_end; p++) { + const i_t i = lp.A.i[p]; + const f_t value = lp.A.x[p]; + b_reduced[i] -= value * soln.x[j]; + } + } + } + lp_reduced.rhs = b_reduced; + lp_reduced.obj_scale = 1.0; + + settings_.log.printf( + "Constructed dual degenerate feasibility pump LP with %d rows and %d columns\n", m, n); + + std::vector reduced_basic_list(m); + std::vector reduced_nonbasic_list(zero_reduced_costs_vars.size()); + std::vector reduced_vstatus(n); + i_t num_basic = 0; + i_t num_nonbasic = 0; + reduced_col = 0; + for (i_t j = 0; j < lp.num_cols; j++) { + if (vstatus[j] == variable_status_t::BASIC) { + reduced_vstatus[reduced_col++] = variable_status_t::BASIC; + } else if (std::abs(soln.z[j]) <= 1e-10) { + reduced_nonbasic_list[num_nonbasic++] = + reduced_col; // Does ordering of nonbasic variables matter? + reduced_vstatus[reduced_col++] = vstatus[j]; + } + } + + simplex::lp_solution_t reduced_solution(m, n); + reduced_col = 0; + for (i_t j = 0; j < lp.num_cols; j++) { + if (vstatus[j] == variable_status_t::BASIC || std::abs(soln.z[j]) <= settings_.tight_tol) { + reduced_solution.x[reduced_col++] = soln.x[j]; + } + } + + std::vector reduced_edge_norms(n); + reduced_col = 0; + for (i_t j = 0; j < lp.num_cols; j++) { + if (vstatus[j] == variable_status_t::BASIC || std::abs(soln.z[j]) <= settings_.tight_tol) { + reduced_edge_norms[reduced_col++] = edge_norms_[j]; + } + } + + simplex::basis_update_mpf_t reduced_basis_update = basis_update; + for (i_t k = 0; k < m; k++) { + reduced_basic_list[k] = original_col_to_reduced_col[basic_list[k]]; + } + + f_t primal_work_estimate = 0.0; + i_t iter = 0; + i_t max_pump_iter = 10; + simplex::random_t rng(settings_.random_seed); + i_t best_num_fractional = num_fractional; + std::vector best_reduced_vstatus(n); + bool stalled = false; + for (i_t pump_iter = 0; pump_iter < max_pump_iter; pump_iter++) { + reduced_col = 0; + for (i_t j = 0; j < lp.num_cols; j++) { + if (vstatus[j] == variable_status_t::BASIC || std::abs(soln.z[j]) <= settings_.tight_tol) { + lp_reduced.objective[reduced_col] = 0; + if (var_types_[j] == variable_type_t::INTEGER) { + if (is_fractional( + reduced_solution.x[reduced_col], var_types_[j], settings_.integer_tol)) { + // Default to the exact nearest-integer rounding. Only perturb the + // rounding direction when the previous pass made no progress (a + // zero-pivot solve), to break out of the stall. + const f_t random_value = + stalled ? 0.25 * (2.0 * rng.random() - 1.0) : 0.0; // [-0.25, 0.25] + if (reduced_solution.x[reduced_col] + random_value < + std::floor(reduced_solution.x[reduced_col]) + 0.5) { + lp_reduced.objective[reduced_col] = 1; + } else { + lp_reduced.objective[reduced_col] = -1; + } + } else if (reduced_vstatus[reduced_col] == variable_status_t::NONBASIC_LOWER) { + lp_reduced.objective[reduced_col] = 1; + } else if (reduced_vstatus[reduced_col] == variable_status_t::NONBASIC_UPPER) { + lp_reduced.objective[reduced_col] = -1; + } + } + reduced_col++; + } + } + + bool recompute_basis = false; + const i_t iter_before = iter; + simplex_solver_settings_t primal_settings = settings_; + primal_settings.log.log = false; + primal_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); + primal_settings.work_limit = root_relax_work_estimate_; + simplex::primal_status_t lp_status = + simplex::primal_phase2_with_advanced_basis(2, + exploration_stats_.start_time, + lp_reduced, + primal_settings, + reduced_vstatus, + reduced_basis_update, + reduced_basic_list, + reduced_nonbasic_list, + reduced_solution, + iter, + primal_work_estimate); + // Detect a stall: the solve made no pivots, so the incumbent vertex was + // already optimal for this objective and x did not move. Perturb next pass. + stalled = (iter == iter_before); + + if (lp_status == simplex::primal_status_t::OPTIMAL) { + std::vector adjusted_solution(lp.num_cols, 0.0); + reduced_col = 0; + for (i_t j = 0; j < lp.num_cols; j++) { + if (vstatus[j] == variable_status_t::BASIC || std::abs(soln.z[j]) <= settings_.tight_tol) { + adjusted_solution[j] = reduced_solution.x[reduced_col++]; + } else { + adjusted_solution[j] = soln.x[j]; + } + } + + // Verify the solution is primal feasible + std::vector residual = lp.rhs; + matrix_vector_multiply(lp.A, 1.0, adjusted_solution, -1.0, residual); + const f_t primal_residual = vector_norm_inf(residual); + + if (primal_residual > 1e-6) { + settings_.log.printf("Reduced LP residual|| A*x - b ||_inf = %.4e\n", primal_residual); + } + + std::vector tmp_fractional; + i_t num_fractional_reduced = + fractional_variables(settings_, adjusted_solution, var_types_, tmp_fractional); + settings_.log.printf( + "Degenerate feasibility pump (%d/%d): primal work estimate %.2e, iter %d, fractional " + "variables %d/%d. Time %.2f\n", + pump_iter, + max_pump_iter, + primal_work_estimate, + iter, + num_fractional_reduced, + num_fractional, + toc(dual_degenerate_feasibility_pump_start_time)); + // Also treat a pass that fails to improve the best as a stall, so we perturb + // the next pass even when the solve pivoted (moved) without reducing the count. + stalled = stalled || (num_fractional_reduced >= best_num_fractional); + if (num_fractional_reduced < best_num_fractional) { + best_num_fractional = num_fractional_reduced; + best_reduced_vstatus = reduced_vstatus; + } + } else { + break; + } + } + + settings_.log.printf( + "Degenerate feasibility pump: Simplex iterations %d, Best number of fractional variables " + "%d/%d. Work estimate %.2e, Time %.2f\n", + iter, + best_num_fractional, + num_fractional, + primal_work_estimate, + toc(dual_degenerate_feasibility_pump_start_time)); + if (best_num_fractional < num_fractional) { + // Translate the vstatus from the reduced problem to the vstatus for the original problem + i_t reduced_cols = 0; + for (i_t j = 0; j < lp.num_cols; j++) { + if (vstatus[j] == variable_status_t::BASIC || std::abs(soln.z[j]) <= settings_.tight_tol) { + vstatus[j] = best_reduced_vstatus[reduced_cols++]; + } + } + + std::vector superbasic_list; + nonbasic_list.clear(); + simplex::get_basis_from_vstatus(m, vstatus, basic_list, nonbasic_list, superbasic_list); + assert(superbasic_list.empty()); + const i_t refactor_status = basis_update.refactor_basis(lp.A, + settings_, + lp.lower, + lp.upper, + exploration_stats_.start_time, + basic_list, + nonbasic_list, + vstatus); + if (refactor_status == CONCURRENT_HALT_RETURN || refactor_status == TIME_LIMIT_RETURN) { + // TODO: On failure vstatus, basic_list, and nonbasic_list are in a bad state. + // We should save copies before the failure and restore them after the failure. + return; + } + if (refactor_status != 0) { + settings_.log.printf( + "Failed to refactor basis after dual degenerate feasibility pump. " + "%d deficient columns.\n", + refactor_status); + return; + } + + // Update the solution + // First set the nonbasic variables on their bounds + for (i_t k = 0; k < lp.num_cols - lp.num_rows; k++) { + const i_t j = nonbasic_list[k]; + if (vstatus[j] == variable_status_t::NONBASIC_LOWER || + vstatus[j] == variable_status_t::NONBASIC_FIXED) { + soln.x[j] = lp.lower[j]; + } else if (vstatus[j] == variable_status_t::NONBASIC_UPPER) { + soln.x[j] = lp.upper[j]; + } else { + soln.x[j] = 0; + } + } + // Then compute the effective rhs + std::vector rhs = lp.rhs; + for (i_t j = 0; j < lp.num_cols; j++) { + if (vstatus[j] == variable_status_t::BASIC) { continue; } + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + + const f_t x_j = soln.x[j]; + for (i_t p = col_start; p < col_end; p++) { + const i_t i = lp.A.i[p]; + const f_t aij = lp.A.x[p]; + rhs[i] -= aij * x_j; + } + } + + // Then solve B xB = rhs + std::vector xB(lp.num_rows); + basis_update.b_solve(rhs, xB); + + // Then update the basic variables + for (i_t k = 0; k < lp.num_rows; k++) { + soln.x[basic_list[k]] = xB[k]; + } + + fractional.clear(); + num_fractional = fractional_variables(settings_, soln.x, var_types_, fractional); + } +} + +template +void branch_and_bound_t::apply_delta_x_for_integer_pivot( + const simplex::lp_problem_t& lp, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& nonbasic_index, + std::vector& vstatus, + i_t entering_index, + i_t nonbasic_entering, + i_t direction, + std::vector& delta_x, + const sparse_vector_t& utilde_sparse, + simplex::lp_solution_t& solution, + simplex::basis_update_mpf_t& basis_update, + f_t& work_estimate) +{ + f_t step_length; + i_t basic_leaving; + const i_t leaving_index = simplex::primal_ratio_test(lp, + settings_, + vstatus, + basic_list, + solution.x, + delta_x, + step_length, + basic_leaving, + entering_index, + direction, + work_estimate); + bool binding_integer = + leaving_index != -1 && + is_fractional(solution.x[leaving_index], var_types_[leaving_index], settings_.integer_tol); + if (!binding_integer) { return; } + + std::vector test_x = solution.x; + i_t integer_destroyed = 0; + for (i_t h = 0; h < lp.num_cols; ++h) { + test_x[h] += step_length * delta_x[h]; + if (var_types_[h] != variable_type_t::INTEGER) { continue; } + const bool was_fractional = is_fractional(solution.x[h], var_types_[h], settings_.integer_tol); + const bool now_fractional = is_fractional(test_x[h], var_types_[h], settings_.integer_tol); + if (now_fractional && !was_fractional) { + integer_destroyed++; + } else if (!now_fractional && was_fractional) { + integer_destroyed--; + } + } + // Require a strict net decrease in fractional integers. + if (integer_destroyed >= 0) { return; } + + solution.x = test_x; + basic_list[basic_leaving] = entering_index; + nonbasic_list[nonbasic_entering] = leaving_index; + vstatus[entering_index] = variable_status_t::BASIC; + if (std::abs(lp.upper[leaving_index] - lp.lower[leaving_index]) < 1e-12) { + vstatus[leaving_index] = variable_status_t::NONBASIC_FIXED; + } else if (delta_x[leaving_index] < 0) { + vstatus[leaving_index] = variable_status_t::NONBASIC_LOWER; + } else { + vstatus[leaving_index] = variable_status_t::NONBASIC_UPPER; + } + + // Keep nonbasic_index consistent with nonbasic_list: entering_index is now basic, + // and leaving_index has taken its slot in nonbasic_list. + nonbasic_index[entering_index] = -1; + nonbasic_index[leaving_index] = nonbasic_entering; + + const i_t m = lp.num_rows; + sparse_vector_t es_sparse(m, 1); + es_sparse.i[0] = basic_leaving; + es_sparse.x[0] = 1.0; + sparse_vector_t UTsol_sparse(m, 1); + sparse_vector_t solution_sparse(m, 1); + basis_update.b_transpose_solve(es_sparse, solution_sparse, UTsol_sparse); + const i_t recommend_refactor = basis_update.update(utilde_sparse, UTsol_sparse, basic_leaving); + if (recommend_refactor == 1) { + csc_matrix_t L(m, m, 1); + csc_matrix_t U(m, m, 1); + std::vector pinv(m); + std::vector p(m); + std::vector q(m); + std::vector deficient; + std::vector slacks_needed; + f_t factorize_work_estimate = 0.0; + const i_t rank = factorize_basis(lp.A, + settings_, + basic_list, + exploration_stats_.start_time, + L, + U, + p, + pinv, + q, + deficient, + slacks_needed, + factorize_work_estimate); + if (rank == CONCURRENT_HALT_RETURN || rank == TIME_LIMIT_RETURN) { return; } + if (rank < 0 || rank != lp.num_rows) { return; } + simplex::reorder_basic_list(q, basic_list); + basis_update.reset(L, U, p); + } +} + +template +void branch_and_bound_t::pivot_out_integer_variables( + const simplex::lp_problem_t& lp, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& vstatus, + simplex::lp_solution_t& solution, + simplex::basis_update_mpf_t& basis_update, + i_t& num_fractional, + std::vector& fractional) +{ + if (num_fractional == 0) { return; } + f_t pivot_out_integer_variables_start_time = tic(); + std::vector zero_reduced_costs_vars; + std::vector zero_reduced_costs_vars_nonbasic_index; + bool dual_degenerate = check_for_dual_degeneracy( + solution, nonbasic_list, zero_reduced_costs_vars, zero_reduced_costs_vars_nonbasic_index); + if (!dual_degenerate) { return; } + + lp_solution_t soln_copy = solution; + std::vector basic_list_copy = basic_list; + std::vector nonbasic_list_copy = nonbasic_list; + std::vector vstatus_copy = vstatus; + simplex::basis_update_mpf_t basis_update_copy = basis_update; + + const i_t start_num_fractional = num_fractional; + + const i_t num_zero_reduced_costs_vars = zero_reduced_costs_vars.size(); + + std::vector row_to_slack(lp.num_rows, -1); + for (i_t j : new_slacks_) { + if (lp.lower[j] != 0 || lp.upper[j] != inf) { continue; } + const i_t p = lp.A.col_start[j]; + row_to_slack[lp.A.i[p]] = j; + } + + f_t work_estimate = 0.0; + + std::vector fast_candidates; + std::vector fast_rows; + std::vector fast_nonbasic_slacks; + for (i_t j : fractional) { + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + const i_t num_rows = col_end - col_start; + i_t num_basic_slacks = 0; + i_t num_nonbasic_slacks_with_reduced_cost_zero = 0; + i_t nonbasic_slack = -1; + i_t slack_row = -1; + for (i_t p = col_start; p < col_end; p++) { + const i_t i = lp.A.i[p]; + const i_t slack = row_to_slack[i]; + if (slack >= 0) { + if (vstatus_copy[slack] == variable_status_t::BASIC) { + num_basic_slacks++; + } else if (std::abs(solution.z[slack]) <= 1e-10) { + num_nonbasic_slacks_with_reduced_cost_zero++; + nonbasic_slack = slack; + slack_row = i; + } + } + } + if (num_basic_slacks == num_rows - 1 && num_nonbasic_slacks_with_reduced_cost_zero == 1) { + fast_candidates.push_back(j); + fast_rows.push_back(slack_row); + fast_nonbasic_slacks.push_back(nonbasic_slack); + } + } + + if (fast_candidates.size() > 0) { + settings_.log.printf("Found %ld fast candidates for pivot out integer variables\n", + fast_candidates.size()); + } + + // Build a reverse index nonbasic_index[v] = position of v in nonbasic_list_copy, or -1 if not + // present. Used to locate the entering variable's slot in the fast-candidate path. + // apply_delta_x_for_integer_pivot keeps this index consistent by applying an O(1) fix-up + // on each successful pivot; the two variables whose (non)basic status changes are the only + // entries that need to be updated. + std::vector nonbasic_index(lp.num_cols, -1); + for (i_t p = 0; p < static_cast(nonbasic_list_copy.size()); ++p) { + nonbasic_index[nonbasic_list_copy[p]] = p; + } + + const i_t num_candidates = fast_candidates.size(); + for (i_t k = 0; k < num_candidates; k++) { + const i_t j = fast_candidates[k]; + const i_t row = fast_rows[k]; + const i_t nonbasic_slack = fast_nonbasic_slacks[k]; + // Skip if state changed by a prior successful pivot. + if (vstatus_copy[j] != variable_status_t::BASIC) { continue; } + if (vstatus_copy[nonbasic_slack] == variable_status_t::BASIC) { continue; } + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + f_t a_ij = 0.0; + for (i_t p = col_start; p < col_end; p++) { + const i_t i = lp.A.i[p]; + if (i == row) { + a_ij = lp.A.x[p]; + break; + } + } + f_t bound = a_ij > 0 ? lp.lower[j] : lp.upper[j]; + if (std::abs(bound) == inf) { continue; } + + const f_t delta_xj = bound - soln_copy.x[j]; + const f_t scale = -delta_xj * a_ij; + if (std::abs(scale) <= 1e-12) { continue; } + + // Build delta_x describing "move x[j] to its bound, let the basic slacks compensate to + // keep A*x = b". This is a feasible direction (A*delta_x = 0). The nonzero pattern lives + // on the entries of column A(:, j) plus j itself. In the nonbasic_slack slot, + // delta_x[nonbasic_slack] = -delta_xj * a_ij > 0, i.e. the entering slack moves up from + // its lower bound 0. We build the sparse version to feed the feasibility scan, then + // normalize so that delta_x[nonbasic_slack] == 1 (the convention primal_ratio_test expects + // for entering variables) and scatter into a dense vector. + sparse_vector_t delta_x_sparse; + delta_x_sparse.n = lp.num_cols; + delta_x_sparse.i.reserve(col_end - col_start + 1); + delta_x_sparse.x.reserve(col_end - col_start + 1); + delta_x_sparse.i.push_back(j); + delta_x_sparse.x.push_back(delta_xj); + for (i_t p = col_start; p < col_end; p++) { + const i_t r = lp.A.i[p]; + const f_t a_rj = lp.A.x[p]; + const f_t delta_slack_r = -delta_xj * a_rj; + delta_x_sparse.i.push_back(row_to_slack[r]); + delta_x_sparse.x.push_back(delta_slack_r); + } + + // Reject if the full unit step would drive any basic slack below zero. + bool ok = true; + const i_t ndx = delta_x_sparse.i.size(); + for (i_t h = 0; h < ndx; h++) { + const i_t jj = delta_x_sparse.i[h]; + if (jj == j) continue; + const f_t val = delta_x_sparse.x[h]; + const f_t slack_value = soln_copy.x[jj]; + if (val < -slack_value) { + ok = false; + break; + } + } + if (!ok) { continue; } + + // Normalize so that delta_x[nonbasic_slack] == 1 (the standard entering-direction + // convention). Done on the sparse vector, after the feasibility scan above, which reads + // the unnormalized values. + for (f_t& val : delta_x_sparse.x) { + val /= scale; + } + + std::vector delta_x(lp.num_cols, 0.0); + delta_x_sparse.to_dense(delta_x); + + // Entering variable is the nonbasic slack, moving up from its lower bound 0. + const i_t entering_index = nonbasic_slack; + const i_t nonbasic_entering = nonbasic_index[nonbasic_slack]; + if (nonbasic_entering < 0) { continue; } + const i_t direction = 1; + + // Recover B^{-1} * abar from the full-vector delta_x. In our sign convention, + // delta_x[basic_list[h]] = -direction * (B^{-1} abar)[h], so + // (B^{-1} abar)[h] = -direction * delta_x[basic_list[h]]. + // Then utilde = L^{-1} P abar = U * (B^{-1} abar). In MPF, U == U0 (rank-1 updates all + // live in L), so u_multiply is a single sparse matvec against U0. + std::vector b_inv_abar(lp.num_rows); + for (i_t h = 0; h < lp.num_rows; ++h) { + b_inv_abar[h] = -direction * delta_x[basic_list_copy[h]]; + } + std::vector utilde_dense; + basis_update_copy.u_multiply(b_inv_abar, utilde_dense); + sparse_vector_t utilde_sparse; + utilde_sparse.from_dense(utilde_dense); + + apply_delta_x_for_integer_pivot(lp, + basic_list_copy, + nonbasic_list_copy, + nonbasic_index, + vstatus_copy, + entering_index, + nonbasic_entering, + direction, + delta_x, + utilde_sparse, + soln_copy, + basis_update_copy, + work_estimate); + // apply_delta_x_for_integer_pivot only mutates vstatus when the pivot actually fires, + // so entering_index transitioning to BASIC is a reliable success signal. + if (vstatus_copy[entering_index] == variable_status_t::BASIC) { + settings_.log.printf( + "Fast candidate pivot succeeded: j=%d entering slack=%d row=%d\n", j, entering_index, row); + } + } + + for (i_t k = 0; k < num_zero_reduced_costs_vars; k++) { + const i_t j = zero_reduced_costs_vars[k]; + if (var_types_[j] == variable_type_t::INTEGER) { continue; } + if (vstatus_copy[j] == variable_status_t::BASIC) { continue; } + + const i_t direction = (vstatus_copy[j] == variable_status_t::NONBASIC_LOWER || + vstatus_copy[j] == variable_status_t::NONBASIC_FIXED) + ? 1 + : -1; + const i_t entering_index = j; + const i_t nonbasic_entering = zero_reduced_costs_vars_nonbasic_index[k]; + if (nonbasic_entering < 0 || nonbasic_entering >= static_cast(nonbasic_list_copy.size()) || + nonbasic_list_copy[nonbasic_entering] != j) { + continue; + } + + // Solve B * dxB = A(:, j) so utilde is valid for the MPF update. + // Apply direction when forming delta_x (same convention as primal_phase2). + sparse_vector_t rhs(lp.A, j); + sparse_vector_t delta_xB; + sparse_vector_t utilde_sparse; + basis_update_copy.b_solve(rhs, delta_xB, utilde_sparse); + + std::vector delta_xB_dense; + delta_xB.to_dense(delta_xB_dense); + std::vector delta_x(lp.num_cols, 0.0); + for (i_t h = 0; h < static_cast(basic_list_copy.size()); h++) { + delta_x[basic_list_copy[h]] = -direction * delta_xB_dense[h]; + } + delta_x[j] = direction; + + apply_delta_x_for_integer_pivot(lp, + basic_list_copy, + nonbasic_list_copy, + nonbasic_index, + vstatus_copy, + entering_index, + nonbasic_entering, + direction, + delta_x, + utilde_sparse, + soln_copy, + basis_update_copy, + work_estimate); + } + + std::vector new_fractional; + const i_t num_new_fractional = + fractional_variables(settings_, soln_copy.x, var_types_, new_fractional); + if (num_new_fractional < start_num_fractional) { + i_t num_integer_increased = start_num_fractional - num_new_fractional; + integer_pivots_.fetch_add(num_integer_increased, std::memory_order_release); + settings_.log.printf("Pivoted out %d integer variables: %d -> %d in %.2f\n", + num_integer_increased, + start_num_fractional, + num_new_fractional, + toc(pivot_out_integer_variables_start_time)); + num_fractional = num_new_fractional; + fractional = new_fractional; + basic_list = basic_list_copy; + nonbasic_list = nonbasic_list_copy; + vstatus = vstatus_copy; + basis_update = basis_update_copy; + solution = soln_copy; + } +} + template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { @@ -3287,8 +4033,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut lp_status_t root_status = lp_status_t::UNSET; solving_root_relaxation_ = true; - f_t root_relax_start_time = tic(); - + f_t root_relax_start_time = tic(); + root_relax_work_estimate_ = 0.0; if (!enable_concurrent_lp_root_solve()) { // RINS/SUBMIP path settings_.log.printf("\n"); @@ -3301,7 +4047,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut basic_list, nonbasic_list, root_vstatus_, - edge_norms_); + edge_norms_, + root_relax_work_estimate_); root_relax_solved_by = DualSimplex; exploration_stats_.total_simplex_iters = root_relax_soln_.iterations; @@ -3314,12 +4061,14 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut basis_update, basic_list, nonbasic_list, - edge_norms_); + edge_norms_, + root_relax_work_estimate_); } solving_root_relaxation_ = false; f_t root_relax_elapsed_time = toc(root_relax_start_time); exploration_stats_.total_lp_solve_time = root_relax_elapsed_time; + i_t root_iterations = exploration_stats_.total_simplex_iters; if (root_status == lp_status_t::INFEASIBLE) { settings_.log.printf("\nThe root LP relaxation is infeasible\n", @@ -3373,6 +4122,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_soln_.iterations, root_relax_elapsed_time, method_to_string(root_relax_solved_by)); + settings_.log.printf("Dual simplex iteration %d work estimate %.2e work per second %.2e\n", + root_iterations, + root_relax_work_estimate_, + root_relax_work_estimate_ / root_relax_elapsed_time); settings_.log.printf("Root relaxation objective %+.8e\n\n", root_relax_soln_.user_objective); assert(root_vstatus_.size() == original_lp_.num_cols); @@ -3415,6 +4168,24 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut is_running_ = true; lower_bound_numerical_ = inf; + pivot_out_integer_variables(original_lp_, + basic_list, + nonbasic_list, + root_vstatus_, + root_relax_soln_, + basis_update, + num_fractional, + fractional); + + dual_degenerate_feasibility_pump(original_lp_, + basic_list, + nonbasic_list, + root_vstatus_, + root_relax_soln_, + basis_update, + num_fractional, + fractional); + if (num_fractional != 0 && settings_.max_cut_passes > 0) { print_table_header(); } cut_pool_t cut_pool(original_lp_.num_cols, settings_); @@ -4304,8 +5075,8 @@ node_status_t branch_and_bound_t::solve_node_deterministic( i_t node_iter = 0; f_t lp_start_time = tic(); std::vector leaf_edge_norms = edge_norms_; - - dual_status_t lp_status = dual_phase2_with_advanced_basis(2, + f_t dual_work_estimate = 0.0; + dual_status_t lp_status = dual_phase2_with_advanced_basis(2, 0, worker.recompute_bounds_and_basis, lp_start_time, @@ -4317,6 +5088,7 @@ node_status_t branch_and_bound_t::solve_node_deterministic( worker.nonbasic_list, worker.leaf_solution, node_iter, + dual_work_estimate, leaf_edge_norms, &worker.work_context); @@ -4332,6 +5104,7 @@ node_status_t branch_and_bound_t::solve_node_deterministic( worker.nonbasic_list, worker.leaf_vstatus, leaf_edge_norms, + dual_work_estimate, &worker.work_context); lp_status = convert_lp_status_to_dual_status(second_status); } @@ -4917,6 +5690,7 @@ void branch_and_bound_t::deterministic_dive( worker.leaf_solution.resize(worker.leaf_problem.num_rows, worker.leaf_problem.num_cols); i_t node_iter = 0; f_t lp_start_time = tic(); + f_t dual_work_estimate = 0.0; std::vector leaf_edge_norms = edge_norms_; decompress_vstatus(node_ptr->packed_vstatus, worker.leaf_problem.num_cols, worker.leaf_vstatus); @@ -4932,6 +5706,7 @@ void branch_and_bound_t::deterministic_dive( worker.nonbasic_list, worker.leaf_solution, node_iter, + dual_work_estimate, leaf_edge_norms, &worker.work_context); @@ -4945,6 +5720,7 @@ void branch_and_bound_t::deterministic_dive( worker.nonbasic_list, worker.leaf_vstatus, leaf_edge_norms, + dual_work_estimate, &worker.work_context); lp_status = convert_lp_status_to_dual_status(second_status); } diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 96b8a6d8fe..eaf622b1e3 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -173,7 +173,8 @@ class branch_and_bound_t { simplex::basis_update_mpf_t& basis_update, std::vector& basic_list, std::vector& nonbasic_list, - std::vector& edge_norms); + std::vector& edge_norms, + f_t& work_estimate); i_t find_reduced_cost_fixings(f_t upper_bound, std::vector& lower_bounds, @@ -250,6 +251,7 @@ class branch_and_bound_t { simplex::lp_solution_t root_relax_soln_; simplex::lp_solution_t root_crossover_soln_; method_t root_relax_solved_by{Unset}; + f_t root_relax_work_estimate_; std::vector edge_norms_; std::atomic root_crossover_solution_set_{false}; omp_atomic_t root_lp_current_lower_bound_; @@ -340,6 +342,44 @@ class branch_and_bound_t { i_t leaf_depth, search_strategy_t thread_type); + omp_atomic_t integer_pivots_{0}; + bool check_for_dual_degeneracy(const simplex::lp_solution_t& solution, + const std::vector& nonbasic_list, + std::vector& zero_reduced_costs_vars, + std::vector& zero_reduced_costs_vars_nonbasic_index); + + void pivot_out_integer_variables(const simplex::lp_problem_t& lp, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& vstatus, + simplex::lp_solution_t& soln, + simplex::basis_update_mpf_t& basis_update, + i_t& num_fractional, + std::vector& fractional); + + void apply_delta_x_for_integer_pivot(const simplex::lp_problem_t& lp, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& nonbasic_index, + std::vector& vstatus, + i_t entering_index, + i_t nonbasic_entering, + i_t direction, + std::vector& delta_x, + const sparse_vector_t& utilde_sparse, + simplex::lp_solution_t& solution, + simplex::basis_update_mpf_t& basis_update, + f_t& work_estimate); + + void dual_degenerate_feasibility_pump(const simplex::lp_problem_t& lp, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& vstatus, + simplex::lp_solution_t& soln, + simplex::basis_update_mpf_t& basis_update, + i_t& num_fractional, + std::vector& fractional); + // Repairs low-quality solutions from the heuristics, if it is applicable. void repair_heuristic_solutions(); diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index cdba90f219..eaa60cf475 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -370,6 +370,7 @@ void strong_branch_helper(i_t start, i_t iter = 0; std::vector vstatus = root_vstatus; std::vector child_edge_norms = edge_norms; + f_t child_work_estimate = 0.0; dual_status_t status = simplex::dual_phase2(2, 0, lp_start_time, @@ -378,6 +379,7 @@ void strong_branch_helper(i_t start, vstatus, solution, iter, + child_work_estimate, child_edge_norms); f_t obj = std::numeric_limits::quiet_NaN(); @@ -506,7 +508,8 @@ std::pair trial_branching(const lp_problem_t& orig // Only refactor the basis if we encounter numerical issues. child_basis_factors.set_refactor_frequency(iter_limit); - dual_status_t status = simplex::dual_phase2_with_advanced_basis(2, + f_t child_work_estimate = 0.0; + dual_status_t status = simplex::dual_phase2_with_advanced_basis(2, 0, initialize_basis, start_time, @@ -518,6 +521,7 @@ std::pair trial_branching(const lp_problem_t& orig child_nonbasic_list, solution, iter, + child_work_estimate, child_edge_norms); settings.log.debug("Trial branching on variable %d. Lo: %e Up: %e. Iter %d. Status %s. Obj %e\n", diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index f81962d054..84468ba097 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -2009,6 +2009,34 @@ i_t basis_update_mpf_t::u_solve(sparse_vector_t& rhs) const return 0; } + +// Compute y = U*x. In the MPF factorization, the rank-1 update factors are absorbed into L, so +// U == U0 and U*x reduces to a sparse matvec against U0. +template +void basis_update_mpf_t::u_multiply(const std::vector& x, std::vector& y) const +{ + const i_t m = L0_.m; + y.assign(m, 0.0); + matrix_vector_multiply(U0_, f_t(1.0), x, f_t(0.0), y); + work_estimate_ += 2 * U0_.col_start[U0_.n]; +} + +// Sparse-in/sparse-out overload of u_multiply. Same semantics as the dense version. +template +void basis_update_mpf_t::u_multiply(const sparse_vector_t& x, + sparse_vector_t& y) const +{ + const i_t m = L0_.m; + // Scatter x into a dense workspace, compute U0 * x, gather back to sparse. + std::vector x_dense; + x.to_dense(x_dense); + std::vector y_dense(m, 0.0); + matrix_vector_multiply(U0_, f_t(1.0), x_dense, f_t(0.0), y_dense); + work_estimate_ += 2 * U0_.col_start[U0_.n]; + y.from_dense(y_dense); + work_estimate_ += m; +} + // Solve for x such that L*x = y template i_t basis_update_mpf_t::l_solve(std::vector& rhs) const @@ -2202,7 +2230,7 @@ i_t basis_update_mpf_t::update(const sparse_vector_t& utilde // Ensure the workspace is sorted. Otherwise, the sparse dot will be incorrect. std::sort(xi_workspace_.begin() + m, xi_workspace_.begin() + m + nz, std::less()); - work_estimate_ += (m + nz) * std::log2(m + nz); + work_estimate_ += nz > 1 ? nz * std::log2(nz) : 0; // Gather the workspace into a column of S i_t S_start; diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp index d1c623db55..bdedcc4a18 100644 --- a/cpp/src/dual_simplex/basis_updates.hpp +++ b/cpp/src/dual_simplex/basis_updates.hpp @@ -353,6 +353,14 @@ class basis_update_mpf_t { // Solve for x such that U'*x = y i_t u_transpose_solve(sparse_vector_t& rhs) const; + // Compute y = U*x. In the MPF factorization the rank-1 update factors are absorbed into L, so + // U is unchanged from the initial factorization (U == U0), and U*x is just a sparse matvec + // against U0. + void u_multiply(const std::vector& x, std::vector& y) const; + + // Sparse-in/sparse-out overload of u_multiply. + void u_multiply(const sparse_vector_t& x, sparse_vector_t& y) const; + // Replace the column B(:, leaving_index) with the vector abar. Pass in utilde such that L*utilde // = abar i_t update(const std::vector& utilde, const std::vector& etilde, i_t leaving_index); diff --git a/cpp/src/dual_simplex/bound_flipping_ratio_test.cpp b/cpp/src/dual_simplex/bound_flipping_ratio_test.cpp index cb0964dc05..3fbfbd1f82 100644 --- a/cpp/src/dual_simplex/bound_flipping_ratio_test.cpp +++ b/cpp/src/dual_simplex/bound_flipping_ratio_test.cpp @@ -229,14 +229,14 @@ void bound_flipping_ratio_test_t::heap_passes(const std::vector& }; std::make_heap(bare_idx.begin(), bare_idx.end(), compare); - work_estimate_ += 3 * bare_idx.size(); + work_estimate_ += 10 * bare_idx.size(); while (bare_idx.size() > 0 && slope > 0) { // Remove minimum ratio from the heap and rebalance i_t heap_index = bare_idx.front(); std::pop_heap(bare_idx.begin(), bare_idx.end(), compare); - work_estimate_ += 2 * std::log2(bare_idx.size()); bare_idx.pop_back(); + work_estimate_ += 7 * std::log2(bare_idx.size() + 1); nonbasic_entering = current_indicies[heap_index]; const i_t j = entering_index = nonbasic_list_[nonbasic_entering]; @@ -264,6 +264,7 @@ void bound_flipping_ratio_test_t::heap_passes(const std::vector& // The variable is not bounded. Stop the search. break; } + work_estimate_ += 10; if (toc(start_time_) > settings_.time_limit) { entering_index = RATIO_TEST_TIME_LIMIT; diff --git a/cpp/src/dual_simplex/bound_flipping_ratio_test.hpp b/cpp/src/dual_simplex/bound_flipping_ratio_test.hpp index 2e73d05eff..2f73069451 100644 --- a/cpp/src/dual_simplex/bound_flipping_ratio_test.hpp +++ b/cpp/src/dual_simplex/bound_flipping_ratio_test.hpp @@ -100,7 +100,7 @@ class bound_flipping_ratio_test_t { i_t n_; i_t m_; - f_t work_estimate_; + f_t work_estimate_{0.0}; }; } // namespace cuopt::mathematical_optimization::simplex diff --git a/cpp/src/dual_simplex/crossover.cpp b/cpp/src/dual_simplex/crossover.cpp index e1ba272adf..977f5e5511 100644 --- a/cpp/src/dual_simplex/crossover.cpp +++ b/cpp/src/dual_simplex/crossover.cpp @@ -168,9 +168,10 @@ f_t primal_infeasibility(const lp_problem_t& lp, f_t primal_inf = 0; constexpr bool verbose = false; constexpr f_t infeas_tol = 1e-3; + const f_t primal_tol = settings.primal_tol; for (i_t j = 0; j < n; ++j) { - if (x[j] < lp.lower[j]) { - // x_j < l_j => -x_j > -l_j => -x_j + l_j > 0 + if (x[j] < lp.lower[j] - primal_tol) { + // x_j < l_j - tol => violation exceeds per-variable threshold const f_t infeas = -x[j] + lp.lower[j]; primal_inf += infeas; if (verbose && infeas > infeas_tol) { @@ -183,8 +184,8 @@ f_t primal_infeasibility(const lp_problem_t& lp, vstatus[j]); } } - if (x[j] > lp.upper[j]) { - // x_j > u_j => x_j - u_j > 0 + if (x[j] > lp.upper[j] + primal_tol) { + // x_j > u_j + tol => violation exceeds per-variable threshold const f_t infeas = x[j] - lp.upper[j]; primal_inf += infeas; if (verbose && infeas > infeas_tol) { @@ -1423,8 +1424,11 @@ crossover_status_t crossover(const lp_problem_t& lp, } else if (dual_feasible && !primal_feasible) { i_t dual_iter = 0; std::vector edge_norms; - dual_status_t status = - dual_phase2(2, 0, start_time, lp, settings, vstatus, solution, dual_iter, edge_norms); + f_t work_estimate = 0.0; + simplex_solver_settings_t dual_settings = settings; + dual_settings.iteration_limit = std::numeric_limits::max(); + dual_status_t status = dual_phase2( + 2, 0, start_time, lp, dual_settings, vstatus, solution, dual_iter, work_estimate, edge_norms); if (toc(start_time) > settings.time_limit) { settings.log.printf("Time limit exceeded\n"); return crossover_status_t::TIME_LIMIT; @@ -1443,7 +1447,33 @@ crossover_status_t crossover(const lp_problem_t& lp, solution.iterations += dual_iter; primal_feasible = primal_infeas <= primal_tol && primal_res <= primal_tol; dual_feasible = dual_infeas <= dual_tol && dual_res <= dual_tol; + } else if (primal_feasible && !dual_feasible) { + i_t primal_iter = 0; + simplex_solver_settings_t primal_settings = settings; + primal_settings.iteration_limit = std::numeric_limits::max(); + primal_status_t primal_status = + primal_phase2(2, start_time, lp, primal_settings, vstatus, solution, primal_iter); + if (toc(start_time) > settings.time_limit) { + settings.log.printf("Time limit exceeded\n"); + return crossover_status_t::TIME_LIMIT; + } + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + if (!settings.inside_mip) { settings.log.printf("Concurrent halt\n"); } + return crossover_status_t::CONCURRENT_LIMIT; + } + primal_infeas = primal_infeasibility(lp, settings, vstatus, solution.x); + dual_infeas = dual_infeasibility(lp, settings, vstatus, solution.z); + primal_res = primal_residual(lp, solution); + dual_res = dual_residual(lp, solution); + if (primal_status != primal_status_t::OPTIMAL) { + print_crossover_info(lp, settings, vstatus, solution, "Primal phase 2 complete"); + } + solution.iterations += primal_iter; + primal_feasible = primal_infeas <= primal_tol && primal_res <= primal_tol; + dual_feasible = dual_infeas <= dual_tol && dual_res <= dual_tol; } else { + simplex_solver_settings_t dual_settings = settings; + dual_settings.iteration_limit = std::numeric_limits::max(); lp_problem_t phase1_problem(lp.handle_ptr, 1, 1, 1); create_phase1_problem(lp, phase1_problem); std::vector phase1_vstatus(n); @@ -1469,8 +1499,17 @@ crossover_status_t crossover(const lp_problem_t& lp, i_t iter = 0; lp_solution_t phase1_solution(phase1_problem.num_rows, phase1_problem.num_cols); std::vector junk; - dual_status_t phase1_status = dual_phase2( - 1, 1, start_time, phase1_problem, settings, phase1_vstatus, phase1_solution, iter, junk); + f_t phase1_work_estimate = 0.0; + dual_status_t phase1_status = dual_phase2(1, + 1, + start_time, + phase1_problem, + dual_settings, + phase1_vstatus, + phase1_solution, + iter, + phase1_work_estimate, + junk); if (phase1_status == dual_status_t::NUMERICAL || phase1_status == dual_status_t::DUAL_UNBOUNDED) { settings.log.printf("Failed in Phase 1\n"); @@ -1585,8 +1624,17 @@ crossover_status_t crossover(const lp_problem_t& lp, dual_status_t status = dual_status_t::NUMERICAL; if (dual_infeas <= settings.dual_tol) { std::vector edge_norms; - status = dual_phase2( - 2, iter == 0 ? 1 : 0, start_time, lp, settings, vstatus, solution, iter, edge_norms); + f_t phase2_work_estimate = 0.0; + status = dual_phase2(2, + iter == 0 ? 1 : 0, + start_time, + lp, + dual_settings, + vstatus, + solution, + iter, + phase2_work_estimate, + edge_norms); if (toc(start_time) > settings.time_limit) { settings.log.printf("Time limit exceeded\n"); return crossover_status_t::TIME_LIMIT; diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index a5f10c3229..6e8ef4bbdd 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -160,7 +161,7 @@ void compute_delta_z(const csr_matrix_t& Arow, } } work_estimate += 4 * nz_delta_y; - work_estimate += 4 * nnz_processed; + work_estimate += 5 * nnz_processed; work_estimate += 2 * delta_z_indices.size(); // delta_zB = sigma*ei @@ -469,7 +470,7 @@ void initial_perturbation(const lp_problem_t& lp, f_t sum_perturb = 0.0; i_t num_perturb = 0; - random_t random(settings.seed); + random_t random(settings.random_seed); for (i_t j = 0; j < n; ++j) { f_t obj = objective[j] = lp.objective[j]; @@ -904,7 +905,7 @@ bool update_primal_infeasibilities(const lp_problem_t& lp, primal_inf); if (old_val != 0.0 && squared_infeasibilities[j] == 0.0) { became_feasible = true; } } - work_estimate += 8 * nz; + work_estimate += 9 * nz; return became_feasible; } @@ -1256,6 +1257,7 @@ i_t flip_bounds(const lp_problem_t& lp, num_flipped++; } } + work_estimate += 4 * delta_z_indices.size(); return num_flipped; } @@ -2331,21 +2333,23 @@ void prepare_optimality(i_t info, const simplex_solver_settings_t& settings, basis_update_mpf_t& ft, const std::vector& objective, - const std::vector& basic_list, - const std::vector& nonbasic_list, - const std::vector& vstatus, + // Primal cleanup below pivots, so the basis, the statuses + // and the iteration count are updated in place. + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& vstatus, int phase, f_t start_time, f_t max_val, - i_t iter, + f_t& work_estimate, + i_t& iter, const std::vector& x, std::vector& y, std::vector& z, lp_solution_t& sol) { - const i_t m = lp.num_rows; - const i_t n = lp.num_cols; - f_t work_estimate = 0; // Work in this function is not captured + const i_t m = lp.num_rows; + const i_t n = lp.num_cols; sol.objective = compute_objective(lp, sol.x); sol.user_objective = compute_user_objective(lp, sol.objective); @@ -2367,6 +2371,60 @@ void prepare_optimality(i_t info, perturbation = 0.0; } else { settings.log.printf("Failed to remove perturbation of %.2e.\n", perturbation); + settings.log.printf("Unperturbed dual infeasibility: %.2e\n", dual_infeas); + settings.log.printf("Objective: %+.16e\n", sol.user_objective); + settings.log.printf("Num updates: %d\n", ft.num_updates()); + settings.log.printf("Iterations: %d\n", iter); + + i_t dual_iter = iter; + + // Primal pivots in place, so keep the perturbed solution to fall back on. + // The factor is snapshot rather than refactorized on failure: the copy is + // exact, keeps ft consistent with the restored basis, and cannot itself + // fail the way a refactorization can. + const basis_update_mpf_t saved_ft = ft; + const std::vector saved_x = sol.x; + const std::vector saved_y = sol.y; + const std::vector saved_z = sol.z; + const std::vector saved_vstatus = vstatus; + const std::vector saved_basic_list = basic_list; + const std::vector saved_nonbasic_list = nonbasic_list; + + // Reoptimize the unperturbed objective from this basis. The point is + // primal feasible, so primal simplex stays in phase 2 and pivots only to + // restore dual feasibility. It writes through sol, so x, y and z here see + // the cleaned up solution. It prints no summary; the one below reports the + // final result. + primal_status_t primal_status = primal_phase2_with_advanced_basis(2, + start_time, + lp, + settings, + vstatus, + ft, + basic_list, + nonbasic_list, + sol, + iter, + work_estimate, + false); + if (primal_status == primal_status_t::OPTIMAL) { + // z now prices the original objective, so no perturbation remains. + settings.log.printf("Primal cleanup successful. Iterations %d\n", iter - dual_iter); + perturbation = 0.0; + sol.objective = compute_objective(lp, sol.x); + sol.user_objective = compute_user_objective(lp, sol.objective); + } else { + // Restore the perturbed optimum; a partially pivoted basis is worse than + // the dual feasible point we started from. + settings.log.printf("Primal cleanup failed. Reporting the perturbed solution.\n"); + ft = saved_ft; + sol.x = saved_x; + sol.y = saved_y; + sol.z = saved_z; + vstatus = saved_vstatus; + basic_list = saved_basic_list; + nonbasic_list = saved_nonbasic_list; + } } } } @@ -2379,6 +2437,9 @@ void prepare_optimality(i_t info, settings.log.printf("Dual phase I complete. Iterations %d. Time %.2f\n", iter, toc(start_time)); } if (phase == 2) { + if (settings.inside_mip == 0 || settings.inside_mip == 1) { + settings.log.printf("Work estimate: %.2e\n", work_estimate); + } if (!settings.inside_mip) { settings.log.printf("\n"); settings.log.printf( @@ -2413,6 +2474,21 @@ void prepare_optimality(i_t info, #endif } +template +struct work_timer_t { + work_timer_t(f_t t) : time(t) {} + f_t time{0.0}; + f_t work{0.0}; +}; + +template +work_timer_t& operator+=(work_timer_t& lhs, const work_timer_t& rhs) +{ + lhs.time += rhs.time; + lhs.work += rhs.work; + return lhs; +} + template class phase2_timers_t { public: @@ -2435,60 +2511,89 @@ class phase2_timers_t { { } - void start_timer() + void start_timer(f_t work) { if (!record_time) { return; } start_time = tic(); + start_work = work; + } + + work_timer_t stop_timer(f_t stop_work) + { + if (!record_time) { return work_timer_t(0.0); } + work_timer_t result(toc(start_time)); + result.work = stop_work - start_work; + return result; } - f_t stop_timer() + void print_one(const simplex_solver_settings_t& settings, + const char* name, + const work_timer_t& t, + f_t total_time, + f_t total_work) const { - if (!record_time) { return 0.0; } - return toc(start_time); + const f_t work_per_sec = t.time > 0.0 ? t.work / t.time : f_t(0); + settings.log.printf("%-15s %.2fs %4.1f%% (%.2e work %4.1f%% %.2e/s)\n", + name, + t.time, + total_time > 0.0 ? 100.0 * t.time / total_time : 0.0, + t.work, + total_work > 0.0 ? 100.0 * t.work / total_work : 0.0, + work_per_sec); } void print_timers(const simplex_solver_settings_t& settings) const { if (!record_time) { return; } - const f_t total_time = bfrt_time + pricing_time + btran_time + ftran_time + flip_time + - delta_z_time + lu_update_time + lu_factorization_time + se_norms_time + - se_entering_time + perturb_time + vector_time + objective_time + - update_infeasibility_time; + const f_t total_time = bfrt_time.time + pricing_time.time + btran_time.time + ftran_time.time + + flip_time.time + delta_z_time.time + lu_update_time.time + + lu_factorization_time.time + se_norms_time.time + se_entering_time.time + + perturb_time.time + vector_time.time + objective_time.time + + update_infeasibility_time.time; + const f_t total_work = bfrt_time.work + pricing_time.work + btran_time.work + ftran_time.work + + flip_time.work + delta_z_time.work + lu_update_time.work + + lu_factorization_time.work + se_norms_time.work + se_entering_time.work + + perturb_time.work + vector_time.work + objective_time.work + + update_infeasibility_time.work; // clang-format off - settings.log.printf("BFRT time %.2fs %4.1f%\n", bfrt_time, 100.0 * bfrt_time / total_time); - settings.log.printf("Pricing time %.2fs %4.1f%\n", pricing_time, 100.0 * pricing_time / total_time); - settings.log.printf("BTran time %.2fs %4.1f%\n", btran_time, 100.0 * btran_time / total_time); - settings.log.printf("FTran time %.2fs %4.1f%\n", ftran_time, 100.0 * ftran_time / total_time); - settings.log.printf("Flip time %.2fs %4.1f%\n", flip_time, 100.0 * flip_time / total_time); - settings.log.printf("Delta_z time %.2fs %4.1f%\n", delta_z_time, 100.0 * delta_z_time / total_time); - settings.log.printf("LU update time %.2fs %4.1f%\n", lu_update_time, 100.0 * lu_update_time / total_time); - settings.log.printf("LU factor time %.2fs %4.1f%\n", lu_factorization_time, 100.0 * lu_factorization_time / total_time); - settings.log.printf("SE norms time %.2fs %4.1f%\n", se_norms_time, 100.0 * se_norms_time / total_time); - settings.log.printf("SE enter time %.2fs %4.1f%\n", se_entering_time, 100.0 * se_entering_time / total_time); - settings.log.printf("Perturb time %.2fs %4.1f%\n", perturb_time, 100.0 * perturb_time / total_time); - settings.log.printf("Vector time %.2fs %4.1f%\n", vector_time, 100.0 * vector_time / total_time); - settings.log.printf("Objective time %.2fs %4.1f%\n", objective_time, 100.0 * objective_time / total_time); - settings.log.printf("Inf update time %.2fs %4.1f%\n", update_infeasibility_time, 100.0 * update_infeasibility_time / total_time); - settings.log.printf("Sum %.2fs\n", total_time); + print_one(settings, "BFRT time", bfrt_time, total_time, total_work); + print_one(settings, "Pricing time", pricing_time, total_time, total_work); + print_one(settings, "BTran time", btran_time, total_time, total_work); + print_one(settings, "FTran time", ftran_time, total_time, total_work); + print_one(settings, "Flip time", flip_time, total_time, total_work); + print_one(settings, "Delta_z time", delta_z_time, total_time, total_work); + print_one(settings, "LU update time", lu_update_time, total_time, total_work); + print_one(settings, "LU factor time", lu_factorization_time, total_time, total_work); + print_one(settings, "SE norms time", se_norms_time, total_time, total_work); + print_one(settings, "SE enter time", se_entering_time, total_time, total_work); + print_one(settings, "Perturb time", perturb_time, total_time, total_work); + print_one(settings, "Vector time", vector_time, total_time, total_work); + print_one(settings, "Objective time", objective_time, total_time, total_work); + print_one(settings, "Inf update time", update_infeasibility_time, total_time, total_work); + settings.log.printf("Sum %.2fs (%.2e work %.2e/s)\n", + total_time, + total_work, + total_time > 0.0 ? total_work / total_time : f_t(0)); // clang-format on } - f_t bfrt_time; - f_t pricing_time; - f_t btran_time; - f_t ftran_time; - f_t flip_time; - f_t delta_z_time; - f_t se_norms_time; - f_t se_entering_time; - f_t lu_update_time; - f_t lu_factorization_time; - f_t perturb_time; - f_t vector_time; - f_t objective_time; - f_t update_infeasibility_time; + work_timer_t bfrt_time; + work_timer_t pricing_time; + work_timer_t btran_time; + work_timer_t ftran_time; + work_timer_t flip_time; + work_timer_t delta_z_time; + work_timer_t se_norms_time; + work_timer_t se_entering_time; + work_timer_t lu_update_time; + work_timer_t lu_factorization_time; + work_timer_t perturb_time; + work_timer_t vector_time; + work_timer_t objective_time; + work_timer_t update_infeasibility_time; private: f_t start_time; + f_t start_work; bool record_time; }; @@ -2503,6 +2608,7 @@ dual_status_t dual_phase2(i_t phase, std::vector& vstatus, lp_solution_t& sol, i_t& iter, + f_t& work_estimate, std::vector& delta_y_steepest_edge, work_limit_context_t* work_unit_context) { @@ -2525,6 +2631,7 @@ dual_status_t dual_phase2(i_t phase, nonbasic_list, sol, iter, + work_estimate, delta_y_steepest_edge, work_unit_context); } @@ -2542,6 +2649,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, std::vector& nonbasic_list, lp_solution_t& sol, i_t& iter, + f_t& phase2_work_estimate, std::vector& delta_y_steepest_edge, work_limit_context_t* work_unit_context) { @@ -2556,7 +2664,6 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, assert(lp.lower.size() == n); assert(lp.upper.size() == n); assert(lp.rhs.size() == m); - f_t phase2_work_estimate = 0.0; ft.clear_work_estimate(); std::vector& x = sol.x; @@ -2609,6 +2716,10 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, if (toc(start_time) > settings.time_limit) { return dual_status_t::TIME_LIMIT; } } + if (settings.initial_perturbation == 1 && phase == 2) { + phase2::initial_perturbation(lp, settings, vstatus, objective); + } + // Populate c_basic after basis is initialized for (i_t k = 0; k < m; ++k) { const i_t j = basic_list[k]; @@ -2835,7 +2946,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, i_t basic_leaving_index = -1; i_t leaving_index = -1; f_t max_val; - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); { PHASE2_NVTX_RANGE("DualSimplex::pricing"); if (settings.use_steepest_edge_pricing) { @@ -2856,7 +2967,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, lp, settings, x, basic_list, direction, basic_leaving_index, primal_infeasibility); } } - timers.pricing_time += timers.stop_timer(); + timers.pricing_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); if (leaving_index == -1) { #ifdef CHECK_BASIS_UPDATE for (i_t k = 0; k < basic_list.size(); k++) { @@ -2981,6 +3092,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, phase, start_time, max_val, + phase2_work_estimate, iter, x, y, @@ -2998,7 +3110,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, // BTran // BT*delta_y = -delta_zB = -sigma*ei - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); delta_y_sparse.clear(); UTsol_sparse.clear(); f_t btran_start_work = ft.work_estimate(); @@ -3006,7 +3118,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, PHASE2_NVTX_RANGE("DualSimplex::btran"); phase2::compute_delta_y(ft, basic_leaving_index, direction, delta_y_sparse, UTsol_sparse); } - timers.btran_time += timers.stop_timer(); + timers.btran_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); solve_work += (ft.work_estimate() - btran_start_work); if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { @@ -3030,7 +3142,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, continue; } - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); i_t delta_y_nz0 = 0; const i_t nz_delta_y = delta_y_sparse.i.size(); for (i_t k = 0; k < nz_delta_y; k++) { @@ -3069,7 +3181,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, phase2_work_estimate); } } - timers.delta_z_time += timers.stop_timer(); + timers.delta_z_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return dual_status_t::CONCURRENT_LIMIT; } @@ -3105,7 +3217,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, step_length, nonbasic_entering_index); } else if (bound_flip_ratio) { - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); f_t slope = direction == 1 ? (lp.lower[leaving_index] - x[leaving_index]) : (x[leaving_index] - lp.upper[leaving_index]); bound_flipping_ratio_test_t bfrt(settings, @@ -3128,7 +3240,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, settings.log.printf("Numerical issues encountered in ratio test.\n"); return dual_status_t::NUMERICAL; } - timers.bfrt_time += timers.stop_timer(); + timers.bfrt_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); } else { entering_index = phase2::phase2_ratio_test( lp, settings, vstatus, nonbasic_list, z, delta_z, step_length, nonbasic_entering_index); @@ -3201,6 +3313,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, phase, start_time, max_val, + phase2_work_estimate, iter, x, y, @@ -3257,6 +3370,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, phase, start_time, max_val, + phase2_work_estimate, iter, x, y, @@ -3317,7 +3431,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, return dual_status_t::DUAL_UNBOUNDED; } - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); // Update dual variables // y <- y + steplength * delta_y // z <- z + steplength * delta_z @@ -3333,7 +3447,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, settings.log.printf("Numerical issues encountered in update_dual_variables.\n"); return dual_status_t::NUMERICAL; } - timers.vector_time += timers.stop_timer(); + timers.vector_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); #ifdef COMPUTE_DUAL_RESIDUAL std::vector dual_res1; @@ -3344,7 +3458,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, } #endif - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); // Update primal variable const i_t num_flipped = phase2::flip_bounds(lp, settings, @@ -3361,12 +3475,12 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, atilde_index, phase2_work_estimate); - timers.flip_time += timers.stop_timer(); + timers.flip_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); total_bound_flips += num_flipped; delta_xB_0_sparse.clear(); if (num_flipped > 0) { - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); phase2::adjust_for_flips(ft, basic_list, delta_z_indices, @@ -3378,10 +3492,10 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, delta_x_flip, x, phase2_work_estimate); - timers.ftran_time += timers.stop_timer(); + timers.ftran_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); } - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); utilde_sparse.clear(); scaled_delta_xB_sparse.clear(); rhs_sparse.from_csc_column(lp.A, entering_index); @@ -3408,7 +3522,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, } } solve_work += (ft.work_estimate() - ftran_start_work); - timers.ftran_time += timers.stop_timer(); + timers.ftran_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return dual_status_t::CONCURRENT_LIMIT; } @@ -3420,7 +3534,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, if (primal_step_err > 1e-4) { settings.log.printf("|| A * dx || %e\n", primal_step_err); } #endif - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); f_t se_norms_start_work = ft.work_estimate(); const i_t steepest_edge_status = phase2::update_steepest_edge_norms(settings, basic_list, @@ -3442,18 +3556,18 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, } #endif assert(steepest_edge_status == 0); - timers.se_norms_time += timers.stop_timer(); + timers.se_norms_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); solve_work += (ft.work_estimate() - se_norms_start_work); if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return dual_status_t::CONCURRENT_LIMIT; } - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); // x <- x + delta_x phase2::update_primal_variables( scaled_delta_xB_sparse, basic_list, delta_x, entering_index, x, phase2_work_estimate); - timers.vector_time += timers.stop_timer(); + timers.vector_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); #ifdef COMPUTE_PRIMAL_RESIDUAL residual = lp.rhs; @@ -3464,7 +3578,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, } #endif - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); // TODO(CMM): Do I also need to update the objective due to the bound flips? // TODO(CMM): I'm using the unperturbed objective here, should this be the perturbed objective? phase2::update_objective(basic_list, @@ -3474,9 +3588,9 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, entering_index, obj, phase2_work_estimate); - timers.objective_time += timers.stop_timer(); + timers.objective_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); // Update primal infeasibilities due to changes in basic variables // from flipping bounds #ifdef CHECK_BASIC_INFEASIBILITIES @@ -3529,17 +3643,17 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, phase2::check_primal_infeasibilities( lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices); #endif - timers.update_infeasibility_time += timers.stop_timer(); + timers.update_infeasibility_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); // Clear delta_x phase2::clear_delta_x( basic_list, entering_index, scaled_delta_xB_sparse, delta_x, phase2_work_estimate); - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); f_t sum_perturb = 0.0; phase2::compute_perturbation( lp, settings, delta_z_indices, z, objective, sum_perturb, phase2_work_estimate); - timers.perturb_time += timers.stop_timer(); + timers.perturb_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); // Update basis information vstatus[entering_index] = variable_status_t::BASIC; @@ -3562,7 +3676,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 5); #endif - timers.start_timer(); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); // Refactor or update the basis factorization { PHASE2_NVTX_RANGE("DualSimplex::basis_update"); @@ -3578,8 +3692,8 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, phase2::check_update(lp, settings, ft, basic_list, basic_leaving_index); #endif should_refactor = recommend_refactor == 1; - timers.lu_update_time += timers.stop_timer(); - timers.start_timer(); + timers.lu_update_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); + timers.start_timer(phase2_work_estimate + ft.work_estimate()); } #ifdef CHECK_BASIC_INFEASIBILITIES @@ -3657,7 +3771,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 7); #endif } - timers.lu_factorization_time += timers.stop_timer(); + timers.lu_factorization_time += timers.stop_timer(phase2_work_estimate + ft.work_estimate()); #ifdef STEEPEST_EDGE_DEBUG if (iter < 100 || iter % 100 == 0)) @@ -3737,10 +3851,6 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, 100.0 * dense_delta_z / (sparse_delta_z + dense_delta_z)); ft.print_stats(); } - if (settings.inside_mip == 1 && settings.concurrent_halt != nullptr) { - settings.log.debug("Setting concurrent halt in Dual Simplex Phase 2\n"); - *settings.concurrent_halt = 1; - } } return status; } @@ -3756,6 +3866,7 @@ template dual_status_t dual_phase2( std::vector& vstatus, lp_solution_t& sol, int& iter, + double& work_estimate, std::vector& steepest_edge_norms, work_limit_context_t* work_unit_context); @@ -3772,6 +3883,7 @@ template dual_status_t dual_phase2_with_advanced_basis( std::vector& nonbasic_list, lp_solution_t& sol, int& iter, + double& work_estimate, std::vector& steepest_edge_norms, work_limit_context_t* work_unit_context); diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp index daa946e019..e5a4bacf62 100644 --- a/cpp/src/dual_simplex/phase2.hpp +++ b/cpp/src/dual_simplex/phase2.hpp @@ -60,6 +60,7 @@ dual_status_t dual_phase2(i_t phase, std::vector& vstatus, lp_solution_t& sol, i_t& iter, + f_t& work_estimate, std::vector& steepest_edge_norms, work_limit_context_t* work_unit_context = nullptr); @@ -76,6 +77,7 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, std::vector& nonbasic_list, lp_solution_t& sol, i_t& iter, + f_t& work_estimate, std::vector& delta_y_steepest_edge, work_limit_context_t* work_unit_context = nullptr); diff --git a/cpp/src/dual_simplex/primal.cpp b/cpp/src/dual_simplex/primal.cpp index 78c7107ca3..8867c8c9b4 100644 --- a/cpp/src/dual_simplex/primal.cpp +++ b/cpp/src/dual_simplex/primal.cpp @@ -14,18 +14,128 @@ #include #include +#include + namespace cuopt::mathematical_optimization::simplex { +template +struct primal_work_timer_t { + primal_work_timer_t(f_t t) : time(t) {} + f_t time{0.0}; + f_t work{0.0}; +}; + +template +primal_work_timer_t& operator+=(primal_work_timer_t& lhs, + const primal_work_timer_t& rhs) +{ + lhs.time += rhs.time; + lhs.work += rhs.work; + return lhs; +} + +template +class primal_timers_t { + public: + primal_timers_t(bool should_time) + : record_time(should_time), + pricing_time(0), + ftran_time(0), + ratio_test_time(0), + btran_time(0), + delta_z_time(0), + update_duals_time(0), + lu_update_time(0), + lu_factorization_time(0), + update_x_time(0) + { + } + + void start_timer(f_t work) + { + if (!record_time) { return; } + start_time_ = tic(); + start_work_ = work; + } + + primal_work_timer_t stop_timer(f_t stop_work) + { + if (!record_time) { return primal_work_timer_t(0.0); } + primal_work_timer_t result(toc(start_time_)); + result.work = stop_work - start_work_; + return result; + } + + void print_one(const simplex_solver_settings_t& settings, + const char* name, + const primal_work_timer_t& t, + f_t total_time, + f_t total_work) const + { + const f_t work_per_sec = t.time > 0.0 ? t.work / t.time : f_t(0); + settings.log.printf("%-15s %.2fs %4.1f%% (%.2e work %4.1f%% %.2e/s)\n", + name, + t.time, + total_time > 0.0 ? 100.0 * t.time / total_time : 0.0, + t.work, + total_work > 0.0 ? 100.0 * t.work / total_work : 0.0, + work_per_sec); + } + + void print_timers(const simplex_solver_settings_t& settings) const + { + if (!record_time) { return; } + const f_t total_time = pricing_time.time + ftran_time.time + ratio_test_time.time + + btran_time.time + delta_z_time.time + update_duals_time.time + + lu_update_time.time + lu_factorization_time.time + update_x_time.time; + const f_t total_work = pricing_time.work + ftran_time.work + ratio_test_time.work + + btran_time.work + delta_z_time.work + update_duals_time.work + + lu_update_time.work + lu_factorization_time.work + update_x_time.work; + // clang-format off + print_one(settings, "Pricing time", pricing_time, total_time, total_work); + print_one(settings, "FTran time", ftran_time, total_time, total_work); + print_one(settings, "Ratio test", ratio_test_time, total_time, total_work); + print_one(settings, "BTran time", btran_time, total_time, total_work); + print_one(settings, "Delta_z time", delta_z_time, total_time, total_work); + print_one(settings, "Update duals", update_duals_time, total_time, total_work); + print_one(settings, "LU update time", lu_update_time, total_time, total_work); + print_one(settings, "LU factor time", lu_factorization_time, total_time, total_work); + print_one(settings, "Update x time", update_x_time, total_time, total_work); + settings.log.printf("Sum %.2fs (%.2e work %.2e/s)\n", + total_time, + total_work, + total_time > 0.0 ? total_work / total_time : f_t(0)); + // clang-format on + } + + primal_work_timer_t pricing_time; + primal_work_timer_t ftran_time; + primal_work_timer_t ratio_test_time; + primal_work_timer_t btran_time; + primal_work_timer_t delta_z_time; + primal_work_timer_t update_duals_time; + primal_work_timer_t lu_update_time; + primal_work_timer_t lu_factorization_time; + primal_work_timer_t update_x_time; + + private: + f_t start_time_; + f_t start_work_; + bool record_time; +}; + namespace { template void set_primal_variables_on_bounds(const lp_problem_t& lp, const simplex_solver_settings_t& settings, - const std::vector& z, std::vector& vstatus, - std::vector& x) + std::vector& x, + f_t& work_estimate) { - const i_t n = lp.num_cols; + const i_t m = lp.num_rows; + const i_t n = lp.num_cols; + constexpr f_t diff_tol = 1e-6; for (i_t j = 0; j < n; ++j) { if (vstatus[j] == variable_status_t::BASIC) { continue; } @@ -53,18 +163,20 @@ void set_primal_variables_on_bounds(const lp_problem_t& lp, assert(1 == 0); } } + work_estimate += n + 3.0 * (n - m); } template f_t dual_infeasibility(const lp_problem_t& lp, const std::vector& vstatus, - const std::vector& z) + const std::vector& z, + f_t tight_tol, + i_t& num_infeasible, + f_t& work_estimate) { const i_t n = lp.num_cols; - const i_t m = lp.num_rows; - i_t num_infeasible = 0; + num_infeasible = 0; f_t sum_infeasible = 0.0; - constexpr f_t tight_tol = 0; i_t lower_bound_inf = 0; i_t upper_bound_inf = 0; i_t free_inf = 0; @@ -102,6 +214,7 @@ f_t dual_infeasibility(const lp_problem_t& lp, non_basic_upper_inf++; } } + work_estimate += 8 * n; return sum_infeasible; } @@ -111,9 +224,11 @@ i_t phase2_pricing(const lp_problem_t& lp, const std::vector& z, const std::vector& nonbasic_list, const std::vector& vstatus, + f_t dual_tol, i_t& direction, i_t& basic_entering, - f_t& dual_inf) + f_t& dual_inf, + f_t& work_estimate) { const i_t m = lp.num_rows; const i_t n = lp.num_cols; @@ -121,8 +236,7 @@ i_t phase2_pricing(const lp_problem_t& lp, f_t max_infeas = 0.0; dual_inf = 0.0; for (i_t k = 0; k < n - m; ++k) { - const i_t j = nonbasic_list[k]; - constexpr f_t dual_tol = 1e-6; + const i_t j = nonbasic_list[k]; if (vstatus[j] == variable_status_t::NONBASIC_FIXED) { continue; } if ((vstatus[j] == variable_status_t::NONBASIC_LOWER || vstatus[j] == variable_status_t::NONBASIC_FREE) && @@ -148,69 +262,30 @@ i_t phase2_pricing(const lp_problem_t& lp, } } } + work_estimate += 5 * (n - m); return entering_index; } -template -i_t ratio_test(const lp_problem_t& lp, - const std::vector& vstatus, - const std::vector& basic_list, - std::vector& x, - std::vector& delta_x, - f_t& step_length, - i_t& basic_leaving) -{ - const i_t m = lp.num_rows; - const i_t n = lp.num_cols; - basic_leaving = -1; - i_t leaving_index = -1; - f_t min_val = inf; - constexpr f_t pivot_tol = 1e-8; - for (i_t k = 0; k < m; ++k) { - const i_t j = basic_list[k]; - if (delta_x[j] == 0.0) { continue; } - if (lp.lower[j] > -inf && x[j] >= lp.lower[j] && delta_x[j] < -pivot_tol) { - // xj + step * delta_x[j] >= lp.lower[j] - // step * delta_x[j] >= lp.lower[j] - x[j] - // step <= (lp.lower[j] - x[j]) / delta_x[j], delta_x[j] < 0 - const f_t neum = lp.lower[j] - x[j]; - f_t ratio = neum / delta_x[j]; - if (ratio < min_val) { - min_val = ratio; - basic_leaving = k; - leaving_index = j; - } - } - if (lp.upper[j] < inf && x[j] <= lp.upper[j] && delta_x[j] > pivot_tol) { - // xj + step * delta_x[j] <= lp.upper[j] - // step * delta_x[j] <= lp.upper[j] - x[j] - // step <= (lp.upper[j] - x[j]) / delta_x[j], delta_x[j] > 0 - const f_t neum = lp.upper[j] - x[j]; - f_t ratio = neum / delta_x[j]; - if (ratio < min_val) { - min_val = ratio; - basic_leaving = k; - leaving_index = j; - } - } - } - step_length = min_val; - return leaving_index; -} - template f_t primal_infeasibility(const lp_problem_t& lp, const simplex_solver_settings_t& settings, const std::vector& vstatus, - const std::vector& x) + const std::vector& x, + i_t& num_infeasible, + f_t& work_estimate) { + const i_t m = lp.num_rows; const i_t n = lp.num_cols; f_t primal_inf = 0; + num_infeasible = 0; for (i_t j = 0; j < n; ++j) { - if (x[j] < lp.lower[j]) { + // Nonbasics are pinned to a bound; only basics can be (legitimately) infeasible. + if (vstatus[j] != variable_status_t::BASIC) { continue; } + if (x[j] < lp.lower[j] - settings.primal_tol) { // x_j < l_j => -x_j > -l_j => -x_j + l_j > 0 const f_t infeas = -x[j] + lp.lower[j]; primal_inf += infeas; + num_infeasible++; if (infeas > 1e-6) { settings.log.debug("x %d infeas %e lo %e val %e up %e vstatus %hhd\n", j, @@ -221,10 +296,11 @@ f_t primal_infeasibility(const lp_problem_t& lp, vstatus[j]); } } - if (x[j] > lp.upper[j]) { + if (x[j] > lp.upper[j] + settings.primal_tol) { // x_j > u_j => x_j - u_j > 0 const f_t infeas = x[j] - lp.upper[j]; primal_inf += infeas; + num_infeasible++; if (infeas > 1e-6) { settings.log.debug("x %d infeas %e lo %e val %e up %e vstatus %hhd\n", j, @@ -236,15 +312,366 @@ f_t primal_infeasibility(const lp_problem_t& lp, } } } + work_estimate += n + 4 * m; return primal_inf; } +template +f_t primal_infeasibility(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& vstatus, + const std::vector& x, + f_t& work_estimate) +{ + i_t num_infeasible = 0; + return primal_infeasibility(lp, settings, vstatus, x, num_infeasible, work_estimate); +} + +// work estimate: n-m + 4 * m +template +void compute_phase1_objective(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& vstatus, + const std::vector& x, + std::vector& objective, + f_t& work_estimate) +{ + const i_t m = lp.num_rows; + const i_t n = lp.num_cols; + for (i_t j = 0; j < n; ++j) { + if (vstatus[j] != variable_status_t::BASIC) { + objective[j] = 0.0; + } else if (x[j] < lp.lower[j] - settings.primal_tol) { + objective[j] = -1.0; + } else if (x[j] > lp.upper[j] + settings.primal_tol) { + objective[j] = 1.0; + } else { + objective[j] = 0.0; + } + } + work_estimate += n - m + 4 * m; +} + +template +void compute_delta_y(const basis_update_mpf_t& basis_update, + i_t basic_leaving, + sparse_vector_t& delta_y, + sparse_vector_t& etilde) +{ + const i_t m = delta_y.n; + sparse_vector_t ei(m, 1); + ei.i[0] = basic_leaving; + ei.x[0] = 1.0; + delta_y.clear(); + etilde.clear(); + basis_update.b_transpose_solve(ei, delta_y, etilde); +} + +template +void compute_delta_z(const csr_matrix_t& Arow, + const std::vector& vstatus, + const sparse_vector_t& delta_y, + std::vector& delta_z, + f_t& work_estimate) +{ + // A^T delta_y + delta_z = 0 + // delta_z = -A^T delta_y = - sum_i A(i, :) * delta_y_i + std::fill(delta_z.begin(), delta_z.end(), 0.0); + work_estimate += delta_z.size(); + for (i_t k = 0; k < static_cast(delta_y.i.size()); ++k) { + const i_t i = delta_y.i[k]; + const f_t delta_y_i = delta_y.x[k]; + const i_t row_start = Arow.row_start[i]; + const i_t row_end = Arow.row_start[i + 1]; + for (i_t p = row_start; p < row_end; ++p) { + const i_t j = Arow.j[p]; + if (vstatus[j] != variable_status_t::BASIC) { delta_z[j] -= Arow.x[p] * delta_y_i; } + } + work_estimate += 5 * (row_end - row_start); + } + work_estimate += 4 * delta_y.i.size(); +} + +template +f_t compute_dual_step_length(f_t entering_reduced_cost, f_t pivot) +{ + assert(pivot != 0.0); + return entering_reduced_cost / pivot; +} + +template +void update_y(f_t dual_step_length, + const sparse_vector_t& delta_y, + std::vector& y, + f_t& work_estimate) +{ + for (i_t k = 0; k < static_cast(delta_y.i.size()); ++k) { + const i_t i = delta_y.i[k]; + y[i] += dual_step_length * delta_y.x[k]; + } + work_estimate += 3 * delta_y.i.size(); +} + +template +void update_z(f_t dual_step_length, + const std::vector& nonbasic_list, + i_t entering_index, + const std::vector& delta_z, + std::vector& z, + f_t& work_estimate) +{ + for (i_t k = 0; k < static_cast(nonbasic_list.size()); ++k) { + const i_t j = nonbasic_list[k]; + z[j] += dual_step_length * delta_z[j]; + } + work_estimate += 3 * nonbasic_list.size(); + z[entering_index] = 0.0; +} + +template +void compute_dual_variables(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& objective, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& ft, + std::vector& c_basic, + std::vector& y, + std::vector& z, + f_t& work_estimate) +{ + const i_t m = lp.num_rows; + const i_t n = lp.num_cols; + // Solve for y such that B'*y = c_B + for (i_t k = 0; k < m; ++k) { + const i_t j = basic_list[k]; + c_basic[k] = objective[j]; + } + work_estimate += 3 * m; + ft.b_transpose_solve(c_basic, y); + // zN = cN - N'*y + for (i_t k = 0; k < n - m; k++) { + const i_t j = nonbasic_list[k]; + // z_j <- c_j + z[j] = objective[j]; + + // z_j <- z_j - A(:, j)'*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] * y[lp.A.i[p]]; + } + work_estimate += 3.0 * (col_end - col_start); + z[j] -= dot; + } + work_estimate += 6 * (n - m); + // zB = 0 + for (i_t k = 0; k < m; ++k) { + z[basic_list[k]] = 0.0; + } + work_estimate += 2 * m; +} + +template +void compute_basic_primal_variables(const lp_problem_t& lp, + const basis_update_mpf_t& basis_update, + const std::vector& basic_list, + const std::vector& nonbasic_list, + std::vector& x, + f_t& work_estimate) +{ + const i_t m = lp.num_rows; + const i_t n = lp.num_cols; + std::vector rhs = lp.rhs; + for (i_t k = 0; k < n - m; ++k) { + const i_t j = nonbasic_list[k]; + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + const f_t xj = x[j]; + for (i_t p = col_start; p < col_end; ++p) { + rhs[lp.A.i[p]] -= xj * lp.A.x[p]; + } + work_estimate += 3.0 * (col_end - col_start); + } + work_estimate += 4 * (n - m); + std::vector xB(m); + work_estimate += m; + basis_update.b_solve(rhs, xB); + for (i_t k = 0; k < m; ++k) { + x[basic_list[k]] = xB[k]; + } + work_estimate += 3 * m; +} + +template +f_t primal_constraint_residual(const lp_problem_t& lp, const std::vector& x) +{ + std::vector residual = lp.rhs; + matrix_vector_multiply(lp.A, 1.0, x, -1.0, residual); + return vector_norm_inf(residual); +} + } // namespace -// Note this implementation of primal simplex is experimental -// It is meant only to serve as a method to remove the perturbation to the objective -// after dual simplex has found a primal feasible solution -// The implementation currently cycles. So is not enabled at this time. +template +i_t primal_ratio_test(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& vstatus, + const std::vector& basic_list, + std::vector& x, + std::vector& delta_x, + f_t& step_length, + i_t& basic_leaving, + i_t entering_index, + i_t direction, + f_t& work_estimate) +{ + const i_t m = lp.num_rows; + basic_leaving = -1; + i_t leaving_index = -1; + constexpr f_t pivot_tol = 1e-8; + constexpr f_t harris_tol = 1e-8; + + // Harris ratio test: two passes. + // Pass 1: find the maximum step length alpha_1 such that no variable + // moves more than harris_tol past its bound. + // Pass 2: among all candidates with ratio <= alpha_1, pick the one + // with the largest pivot (|delta_x[j]|). + + f_t alpha_1 = inf; + + // Entering variable can hit its opposite bound: limit step by that + if (direction > 0 && lp.upper[entering_index] < inf) { + const f_t limit = lp.upper[entering_index] - x[entering_index]; + if (limit >= 0 && limit < alpha_1) { alpha_1 = limit; } + } else if (direction < 0 && lp.lower[entering_index] > -inf) { + const f_t limit = x[entering_index] - lp.lower[entering_index]; + if (limit >= 0 && limit < alpha_1) { alpha_1 = limit; } + } + + // Pass 1: compute alpha_1 (Harris step) + for (i_t k = 0; k < m; ++k) { + const i_t j = basic_list[k]; + if (std::abs(delta_x[j]) <= pivot_tol) { continue; } + + // Already below lower and moving back up: stop exactly at the bound. + // No harris tolerance here — these variables are already infeasible + // and must not overshoot their bound (needed for Phase I correctness). + if (x[j] < lp.lower[j] && delta_x[j] > pivot_tol && lp.lower[j] > -inf) { + const f_t ratio = (lp.lower[j] - x[j]) / delta_x[j]; + if (ratio >= 0 && ratio < alpha_1) { alpha_1 = ratio; } + } + // Already above upper and moving back down: stop exactly at the bound. + if (x[j] > lp.upper[j] && delta_x[j] < -pivot_tol && lp.upper[j] < inf) { + const f_t ratio = (lp.upper[j] - x[j]) / delta_x[j]; + if (ratio >= 0 && ratio < alpha_1) { alpha_1 = ratio; } + } + + if (lp.lower[j] > -inf && delta_x[j] < -pivot_tol) { + // xj + step * delta_x[j] >= lp.lower[j] - harris_tol + f_t neum = lp.lower[j] - x[j] - harris_tol; + if (neum > 0 && neum <= settings.primal_tol) { neum = 0.0; } + f_t ratio = neum / delta_x[j]; + if (ratio >= 0 && ratio < alpha_1) { alpha_1 = ratio; } + } + if (lp.upper[j] < inf && delta_x[j] > pivot_tol) { + // xj + step * delta_x[j] <= lp.upper[j] + harris_tol + f_t neum = lp.upper[j] - x[j] + harris_tol; + if (neum < 0 && -neum <= settings.primal_tol) { neum = 0.0; } + f_t ratio = neum / delta_x[j]; + if (ratio >= 0 && ratio < alpha_1) { alpha_1 = ratio; } + } + } + + // Pass 2: among candidates with exact ratio <= alpha_1, pick largest pivot + f_t best_pivot = 0.0; + step_length = alpha_1; + + // Check entering variable bound (no pivot selection needed — it's fixed at direction) + if (direction > 0 && lp.upper[entering_index] < inf) { + const f_t limit = lp.upper[entering_index] - x[entering_index]; + if (limit >= 0 && limit <= alpha_1) { + // Entering hits its own bound — this is always pivot = 1.0 effectively + step_length = limit; + leaving_index = -1; + basic_leaving = -1; + best_pivot = inf; // Always prefer this if it's within alpha_1 + } + } else if (direction < 0 && lp.lower[entering_index] > -inf) { + const f_t limit = x[entering_index] - lp.lower[entering_index]; + if (limit >= 0 && limit <= alpha_1) { + step_length = limit; + leaving_index = -1; + basic_leaving = -1; + best_pivot = inf; + } + } + + for (i_t k = 0; k < m; ++k) { + const i_t j = basic_list[k]; + if (std::abs(delta_x[j]) <= pivot_tol) { continue; } + + const f_t abs_dx = std::abs(delta_x[j]); + + // Already below lower and moving back up: stop when we reach the lower bound. + // Without this, phase I can take an unbounded step (false unbounded) or skip the + // breakpoint of the piecewise phase-I objective and stall still infeasible. + if (x[j] < lp.lower[j] && delta_x[j] > pivot_tol && lp.lower[j] > -inf) { + const f_t ratio = (lp.lower[j] - x[j]) / delta_x[j]; + if (ratio >= 0 && ratio <= alpha_1 && abs_dx > best_pivot) { + best_pivot = abs_dx; + step_length = ratio; + basic_leaving = k; + leaving_index = j; + } + } + // Already above upper and moving back down + if (x[j] > lp.upper[j] && delta_x[j] < -pivot_tol && lp.upper[j] < inf) { + const f_t ratio = (lp.upper[j] - x[j]) / delta_x[j]; + if (ratio >= 0 && ratio <= alpha_1 && abs_dx > best_pivot) { + best_pivot = abs_dx; + step_length = ratio; + basic_leaving = k; + leaving_index = j; + } + } + + if (lp.lower[j] > -inf && delta_x[j] < -pivot_tol) { + // xj + step * delta_x[j] >= lp.lower[j] + // step <= (lp.lower[j] - x[j]) / delta_x[j], delta_x[j] < 0 + f_t neum = lp.lower[j] - x[j]; + // A basic sitting below its bound (within the primal tolerance) is on + // the bound numerically. Treat it as a zero-length block. + if (neum > 0 && neum <= settings.primal_tol) { neum = 0.0; } + f_t ratio = neum / delta_x[j]; + if (ratio >= 0 && ratio <= alpha_1 && abs_dx > best_pivot) { + best_pivot = abs_dx; + step_length = ratio; + basic_leaving = k; + leaving_index = j; + } + } + if (lp.upper[j] < inf && delta_x[j] > pivot_tol) { + // xj + step * delta_x[j] <= lp.upper[j] + // step <= (lp.upper[j] - x[j]) / delta_x[j], delta_x[j] > 0 + f_t neum = lp.upper[j] - x[j]; + // Mirror of the lower bound case: slightly above the bound is on the bound. + if (neum < 0 && -neum <= settings.primal_tol) { neum = 0.0; } + f_t ratio = neum / delta_x[j]; + if (ratio >= 0 && ratio <= alpha_1 && abs_dx > best_pivot) { + best_pivot = abs_dx; + step_length = ratio; + basic_leaving = k; + leaving_index = j; + } + } + } + + work_estimate += 10 * m; + return leaving_index; +} + template primal_status_t primal_phase2(i_t phase, f_t start_time, @@ -256,34 +683,14 @@ primal_status_t primal_phase2(i_t phase, { const i_t m = lp.num_rows; const i_t n = lp.num_cols; - assert(m <= n); - assert(vstatus.size() == n); - assert(lp.A.m == m); - assert(lp.A.n == n); - assert(lp.objective.size() == n); - assert(lp.lower.size() == n); - assert(lp.upper.size() == n); - assert(lp.rhs.size() == m); + f_t work_estimate = 0; std::vector basic_list(m); std::vector nonbasic_list; std::vector superbasic_list; - std::vector bound_info(n - m); - - std::vector& x = sol.x; - std::vector& y = sol.y; - std::vector& z = sol.z; - - std::vector incoming_x = x; - std::vector incoming_vstatus = vstatus; - - settings.log.printf("Primal Simplex Phase %d\n", phase); - settings.log.printf("Solving a problem with %d constraints %d variables %d nonzeros\n", - lp.num_rows, - lp.num_cols, - lp.A.col_start[lp.num_cols]); get_basis_from_vstatus(m, vstatus, basic_list, nonbasic_list, superbasic_list); + work_estimate += 2 * n; assert(superbasic_list.size() == 0); assert(nonbasic_list.size() == n - m); @@ -308,6 +715,7 @@ primal_status_t primal_phase2(i_t phase, slacks_needed, work_estimate); if (rank == CONCURRENT_HALT_RETURN) { + settings.log.printf("Concurrent halt in primal phase2\n"); return primal_status_t::CONCURRENT_LIMIT; } else if (rank == TIME_LIMIT_RETURN) { return primal_status_t::TIME_LIMIT; @@ -352,47 +760,63 @@ primal_status_t primal_phase2(i_t phase, } } reorder_basic_list(q, basic_list); - reorder_basic_list(q, basic_list); - basis_update_t ft(L, U, p); - - std::vector c_basic(m); - for (i_t k = 0; k < m; ++k) { - const i_t j = basic_list[k]; - c_basic[k] = lp.objective[j]; - } - - // Solve B'*y = cB - ft.b_transpose_solve(c_basic, y); - settings.log.printf( - "|| y || %e || cB || %e\n", vector_norm_inf(y), vector_norm_inf(c_basic)); - - // zN = cN - N'*y - for (i_t k = 0; k < n - m; k++) { - const i_t j = nonbasic_list[k]; - // z_j <- c_j - z[j] = lp.objective[j]; + basis_update_mpf_t ft(L, U, p, settings.refactor_frequency); - // z_j <- z_j - A(:, j)'*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] * y[lp.A.i[p]]; - } - z[j] -= dot; - } - // zB = 0 - for (i_t k = 0; k < m; ++k) { - z[basic_list[k]] = 0.0; - } - settings.log.printf("|| z || %e\n", vector_norm_inf(z)); + return primal_phase2_with_advanced_basis(phase, + start_time, + lp, + settings, + vstatus, + ft, + basic_list, + nonbasic_list, + sol, + iter, + work_estimate); +} +// Note this implementation of primal simplex is experimental +// It is meant only to serve as a method to remove the perturbation to the objective +// after dual simplex has found a primal feasible solution +template +primal_status_t primal_phase2_with_advanced_basis( + i_t phase, + f_t start_time, + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + std::vector& vstatus, + basis_update_mpf_t& basis_update, + std::vector& basic_list, + std::vector& nonbasic_list, + lp_solution_t& sol, + i_t& iter, + f_t& work_estimate, + bool print_summary) +{ + const i_t m = lp.num_rows; + const i_t n = lp.num_cols; + assert(m <= n); + assert(vstatus.size() == n); + assert(lp.A.m == m); + assert(lp.A.n == n); + assert(lp.objective.size() == n); + assert(lp.lower.size() == n); + assert(lp.upper.size() == n); + assert(lp.rhs.size() == m); - set_primal_variables_on_bounds(lp, settings, z, vstatus, x); + std::vector& x = sol.x; + std::vector& y = sol.y; + std::vector& z = sol.z; - const f_t init_dual_inf = dual_infeasibility(lp, vstatus, z); - settings.log.printf("Initial dual infeasibility %e\n", init_dual_inf); + std::vector incoming_x = x; + std::vector incoming_vstatus = vstatus; + work_estimate += 2.0 * n; + settings.log.printf("Primal Simplex\n"); + // Nonbasics must be on their bounds before forming B x_B = b - A_N x_N. + // Setting them after the solve leaves ||A*x - b|| large whenever x_N != 0. + set_primal_variables_on_bounds(lp, settings, vstatus, x, work_estimate); std::vector rhs = lp.rhs; + work_estimate += m; // rhs = b - sum_{j : x_j = l_j} A(:, j) l(j) - sum_{j : x_j = u_j} A(:, j) * // u(j) for (i_t k = 0; k < n - m; ++k) { @@ -403,150 +827,556 @@ primal_status_t primal_phase2(i_t phase, for (i_t p = col_start; p < col_end; ++p) { rhs[lp.A.i[p]] -= xj * lp.A.x[p]; } + work_estimate += 3.0 * (col_end - col_start); } + work_estimate += 4 * (n - m); std::vector xB(m); - ft.b_solve(rhs, xB); + work_estimate += m; + + basis_update.b_solve(rhs, xB); for (i_t k = 0; k < m; ++k) { const i_t j = basic_list[k]; x[j] = xB[k]; } - settings.log.printf("|| x || %e\n", vector_norm2(x)); + work_estimate += 3 * m; + + constexpr bool print_norms = false; + if constexpr (print_norms) { settings.log.printf("|| x || %e\n", vector_norm2(x)); } std::vector residual = lp.rhs; + work_estimate += m; matrix_vector_multiply(lp.A, 1.0, x, -1.0, residual); + work_estimate += m + 2 * n + 4.0 * lp.A.col_start[lp.A.n]; f_t primal_residual = vector_norm_inf(residual); - if (primal_residual > 1e-6) { settings.log.printf("|| A*x - b || %e\n", primal_residual); } - f_t primal_inf = primal_infeasibility(lp, settings, vstatus, x); - settings.log.printf("Initial primal infeasibility %e\n", primal_inf); + work_estimate += m; + if (primal_residual > settings.primal_tol) { + settings.log.printf("|| A*x - b || %e\n", primal_residual); + } + + std::vector objective = lp.objective; + work_estimate += 2 * n; + const f_t primal_tol = settings.primal_tol; + f_t primal_inf = primal_infeasibility(lp, settings, vstatus, x, work_estimate); + if (primal_inf > primal_tol) { + // We are primal infeasible. Switch to phase 1 + compute_phase1_objective(lp, settings, vstatus, x, objective, work_estimate); + settings.log.printf("Phase 1\n"); + settings.log.printf("Initial primal infeasibility %e\n", primal_inf); + phase = 1; + } else { + settings.log.printf("Phase 2\n"); + phase = 2; + } + + std::vector c_basic(m); + work_estimate += m; + compute_dual_variables( + lp, settings, objective, basic_list, nonbasic_list, basis_update, c_basic, y, z, work_estimate); + if constexpr (print_norms) { settings.log.printf("|| z || %e\n", vector_norm_inf(z)); } + + i_t num_dual_inf = 0; + i_t num_primal_inf = 0; + const f_t init_dual_inf = + dual_infeasibility(lp, vstatus, z, settings.dual_tol, num_dual_inf, work_estimate); + if (num_dual_inf > 0) { settings.log.printf("Initial dual infeasibility %e\n", init_dual_inf); } - const i_t iter_limit = iter + 1000; - std::vector delta_y(m); + csr_matrix_t Arow(m, n, lp.A.nnz()); + work_estimate += n + 2 * lp.A.nnz(); + lp.A.to_compressed_row(Arow); + work_estimate += m + 6 * lp.A.nnz(); + + const i_t iter_limit = settings.iteration_limit; + const i_t start_iter = iter; + sparse_vector_t delta_y(m, 0); + sparse_vector_t etilde(m, 0); std::vector delta_z(n); std::vector delta_x(n); + work_estimate += 2 * m + 2 * n; + + f_t dual_inf = init_dual_inf; + f_t obj = compute_objective(lp, x); + work_estimate += 2 * n; + f_t pricing_dual_tol = settings.dual_tol; + primal_inf = primal_infeasibility(lp, settings, vstatus, x, num_primal_inf, work_estimate); + settings.log.printf(" Iter Objective Num Inf. Sum Inf. Time\n"); + settings.log.printf("%5d %+.16e %7d %.8e %.2f\n", + iter, + compute_user_objective(lp, obj), + phase == 1 ? num_primal_inf : num_dual_inf, + phase == 1 ? primal_inf : dual_inf, + toc(start_time)); + bool switched_phase = false; + + work_estimate += basis_update.work_estimate(); + basis_update.clear_work_estimate(); + + if (work_estimate > settings.work_limit) { + return primal_status_t::WORK_LIMIT; + } + + primal_timers_t timers(false); - settings.log.printf("Iter Objective Primal inf Dual Inf. Step Entering Leaving\n"); while (iter < iter_limit) { + timers.start_timer(work_estimate + basis_update.work_estimate()); i_t nonbasic_entering = -1; - f_t dual_inf; i_t direction; - i_t entering_index = - phase2_pricing(lp, z, nonbasic_list, vstatus, direction, nonbasic_entering, dual_inf); + i_t entering_index = phase2_pricing(lp, + z, + nonbasic_list, + vstatus, + pricing_dual_tol, + direction, + nonbasic_entering, + dual_inf, + work_estimate); + timers.pricing_time += timers.stop_timer(work_estimate + basis_update.work_estimate()); if (entering_index == -1) { - f_t obj = compute_objective(lp, x); - f_t primal_inf = primal_infeasibility(lp, settings, vstatus, x); - settings.log.printf( - "Optimal solution found. Objective %e. Dual infeas %e. Primal " - "infeasibility %e. Iterations %d\n", - compute_user_objective(lp, obj), - dual_inf, - primal_inf, - iter); - return primal_status_t::OPTIMAL; + if (phase == 2) { + // Verify optimality with a consistent basic solution: refactor, put + // nonbasics exactly on their status bounds, rebuild x_B so Ax = b, and + // refresh duals. If that point is not primal/dual feasible, continue. + if (basis_update.num_updates() > 0) { + i_t rank = basis_update.refactor_basis( + lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + if (rank == CONCURRENT_HALT_RETURN) { return primal_status_t::CONCURRENT_LIMIT; } + if (rank == TIME_LIMIT_RETURN) { return primal_status_t::TIME_LIMIT; } + if (rank != 0) { + settings.log.printf("Failed to refactor basis at optimality check. Iteration %d\n", + iter); + return primal_status_t::NUMERICAL; + } + work_estimate += basis_update.work_estimate(); + basis_update.clear_work_estimate(); + } + set_primal_variables_on_bounds(lp, settings, vstatus, x, work_estimate); + compute_basic_primal_variables( + lp, basis_update, basic_list, nonbasic_list, x, work_estimate); + compute_dual_variables(lp, + settings, + objective, + basic_list, + nonbasic_list, + basis_update, + c_basic, + y, + z, + work_estimate); + primal_inf = primal_infeasibility(lp, settings, vstatus, x, num_primal_inf, work_estimate); + dual_inf = + dual_infeasibility(lp, vstatus, z, pricing_dual_tol, num_dual_inf, work_estimate); + if (primal_inf > primal_tol) { + compute_phase1_objective(lp, settings, vstatus, x, objective, work_estimate); + phase = 1; + pricing_dual_tol = settings.dual_tol; + compute_dual_variables(lp, + settings, + objective, + basic_list, + nonbasic_list, + basis_update, + c_basic, + y, + z, + work_estimate); + settings.log.printf( + "Switching to Primal Simplex Phase 1 after near optimality. " + "Primal infeasibility %e\n", + primal_inf); + settings.log.printf(" Iter Objective Num Inf. Sum Inf. Time\n"); + switched_phase = true; + continue; + } + if (num_dual_inf > 0) { + // The refreshed reduced costs contain a candidate visible at the active + // pricing tolerance. + continue; + } + + i_t num_tight_dual_inf = 0; + const f_t tight_dual_inf = + dual_infeasibility(lp, vstatus, z, f_t(0.0), num_tight_dual_inf, work_estimate); + if (tight_dual_inf > settings.dual_tol) { + // No candidate is visible at the active pricing tolerance, but the + // zero-tolerance residual is still material. Try tighter pricing before + // accepting optimality. This is needed for problems such as cycle, + // where many small reduced-cost violations lead to improving pivots. + f_t retry_dual_tol = pricing_dual_tol; + f_t retry_dual_inf = 0.0; + i_t retry_entering = -1; + while (retry_entering == -1 && retry_dual_tol > f_t(1e-10)) { + retry_dual_tol *= f_t(0.1); + retry_entering = phase2_pricing(lp, + z, + nonbasic_list, + vstatus, + retry_dual_tol, + direction, + nonbasic_entering, + retry_dual_inf, + work_estimate); + } + if (retry_entering != -1) { + pricing_dual_tol = retry_dual_tol; + continue; + } + } + // Report the unfiltered residual at the accepted solution. + dual_inf = tight_dual_inf; + num_dual_inf = num_tight_dual_inf; + obj = compute_objective(lp, x); + work_estimate += 2 * n; + sol.objective = obj; + sol.user_objective = compute_user_objective(lp, obj); + if (!settings.inside_mip && print_summary) { + settings.log.printf("\n"); + settings.log.printf( + "Optimal solution found in %d iterations and %.2fs\n", iter, toc(start_time)); + settings.log.printf("Objective %+.8e\n", sol.user_objective); + settings.log.printf("\n"); + settings.log.printf("Primal infeasibility (abs): %.2e\n", primal_inf); + settings.log.printf("Dual infeasibility (abs): %.2e\n", dual_inf); + settings.log.printf("Primal residual ||Ax-b||: %.2e\n", + primal_constraint_residual(lp, x)); + } + timers.print_timers(settings); + return primal_status_t::OPTIMAL; + } else { + primal_inf = primal_infeasibility(lp, settings, vstatus, x, num_primal_inf, work_estimate); + + if (primal_inf > primal_tol) { + // Incremental duals may be stale relative to the current phase-I + // objective. Refresh objective and duals, then retry pricing with + // successively tighter dual tolerances. + settings.log.printf("Refreshing phase-I objective and duals. Num updates %d. Iter %d\n", + basis_update.num_updates(), + iter); + compute_phase1_objective(lp, settings, vstatus, x, objective, work_estimate); + compute_dual_variables(lp, + settings, + objective, + basic_list, + nonbasic_list, + basis_update, + c_basic, + y, + z, + work_estimate); + f_t retry_dual_tol = pricing_dual_tol; + while (entering_index == -1 && retry_dual_tol > f_t(1e-10)) { + retry_dual_tol *= f_t(0.1); + settings.log.printf("Retrying phase-I pricing with dual_tol %e\n", retry_dual_tol); + entering_index = phase2_pricing(lp, + z, + nonbasic_list, + vstatus, + retry_dual_tol, + direction, + nonbasic_entering, + dual_inf, + work_estimate); + } + if (entering_index == -1) { + settings.log.printf( + "No entering variable found with large " + "infeasibility %e (%d).\n", + primal_inf, + num_primal_inf); + return primal_status_t::PRIMAL_INFEASIBLE; + } + pricing_dual_tol = retry_dual_tol; + } else { + // Restore the objective to the original objective + objective = lp.objective; + phase = 2; + pricing_dual_tol = settings.dual_tol; + settings.log.printf( + "Primal phase I complete. Iterations %d. Time %.2f\n", iter, toc(start_time)); + settings.log.printf(" Iter Objective Num Inf. Sum Inf. Time\n"); + compute_dual_variables(lp, + settings, + objective, + basic_list, + nonbasic_list, + basis_update, + c_basic, + y, + z, + work_estimate); + obj = compute_objective(lp, x); + work_estimate += 2 * n; + dual_inf = + dual_infeasibility(lp, vstatus, z, settings.dual_tol, num_dual_inf, work_estimate); + iter++; + // Print here: continue may hit dual-optimal Phase 2 and return before + // the end-of-loop log checks switched_phase. + settings.log.printf("%5d %+.16e %7d %.8e %.2f\n", + iter, + compute_user_objective(lp, obj), + num_dual_inf, + dual_inf, + toc(start_time)); + continue; + } + } } + sparse_vector_t rhs_sparse(lp.A, entering_index); + work_estimate += 3 * rhs_sparse.i.size(); + sparse_vector_t scaled_delta_xB_sparse(m, 0); + sparse_vector_t utilde_sparse(m, 0); + timers.start_timer(work_estimate + basis_update.work_estimate()); + basis_update.b_solve(rhs_sparse, scaled_delta_xB_sparse, utilde_sparse); std::vector scaled_delta_xB(m); - std::vector rhs(m); - const i_t col_start = lp.A.col_start[entering_index]; - const i_t col_end = lp.A.col_start[entering_index + 1]; - for (i_t p = col_start; p < col_end; ++p) { - rhs[lp.A.i[p]] = lp.A.x[p]; - } - std::vector utilde(m); - ft.b_solve(rhs, scaled_delta_xB, utilde); + scaled_delta_xB_sparse.to_dense(scaled_delta_xB); + work_estimate += m + scaled_delta_xB_sparse.i.size(); for (i_t k = 0; k < m; ++k) { const i_t j = basic_list[k]; delta_x[j] = -direction * scaled_delta_xB[k]; } + work_estimate += 3 * m; for (i_t k = 0; k < n - m; ++k) { const i_t j = nonbasic_list[k]; delta_x[j] = 0.0; } + work_estimate += 2 * (n - m); delta_x[entering_index] = direction; + timers.ftran_time += timers.stop_timer(work_estimate + basis_update.work_estimate()); - std::vector residual(m); - matrix_vector_multiply(lp.A, 1.0, delta_x, 1.0, residual); +#ifdef CHECK_NULLSPACE + std::vector residual(m, 0.0); + matrix_vector_multiply(lp.A, 1.0, delta_x, 0.0, residual); f_t primal_step_err = vector_norm_inf(residual); - if (primal_step_err > 1e-3) { printf("|| A * dx || %e\n", primal_step_err); } + if (primal_step_err > 1e-3) { + settings.log.printf("|| A * dx || %e at iter %d (updates %d)\n", + primal_step_err, + iter, + basis_update.num_updates()); + } +#endif + timers.start_timer(work_estimate + basis_update.work_estimate()); i_t basic_leaving; f_t step_length; - i_t leaving_index = ratio_test(lp, vstatus, basic_list, x, delta_x, step_length, basic_leaving); - if (leaving_index == -1) { + i_t leaving_index = primal_ratio_test(lp, + settings, + vstatus, + basic_list, + x, + delta_x, + step_length, + basic_leaving, + entering_index, + direction, + work_estimate); + timers.ratio_test_time += timers.stop_timer(work_estimate + basis_update.work_estimate()); + if (leaving_index == -1 && step_length >= inf) { settings.log.printf("No leaving variable. Primal unbounded?\n"); return primal_status_t::PRIMAL_UNBOUNDED; } - assert(step_length >= 0.0); - // Update the primal variables + const bool basis_updated = (leaving_index != -1); + bool recompute_duals = false; + timers.start_timer(work_estimate + basis_update.work_estimate()); for (i_t j = 0; j < n; ++j) { x[j] += step_length * delta_x[j]; } + work_estimate += 2 * n; + timers.update_x_time += timers.stop_timer(work_estimate + basis_update.work_estimate()); + +#ifdef COMPUTE_RESIDUAL + f_t debug_primal_residual = primal_constraint_residual(lp, x); + if (debug_primal_residual > 1e-6) { + settings.log.printf("|| A * x - b || %e at iteration %d (updates %d)\n", + debug_primal_residual, + iter, + basis_update.num_updates()); + } +#endif + + if (basis_updated) { + assert(step_length >= 0.0); + + bool should_refactor = basis_update.num_updates() > settings.refactor_frequency; + f_t dual_step_length = 0.0; + if (!should_refactor) { + timers.start_timer(work_estimate + basis_update.work_estimate()); + compute_delta_y(basis_update, basic_leaving, delta_y, etilde); + timers.btran_time += timers.stop_timer(work_estimate + basis_update.work_estimate()); + const f_t pivot = scaled_delta_xB[basic_leaving]; + dual_step_length = compute_dual_step_length(z[entering_index], pivot); + } - // Update the factorization - ft.update(utilde, basic_leaving); + basic_list[basic_leaving] = entering_index; + nonbasic_list[nonbasic_entering] = leaving_index; + vstatus[entering_index] = variable_status_t::BASIC; + // Place the leaver on its leaving bound. If that bound is far from the + // current value (typical after a zero-step leave of an already-infeasible + // basic), rebuild x_B after the factor matches the new basis so Ax = b; + // phase handling below may then (re)enter Phase I if basics are infeasible. + bool rebuild_x_after_bound_snap = false; + f_t leave_bound = 0.0; + if (std::abs(lp.upper[leaving_index] - lp.lower[leaving_index]) < 1e-12) { + vstatus[leaving_index] = variable_status_t::NONBASIC_FIXED; + leave_bound = lp.lower[leaving_index]; + } else { + // Classify by which bound was hit. Using sign(delta_x) is wrong when the + // variable approached the bound from the infeasible side (phase I). + const f_t x_leave = x[leaving_index]; + const f_t dist_to_lower = std::abs(x_leave - lp.lower[leaving_index]); + const f_t dist_to_upper = std::abs(x_leave - lp.upper[leaving_index]); + if (lp.lower[leaving_index] > -inf && + (lp.upper[leaving_index] >= inf || dist_to_lower <= dist_to_upper)) { + vstatus[leaving_index] = variable_status_t::NONBASIC_LOWER; + leave_bound = lp.lower[leaving_index]; + } else { + vstatus[leaving_index] = variable_status_t::NONBASIC_UPPER; + leave_bound = lp.upper[leaving_index]; + } + } + if (std::abs(x[leaving_index] - leave_bound) > settings.primal_tol) { + rebuild_x_after_bound_snap = true; + } + x[leaving_index] = leave_bound; - // Update the basis - basic_list[basic_leaving] = entering_index; - nonbasic_list[nonbasic_entering] = leaving_index; - vstatus[entering_index] = variable_status_t::BASIC; - if (std::abs(lp.upper[leaving_index] - lp.lower[leaving_index]) < 1e-12) { - vstatus[leaving_index] = variable_status_t::NONBASIC_FIXED; - } else if (direction == 1) { - vstatus[leaving_index] = variable_status_t::NONBASIC_LOWER; + if (!should_refactor) { + timers.start_timer(work_estimate + basis_update.work_estimate()); + compute_delta_z(Arow, vstatus, delta_y, delta_z, work_estimate); + timers.delta_z_time += timers.stop_timer(work_estimate + basis_update.work_estimate()); + timers.start_timer(work_estimate + basis_update.work_estimate()); + update_y(dual_step_length, delta_y, y, work_estimate); + update_z(dual_step_length, nonbasic_list, entering_index, delta_z, z, work_estimate); + timers.update_duals_time += timers.stop_timer(work_estimate + basis_update.work_estimate()); + timers.start_timer(work_estimate + basis_update.work_estimate()); + should_refactor = basis_update.update(utilde_sparse, etilde, basic_leaving) == 1; + timers.lu_update_time += timers.stop_timer(work_estimate + basis_update.work_estimate()); + } + if (should_refactor) { + timers.start_timer(work_estimate + basis_update.work_estimate()); + i_t rank = basis_update.refactor_basis( + lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + if (rank == CONCURRENT_HALT_RETURN) { return primal_status_t::CONCURRENT_LIMIT; } + if (rank == TIME_LIMIT_RETURN) { return primal_status_t::TIME_LIMIT; } + if (rank != 0) { + settings.log.printf("Failed to refactor basis. Iteration %d\n", iter); + return primal_status_t::NUMERICAL; + } + work_estimate += basis_update.work_estimate(); + basis_update.clear_work_estimate(); + recompute_duals = true; + // Factor matches basic_list: rebuild x_B so Ax = b exactly. + set_primal_variables_on_bounds(lp, settings, vstatus, x, work_estimate); + compute_basic_primal_variables( + lp, basis_update, basic_list, nonbasic_list, x, work_estimate); + timers.lu_factorization_time += + timers.stop_timer(work_estimate + basis_update.work_estimate()); + } else if (rebuild_x_after_bound_snap) { + // FT update already matches the new basis; recompute x_B with the leaving variable + // snapped onto its bound. + compute_basic_primal_variables( + lp, basis_update, basic_list, nonbasic_list, x, work_estimate); + } } else { - vstatus[leaving_index] = variable_status_t::NONBASIC_UPPER; + if (direction > 0) { + vstatus[entering_index] = variable_status_t::NONBASIC_UPPER; + x[entering_index] = lp.upper[entering_index]; + } else { + vstatus[entering_index] = variable_status_t::NONBASIC_LOWER; + x[entering_index] = lp.lower[entering_index]; + } } - // Solve for y such that B'*y = c_B - for (i_t k = 0; k < m; ++k) { - const i_t j = basic_list[k]; - c_basic[k] = lp.objective[j]; - } - ft.b_transpose_solve(y, c_basic); - // zN = cN - N'*y - for (i_t k = 0; k < n - m; k++) { - const i_t j = nonbasic_list[k]; - // z_j <- c_j - z[j] = lp.objective[j]; - - // z_j <- z_j - A(:, j)'*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] * y[lp.A.i[p]]; + primal_inf = primal_infeasibility(lp, settings, vstatus, x, num_primal_inf, work_estimate); + if (primal_inf > primal_tol) { + if (phase != 1) { + settings.log.printf( + "Switching to Primal Simplex Phase 1. Iteration %d. Primal infeasibility %e\n", + iter, + primal_inf); + settings.log.printf(" Iter Objective Num Inf. Sum Inf. Time\n"); + switched_phase = true; } - z[j] -= dot; + compute_phase1_objective(lp, settings, vstatus, x, objective, work_estimate); + phase = 1; + recompute_duals = true; + } else if (phase == 1) { + objective = lp.objective; + phase = 2; + pricing_dual_tol = settings.dual_tol; + recompute_duals = true; + settings.log.printf( + "Primal phase I complete. Iterations %d. Time %.2f\n", iter, toc(start_time)); + settings.log.printf(" Iter Objective Num Inf. Sum Inf. Time\n"); + switched_phase = true; } - // zB = 0 - for (i_t k = 0; k < m; ++k) { - z[basic_list[k]] = 0.0; + + if (recompute_duals) { + compute_dual_variables(lp, + settings, + objective, + basic_list, + nonbasic_list, + basis_update, + c_basic, + y, + z, + work_estimate); } - const f_t obj = compute_objective(lp, x); - const f_t primal_inf = primal_infeasibility(lp, settings, vstatus, x); - settings.log.printf("%3d %.10e %.2e %.2e %.2e %d %d\n", - iter, - compute_user_objective(lp, obj), - primal_inf, - dual_inf, - step_length, - entering_index, - leaving_index); + obj = compute_objective(lp, x); + work_estimate += 2 * n; + dual_inf = dual_infeasibility(lp, vstatus, z, pricing_dual_tol, num_dual_inf, work_estimate); iter++; + + f_t now = toc(start_time); + if ((iter - start_iter) < settings.first_iteration_log || + (iter % settings.iteration_log_frequency) == 0 || switched_phase) { + const f_t user_obj = compute_user_objective(lp, obj); + settings.log.printf("%5d %+.16e %7d %.8e %.2f\n", + iter, + user_obj, + phase == 1 ? num_primal_inf : num_dual_inf, + phase == 1 ? primal_inf : dual_inf, + now); + switched_phase = false; + } + + work_estimate += basis_update.work_estimate(); + basis_update.clear_work_estimate(); + + if (now > settings.time_limit) { + timers.print_timers(settings); + return primal_status_t::TIME_LIMIT; + } + if (work_estimate > settings.work_limit) { + timers.print_timers(settings); + return primal_status_t::WORK_LIMIT; + } } - if (iter == iter_limit) { return primal_status_t::ITERATION_LIMIT; } + timers.print_timers(settings); + if (iter >= iter_limit) { return primal_status_t::ITERATION_LIMIT; } return primal_status_t::NUMERICAL; } #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE +template int primal_ratio_test(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& vstatus, + const std::vector& basic_list, + std::vector& x, + std::vector& delta_x, + double& step_length, + int& basic_leaving, + int entering_index, + int direction, + double& work_estimate); + template primal_status_t primal_phase2( int phase, double start_time, @@ -556,6 +1386,20 @@ template primal_status_t primal_phase2( lp_solution_t& sol, int& iter); +template primal_status_t primal_phase2_with_advanced_basis( + int phase, + double start_time, + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + std::vector& vstatus, + basis_update_mpf_t& basis_update, + std::vector& basic_list, + std::vector& nonbasic_list, + lp_solution_t& sol, + int& iter, + double& work_estimate, + bool print_summary); + #endif } // namespace cuopt::mathematical_optimization::simplex diff --git a/cpp/src/dual_simplex/primal.hpp b/cpp/src/dual_simplex/primal.hpp index 930958a802..fc47d90368 100644 --- a/cpp/src/dual_simplex/primal.hpp +++ b/cpp/src/dual_simplex/primal.hpp @@ -7,6 +7,7 @@ #pragma once +#include #include #include #include @@ -18,15 +19,47 @@ namespace cuopt::mathematical_optimization::simplex { enum class primal_status_t { - OPTIMAL = 0, - PRIMAL_UNBOUNDED = 1, - NUMERICAL = 2, - NOT_LOADED = 3, - TIME_LIMIT = 4, - ITERATION_LIMIT = 5, - CONCURRENT_LIMIT = 6 + OPTIMAL = 0, + PRIMAL_UNBOUNDED = 1, + PRIMAL_INFEASIBLE = 2, + NUMERICAL = 3, + TIME_LIMIT = 5, + ITERATION_LIMIT = 6, + CONCURRENT_LIMIT = 7, + WORK_LIMIT = 8, + NOT_LOADED = 9 }; +template +i_t primal_ratio_test(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& vstatus, + const std::vector& basic_list, + std::vector& x, + std::vector& delta_x, + f_t& step_length, + i_t& basic_leaving, + i_t entering_index, + i_t direction, + f_t& work_estimate); + +template +primal_status_t primal_phase2_with_advanced_basis( + i_t phase, + f_t start_time, + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + std::vector& vstatus, + basis_update_mpf_t& basis_update, + std::vector& basic_list, + std::vector& nonbasic_list, + lp_solution_t& sol, + i_t& iter, + f_t& work_estimate, + // Callers that print their own summary (dual simplex perturbation cleanup) + // suppress this one, so optimality is not reported twice. + bool print_summary = true); + template primal_status_t primal_phase2(i_t phase, f_t start_time, diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 9ec57505f9..50e8ca15a7 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -77,6 +77,7 @@ struct simplex_solver_settings_t { augmented(0), dualize(-1), ordering(-1), + initial_perturbation(-1), barrier_dual_initial_point(-1), postsolve_info(-1), qcqp_ruiz_equilibration(-1), @@ -171,6 +172,7 @@ struct simplex_solver_settings_t { i_t augmented; // -1 automatic, 0 to solve with ADAT, 1 to solve with augmented system i_t dualize; // -1 automatic, 0 to not dualize, 1 to dualize i_t ordering; // -1 automatic, 0 to use nested dissection, 1 to use AMD + i_t initial_perturbation; // -1 automatic, 0 to not perturb, 1 to perturb i_t barrier_dual_initial_point; // -1 automatic, 0 to use Lustig, Marsten, and Shanno initial // point, 1 to use initial point form dual least squares problem i_t postsolve_info; // -1 automatic (disabled), 0 disabled, 1 enabled diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 7907abd3b9..fedb9de356 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -61,6 +61,53 @@ void write_matlab(const std::string& filename, const simplex::lp_problem_t +void initialize_slack_basis_vstatus(const lp_problem_t& lp, + std::vector& vstatus) +{ + const i_t m = lp.num_rows; + const i_t n = lp.num_cols; + vstatus.resize(n); + for (i_t j = 0; j < n; ++j) { + if (lp.lower[j] == -inf && lp.upper[j] == inf) { + vstatus[j] = variable_status_t::NONBASIC_FREE; + } else if (std::abs(lp.upper[j] - lp.lower[j]) < 1e-12) { + vstatus[j] = variable_status_t::NONBASIC_FIXED; + } else if (lp.lower[j] > -inf) { + vstatus[j] = variable_status_t::NONBASIC_LOWER; + } else { + vstatus[j] = variable_status_t::NONBASIC_UPPER; + } + } + i_t num_basic = 0; + for (i_t j = n - 1; j >= 0; --j) { + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + const i_t nz = col_end - col_start; + if (nz == 1 && std::abs(lp.A.x[col_start]) == 1.0) { + vstatus[j] = variable_status_t::BASIC; + num_basic++; + } + if (num_basic == m) { break; } + } + assert(num_basic == m); +} + } // namespace template @@ -111,6 +158,7 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t& original lp_solution_t& original_solution, std::vector& vstatus, std::vector& edge_norms, + f_t& work_estimate, work_limit_context_t* work_unit_context) { raft::common::nvtx::range scope("DualSimplex::solve_lp"); @@ -129,6 +177,7 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t& original nonbasic_list, vstatus, edge_norms, + work_estimate, work_unit_context); return result; } @@ -144,6 +193,7 @@ lp_status_t solve_linear_program_with_advanced_basis( std::vector& nonbasic_list, std::vector& vstatus, std::vector& edge_norms, + f_t& work_estimate, work_limit_context_t* work_unit_context) { lp_status_t lp_status = lp_status_t::UNSET; @@ -211,6 +261,7 @@ lp_status_t solve_linear_program_with_advanced_basis( phase1_vstatus, phase1_solution, iter, + work_estimate, edge_norms, work_unit_context); } @@ -249,6 +300,7 @@ lp_status_t solve_linear_program_with_advanced_basis( nonbasic_list, solution, iter, + work_estimate, edge_norms, work_unit_context); if (status == dual_status_t::NUMERICAL) { @@ -269,6 +321,7 @@ lp_status_t solve_linear_program_with_advanced_basis( nonbasic_list, phase1_solution, iter, + work_estimate, edge_norms, work_unit_context); vstatus = phase1_vstatus; @@ -285,12 +338,19 @@ lp_status_t solve_linear_program_with_advanced_basis( nonbasic_list, solution, iter, + work_estimate, edge_norms, work_unit_context); } constexpr bool primal_cleanup = false; if (status == dual_status_t::OPTIMAL && primal_cleanup) { + settings.log.printf("Running primal cleanup\n"); primal_phase2(2, start_time, lp, settings, vstatus, solution, iter); + // TODO: We need to update ft if the basis changed + } + if (settings.inside_mip && settings.concurrent_halt != nullptr) { + settings.log.debug("Setting concurrent halt to 1 inside_mip\n"); + *settings.concurrent_halt = 1; } if (status == dual_status_t::OPTIMAL) { std::vector unscaled_x(lp.num_cols); @@ -698,6 +758,104 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us return solve_linear_program_with_barrier(user_problem, settings, start_time, solution); } +template +lp_status_t solve_linear_program_with_primal(const user_problem_t& user_problem, + const simplex_solver_settings_t& settings, + f_t start_time, + lp_solution_t& solution) +{ + raft::common::nvtx::range scope("PrimalSimplex::solve_lp"); + lp_problem_t original_lp(user_problem.handle_ptr, 1, 1, 1); + std::vector new_slacks; + dualize_info_t dualize_info; + convert_user_problem(user_problem, settings, original_lp, new_slacks, dualize_info); + + solution.resize(user_problem.num_rows, user_problem.num_cols); + lp_solution_t original_solution(original_lp.num_rows, original_lp.num_cols); + + // Presolve adds/retains artificial variables so a full slack basis exists. + lp_problem_t presolved_lp(original_lp.handle_ptr, 1, 1, 1); + presolve_info_t presolve_info; + const i_t 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; } + + lp_problem_t lp(original_lp.handle_ptr, + presolved_lp.num_rows, + presolved_lp.num_cols, + presolved_lp.A.col_start[presolved_lp.num_cols]); + std::vector column_scales; + std::vector row_scales; + scaling(presolved_lp, settings, lp, column_scales, row_scales); + + std::vector vstatus; + initialize_slack_basis_vstatus(lp, vstatus); + + lp_solution_t lp_solution(lp.num_rows, lp.num_cols); + i_t iter = 0; + const primal_status_t primal_status = + primal_phase2(2, start_time, lp, settings, vstatus, lp_solution, iter); + lp_solution.iterations = iter; + original_solution.iterations = iter; + + if (primal_status == primal_status_t::CONCURRENT_LIMIT) { + solution.iterations = iter; + return lp_status_t::CONCURRENT_LIMIT; + } + + if (primal_status == primal_status_t::OPTIMAL) { + lp_solution.objective = compute_objective(lp, lp_solution.x); + lp_solution.user_objective = compute_user_objective(lp, lp_solution.objective); + + std::vector residual = lp.rhs; + matrix_vector_multiply(lp.A, 1.0, lp_solution.x, -1.0, residual); + lp_solution.l2_primal_residual = vector_norm2(residual); + + std::vector dual_residual = lp_solution.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, lp_solution.y, 1.0, dual_residual); + lp_solution.l2_dual_residual = vector_norm2(dual_residual); + + std::vector unscaled_x(lp.num_cols); + std::vector unscaled_y(lp.num_rows); + std::vector unscaled_z(lp.num_cols); + unscale_solution(column_scales, + row_scales, + lp_solution.x, + lp_solution.y, + lp_solution.z, + unscaled_x, + unscaled_y, + unscaled_z); + uncrush_solution(presolve_info, + settings, + original_lp, + unscaled_x, + unscaled_y, + unscaled_z, + original_solution.x, + original_solution.y, + original_solution.z); + original_solution.objective = lp_solution.objective; + original_solution.user_objective = lp_solution.user_objective; + original_solution.l2_primal_residual = lp_solution.l2_primal_residual; + original_solution.l2_dual_residual = lp_solution.l2_dual_residual; + } + + uncrush_primal_solution(user_problem, original_lp, original_solution.x, solution.x); + uncrush_dual_solution( + user_problem, original_lp, original_solution.y, original_solution.z, solution.y, solution.z); + solution.objective = original_solution.objective; + solution.user_objective = original_solution.user_objective; + solution.iterations = original_solution.iterations; + solution.l2_primal_residual = original_solution.l2_primal_residual; + solution.l2_dual_residual = original_solution.l2_dual_residual; + return map_primal_status_to_lp_status(primal_status); +} + template lp_status_t solve_linear_program(const user_problem_t& user_problem, const simplex_solver_settings_t& settings, @@ -712,8 +870,9 @@ lp_status_t solve_linear_program(const user_problem_t& user_problem, lp_solution_t lp_solution(original_lp.num_rows, original_lp.num_cols); std::vector vstatus; std::vector edge_norms; + f_t work_estimate = 0.0; lp_status_t status = solve_linear_program_advanced( - original_lp, start_time, settings, lp_solution, vstatus, edge_norms); + original_lp, start_time, settings, lp_solution, vstatus, edge_norms, work_estimate); if (status == lp_status_t::CONCURRENT_LIMIT) { solution.iterations = lp_solution.iterations; return lp_status_t::CONCURRENT_LIMIT; @@ -765,8 +924,9 @@ i_t solve(const user_problem_t& problem, lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); std::vector vstatus; std::vector edge_norms; + f_t work_estimate = 0.0; lp_status_t lp_status = solve_linear_program_advanced( - original_lp, start_time, settings, solution, vstatus, edge_norms); + original_lp, start_time, settings, solution, vstatus, edge_norms, work_estimate); primal_solution = solution.x; if (lp_status == lp_status_t::OPTIMAL) { status = 0; @@ -820,6 +980,7 @@ template lp_status_t solve_linear_program_advanced( lp_solution_t& original_solution, std::vector& vstatus, std::vector& edge_norms, + double& work_estimate, work_limit_context_t* work_unit_context); template lp_status_t solve_linear_program_with_advanced_basis( @@ -832,6 +993,7 @@ template lp_status_t solve_linear_program_with_advanced_basis( std::vector& nonbasic_list, std::vector& vstatus, std::vector& edge_norms, + double& work_estimate, work_limit_context_t* work_unit_context); template lp_status_t solve_linear_program_with_barrier( @@ -845,6 +1007,12 @@ template lp_status_t solve_linear_program_with_barrier( double start_time, lp_solution_t& solution); +template lp_status_t solve_linear_program_with_primal( + const user_problem_t& user_problem, + const simplex_solver_settings_t& settings, + double start_time, + lp_solution_t& solution); + template lp_status_t solve_linear_program_with_barrier( const user_problem_t& user_problem, const simplex_solver_settings_t& settings, diff --git a/cpp/src/dual_simplex/solve.hpp b/cpp/src/dual_simplex/solve.hpp index 90c2dbd690..291675a67b 100644 --- a/cpp/src/dual_simplex/solve.hpp +++ b/cpp/src/dual_simplex/solve.hpp @@ -70,6 +70,7 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t& original lp_solution_t& original_solution, std::vector& vstatus, std::vector& edge_norms, + f_t& work_estimate, work_limit_context_t* work_unit_context = nullptr); // Solve the LP using dual simplex and keep the `basis_update_mpf_t` @@ -85,6 +86,7 @@ lp_status_t solve_linear_program_with_advanced_basis( std::vector& nonbasic_list, std::vector& vstatus, std::vector& edge_norms, + f_t& work_estimate, work_limit_context_t* work_unit_context = nullptr); template @@ -99,6 +101,11 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us lp_solution_t& solution); template +lp_status_t solve_linear_program_with_primal(const user_problem_t& user_problem, + const simplex_solver_settings_t& settings, + f_t start_time, + lp_solution_t& solution); +template lp_status_t solve_linear_program_with_barrier(const user_problem_t& user_problem, const simplex_solver_settings_t& settings, f_t start_time, diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 125280e549..ff39bbda22 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -132,12 +132,13 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_ITERATION_LIMIT, &pdlp_settings.iteration_limit, 0, std::numeric_limits::max(), std::numeric_limits::max()}, {CUOPT_NODE_LIMIT, &mip_settings.node_limit, 0, std::numeric_limits::max(), std::numeric_limits::max()}, {CUOPT_PDLP_SOLVER_MODE, reinterpret_cast(&pdlp_settings.pdlp_solver_mode), CUOPT_PDLP_SOLVER_MODE_STABLE1, CUOPT_PDLP_SOLVER_MODE_STABLE3, CUOPT_PDLP_SOLVER_MODE_STABLE3}, - {CUOPT_METHOD, reinterpret_cast(&pdlp_settings.method), CUOPT_METHOD_CONCURRENT, CUOPT_METHOD_BARRIER, CUOPT_METHOD_CONCURRENT}, + {CUOPT_METHOD, reinterpret_cast(&pdlp_settings.method), CUOPT_METHOD_CONCURRENT, CUOPT_METHOD_PRIMAL, CUOPT_METHOD_CONCURRENT}, {CUOPT_NUM_CPU_THREADS, &mip_settings.num_cpu_threads, -1, std::numeric_limits::max(), -1}, {CUOPT_AUGMENTED, &pdlp_settings.augmented, -1, 1, -1}, {CUOPT_FOLDING, &pdlp_settings.folding, -1, 1, -1}, {CUOPT_DUALIZE, &pdlp_settings.dualize, -1, 1, -1}, {CUOPT_ORDERING, &pdlp_settings.ordering, -1, 1, -1}, + {CUOPT_INITIAL_PERTURBATION, &pdlp_settings.initial_perturbation, -1, 1, -1}, {CUOPT_BARRIER_DUAL_INITIAL_POINT, &pdlp_settings.barrier_dual_initial_point, -1, 1, -1}, {CUOPT_POSTSOLVE_INFO, &pdlp_settings.postsolve_info, -1, 1, -1}, {CUOPT_MIP_CUT_PASSES, &mip_settings.max_cut_passes, -1, std::numeric_limits::max(), 10}, diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 9e54bb1a11..beb53e8a17 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -433,7 +433,7 @@ optimization_problem_solution_t convert_dual_simplex_sol( termination_status != pdlp_termination_status_t::TimeLimit && termination_status != pdlp_termination_status_t::ConcurrentLimit) { CUOPT_LOG_INFO("%s Solve status %s", - method == method_t::DualSimplex ? "Dual Simplex" : "Barrier", + method_to_string(method).c_str(), sol.get_termination_status_string().c_str()); } @@ -594,9 +594,10 @@ std::tuple, simplex::lp_status_t, f_t, f_t, f_t f_t norm_rhs = vector_norm2(user_problem.rhs); simplex::simplex_solver_settings_t dual_simplex_settings; - dual_simplex_settings.time_limit = settings.time_limit; - dual_simplex_settings.iteration_limit = settings.iteration_limit; - dual_simplex_settings.concurrent_halt = settings.concurrent_halt; + dual_simplex_settings.time_limit = settings.time_limit; + dual_simplex_settings.iteration_limit = settings.iteration_limit; + dual_simplex_settings.concurrent_halt = settings.concurrent_halt; + dual_simplex_settings.initial_perturbation = settings.initial_perturbation; if (dual_simplex_settings.concurrent_halt != nullptr) { // Don't show the dual simplex log in concurrent mode. Show the PDLP log instead dual_simplex_settings.log.log = false; @@ -639,6 +640,60 @@ optimization_problem_solution_t run_dual_simplex( method_t::DualSimplex); } +template +std::tuple, simplex::lp_status_t, f_t, f_t, f_t> run_primal( + simplex::user_problem_t& user_problem, + pdlp_solver_settings_t const& settings, + const timer_t& timer) +{ + f_t norm_user_objective = vector_norm2(user_problem.objective); + f_t norm_rhs = vector_norm2(user_problem.rhs); + + simplex::simplex_solver_settings_t primal_settings; + primal_settings.time_limit = settings.time_limit; + primal_settings.iteration_limit = settings.iteration_limit; + primal_settings.concurrent_halt = settings.concurrent_halt; + if (primal_settings.concurrent_halt != nullptr) { + // Don't show the primal simplex log in concurrent mode. Show the PDLP log instead + primal_settings.log.log = false; + } + + simplex::lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + auto status = simplex::solve_linear_program_with_primal( + user_problem, primal_settings, timer.get_tic_start(), solution); + + CUOPT_LOG_CONDITIONAL_INFO( + !settings.inside_mip, "Primal simplex finished in %.2f seconds", timer.elapsed_time()); + + if (settings.concurrent_halt != nullptr && + (status == simplex::lp_status_t::OPTIMAL || status == simplex::lp_status_t::UNBOUNDED || + status == simplex::lp_status_t::INFEASIBLE || + status == simplex::lp_status_t::UNBOUNDED_OR_INFEASIBLE)) { + // We finished. Tell PDLP to stop if it is still running. + *settings.concurrent_halt = 1; + } + + return {std::move(solution), status, timer.elapsed_time(), norm_user_objective, norm_rhs}; +} + +template +optimization_problem_solution_t run_primal( + mip::problem_t& problem, + pdlp_solver_settings_t const& settings, + const timer_t& timer) +{ + simplex::user_problem_t primal_problem = + cuopt_problem_to_user_problem(problem.handle_ptr, problem); + auto sol_primal = run_primal(primal_problem, settings, timer); + return convert_dual_simplex_sol(problem, + std::get<0>(sol_primal), + std::get<1>(sol_primal), + std::get<2>(sol_primal), + std::get<3>(sol_primal), + std::get<4>(sol_primal), + method_t::Primal); +} + #if PDLP_INSTANTIATE_FLOAT || CUOPT_INSTANTIATE_FLOAT template @@ -1803,19 +1858,27 @@ optimization_problem_solution_t solve_lp_with_method( if constexpr (std::is_same_v) { if (settings.method == method_t::DualSimplex) { return run_dual_simplex(problem, settings, timer); + } else if (settings.method == method_t::Primal) { + return run_primal(problem, settings, timer); } else if (settings.method == method_t::Barrier) { return run_barrier(problem, settings, timer); } else if (settings.method == method_t::Concurrent) { return run_concurrent(problem, settings, timer, is_batch_mode); + } else if (settings.method == method_t::PDLP) { + return run_pdlp(problem, settings, timer, is_batch_mode); } else { + cuopt_expects(false, + error_type_t::ValidationError, + "Invalid LP method. Valid values: Concurrent(0), PDLP(1), DualSimplex(2), " + "Barrier(3), Primal(4)."); return run_pdlp(problem, settings, timer, is_batch_mode); } } else { // Float precision only supports PDLP without presolve/crossover cuopt_expects(settings.method == method_t::PDLP, error_type_t::ValidationError, - "Float precision only supports PDLP method. DualSimplex, Barrier, and Concurrent " - "require double precision."); + "Float precision only supports PDLP method. DualSimplex, Primal, Barrier, and " + "Concurrent require double precision."); return run_pdlp(problem, settings, timer, is_batch_mode); } } diff --git a/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx b/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx index a5dcc78d18..73a2ccccf9 100644 --- a/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx +++ b/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx @@ -1,4 +1,4 @@ -# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # noqa +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # SPDX-License-Identifier: Apache-2.0 # cython: profile=False @@ -62,6 +62,7 @@ class SolverMethod(IntEnum): PDLP = auto() DualSimplex = auto() Barrier = auto() + Primal = auto() Unset = auto() def __str__(self):