From 4c3fa9838b6c5a077a2c53048d2ac41d6f69bfa0 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 22 May 2025 20:04:19 -0700 Subject: [PATCH 1/2] Fix bug: dual simplex does not fill in additional info --- cpp/src/dual_simplex/phase2.cpp | 24 +++++++++++++++ cpp/src/dual_simplex/solution.hpp | 13 ++++++-- cpp/src/dual_simplex/solve.cpp | 20 ++++++++----- cpp/src/linear_programming/solve.cu | 46 ++++++++++++++++++++--------- 4 files changed, 79 insertions(+), 24 deletions(-) diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 4d16d7c076..cbfc66a222 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -34,6 +34,28 @@ namespace cuopt::linear_programming::dual_simplex { namespace phase2 { +template +f_t l2_dual_residual(const lp_problem_t& lp, const lp_solution_t& solution) +{ + std::vector dual_residual = solution.z; + const i_t n = lp.num_cols; + // dual_residual <- z - c + for (i_t j = 0; j < n; j++) { + dual_residual[j] -= lp.objective[j]; + } + // dual_residual <- 1.0*A'*y + 1.0*(z - c) + matrix_transpose_vector_multiply(lp.A, 1.0, solution.y, 1.0, dual_residual); + return vector_norm2(dual_residual); +} + +template +f_t l2_primal_residual(const lp_problem_t& lp, const lp_solution_t& solution) +{ + std::vector primal_residual = lp.rhs; + matrix_vector_multiply(lp.A, 1.0, solution.x, -1.0, primal_residual); + return vector_norm2(primal_residual); +} + template void compute_dual_solution_from_basis(const lp_problem_t& lp, basis_update_t& ft, @@ -976,6 +998,8 @@ void prepare_optimality(const lp_problem_t& lp, } } + sol.l2_primal_residual = l2_primal_residual(lp, sol); + sol.l2_dual_residual = l2_dual_residual(lp, sol); const f_t dual_infeas = phase2::dual_infeasibility(lp, settings, vstatus, z, 0.0, 0.0); const f_t primal_infeas = phase2::primal_infeasibility(lp, settings, vstatus, x); if (phase == 1 && iter > 0) { diff --git a/cpp/src/dual_simplex/solution.hpp b/cpp/src/dual_simplex/solution.hpp index 245355fcce..e4651dd112 100644 --- a/cpp/src/dual_simplex/solution.hpp +++ b/cpp/src/dual_simplex/solution.hpp @@ -26,9 +26,16 @@ namespace cuopt::linear_programming::dual_simplex { template class lp_solution_t { public: - lp_solution_t(i_t m, i_t n) : x(n), y(m), z(n) + lp_solution_t(i_t m, i_t n) + : x(n), + y(m), + z(n), + objective(std::numeric_limits::quiet_NaN()), + user_objective(std::numeric_limits::quiet_NaN()), + iterations(0), + l2_primal_residual(std::numeric_limits::quiet_NaN()), + l2_dual_residual(std::numeric_limits::quiet_NaN()) { - objective = std::numeric_limits::quiet_NaN(); } void resize(i_t m, i_t n) @@ -47,6 +54,8 @@ class lp_solution_t { f_t objective; f_t user_objective; i_t iterations; + f_t l2_primal_residual; + f_t l2_dual_residual; }; template diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 05c1512fa9..a5b3a55911 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -204,10 +204,12 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t& original original_solution.x[j] = solution.x[j] / column_scales[j]; original_solution.z[j] = solution.z[j] / column_scales[j]; } - original_solution.y = solution.y; - original_solution.objective = solution.objective; - original_solution.user_objective = solution.user_objective; - lp_status = lp_status_t::OPTIMAL; + original_solution.y = solution.y; + original_solution.objective = solution.objective; + original_solution.user_objective = solution.user_objective; + original_solution.l2_primal_residual = solution.l2_primal_residual; + original_solution.l2_dual_residual = solution.l2_dual_residual; + lp_status = lp_status_t::OPTIMAL; } if (status == dual::status_t::DUAL_UNBOUNDED) { lp_status = lp_status_t::INFEASIBLE; } if (status == dual::status_t::TIME_LIMIT) { lp_status = lp_status_t::TIME_LIMIT; } @@ -250,10 +252,12 @@ lp_status_t solve_linear_program(const user_problem_t& user_problem, original_lp, start_time, settings, lp_solution, vstatus, edge_norms); uncrush_primal_solution(user_problem, original_lp, lp_solution.x, solution.x); uncrush_primal_solution(user_problem, original_lp, lp_solution.z, solution.z); - solution.y = lp_solution.y; - solution.objective = lp_solution.objective; - solution.user_objective = lp_solution.user_objective; - solution.iterations = lp_solution.iterations; + solution.y = lp_solution.y; + solution.objective = lp_solution.objective; + solution.user_objective = lp_solution.user_objective; + solution.iterations = lp_solution.iterations; + solution.l2_primal_residual = lp_solution.l2_primal_residual; + solution.l2_dual_residual = lp_solution.l2_dual_residual; return status; } diff --git a/cpp/src/linear_programming/solve.cu b/cpp/src/linear_programming/solve.cu index a070d12399..58624733f1 100644 --- a/cpp/src/linear_programming/solve.cu +++ b/cpp/src/linear_programming/solve.cu @@ -226,7 +226,9 @@ optimization_problem_solution_t convert_dual_simplex_sol( detail::problem_t& problem, const dual_simplex::lp_solution_t& solution, dual_simplex::lp_status_t status, - f_t duration) + f_t duration, + f_t norm_user_objective, + f_t norm_rhs) { auto to_termination_status = [](dual_simplex::lp_status_t status) { switch (status) { @@ -253,10 +255,17 @@ optimization_problem_solution_t convert_dual_simplex_sol( // Should be filled with more information from dual simplex typename optimization_problem_solution_t::additional_termination_information_t info; - info.primal_objective = solution.user_objective; - info.solve_time = duration; - info.number_of_steps_taken = solution.iterations; - info.solved_by_pdlp = false; + info.solved_by_pdlp = false; + info.primal_objective = solution.user_objective; + info.dual_objective = solution.user_objective; + info.gap = 0.0; + info.relative_gap = 0.0; + info.solve_time = duration; + info.number_of_steps_taken = solution.iterations; + info.l2_primal_residual = solution.l2_primal_residual; + info.l2_dual_residual = solution.l2_dual_residual; + info.l2_relative_primal_residual = solution.l2_primal_residual / (1.0 + norm_user_objective); + info.l2_relative_dual_residual = solution.l2_dual_residual / (1.0 + norm_rhs); pdlp_termination_status_t termination_status = to_termination_status(status); auto sol = optimization_problem_solution_t(final_primal_solution, @@ -279,12 +288,15 @@ optimization_problem_solution_t convert_dual_simplex_sol( } template -std::tuple, dual_simplex::lp_status_t, f_t> run_dual_simplex( - dual_simplex::user_problem_t& user_problem, - pdlp_solver_settings_t const& settings) +std::tuple, dual_simplex::lp_status_t, f_t, f_t, f_t> +run_dual_simplex(dual_simplex::user_problem_t& user_problem, + pdlp_solver_settings_t const& settings) { auto start_solver = std::chrono::high_resolution_clock::now(); + f_t norm_user_objective = dual_simplex::vector_norm2(user_problem.objective); + f_t norm_rhs = dual_simplex::vector_norm2(user_problem.rhs); + dual_simplex::simplex_solver_settings_t dual_simplex_settings; dual_simplex_settings.time_limit = settings.time_limit; dual_simplex_settings.iteration_limit = settings.iteration_limit; @@ -310,7 +322,7 @@ std::tuple, dual_simplex::lp_status_t, f_t settings.concurrent_halt->store(1, std::memory_order_release); } - return {std::move(solution), status, duration.count() / 1000.0}; + return {std::move(solution), status, duration.count() / 1000.0, norm_user_objective, norm_rhs}; } template @@ -324,7 +336,9 @@ optimization_problem_solution_t run_dual_simplex( return convert_dual_simplex_sol(problem, std::get<0>(sol_dual_simplex), std::get<1>(sol_dual_simplex), - std::get<2>(sol_dual_simplex)); + std::get<2>(sol_dual_simplex), + std::get<3>(sol_dual_simplex), + std::get<4>(sol_dual_simplex)); } template @@ -421,11 +435,12 @@ void run_dual_simplex_thread( dual_simplex::user_problem_t& problem, pdlp_solver_settings_t const& settings, std::unique_ptr< - std::tuple, dual_simplex::lp_status_t, f_t>>& sol_ptr) + std::tuple, dual_simplex::lp_status_t, f_t, f_t, f_t>>& + sol_ptr) { // We will return the solution from the thread as a unique_ptr sol_ptr = std::make_unique< - std::tuple, dual_simplex::lp_status_t, f_t>>( + std::tuple, dual_simplex::lp_status_t, f_t, f_t, f_t>>( run_dual_simplex(problem, settings)); } @@ -452,7 +467,8 @@ optimization_problem_solution_t run_concurrent( dual_simplex::user_problem_t dual_simplex_problem = cuopt_problem_to_simplex_problem(problem); // Create a thread for dual simplex - std::unique_ptr, dual_simplex::lp_status_t, f_t>> + std::unique_ptr< + std::tuple, dual_simplex::lp_status_t, f_t, f_t, f_t>> sol_dual_simplex_ptr; std::thread dual_simplex_thread(run_dual_simplex_thread, std::ref(dual_simplex_problem), @@ -469,7 +485,9 @@ optimization_problem_solution_t run_concurrent( auto sol_dual_simplex = convert_dual_simplex_sol(problem, std::get<0>(*sol_dual_simplex_ptr), std::get<1>(*sol_dual_simplex_ptr), - std::get<2>(*sol_dual_simplex_ptr)); + std::get<2>(*sol_dual_simplex_ptr), + std::get<3>(*sol_dual_simplex_ptr), + std::get<4>(*sol_dual_simplex_ptr)); f_t end_time = dual_simplex::toc(start_time); CUOPT_LOG_INFO("Concurrent time: %.3fs", end_time); From ee8b4c9c0e7f052664ece5dcda7e9af57c9a44b7 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Fri, 23 May 2025 08:33:03 -0700 Subject: [PATCH 2/2] Add more missing fields --- cpp/src/linear_programming/solve.cu | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/cpp/src/linear_programming/solve.cu b/cpp/src/linear_programming/solve.cu index 58624733f1..b4478fbba5 100644 --- a/cpp/src/linear_programming/solve.cu +++ b/cpp/src/linear_programming/solve.cu @@ -255,17 +255,22 @@ optimization_problem_solution_t convert_dual_simplex_sol( // Should be filled with more information from dual simplex typename optimization_problem_solution_t::additional_termination_information_t info; - info.solved_by_pdlp = false; - info.primal_objective = solution.user_objective; - info.dual_objective = solution.user_objective; - info.gap = 0.0; - info.relative_gap = 0.0; - info.solve_time = duration; - info.number_of_steps_taken = solution.iterations; - info.l2_primal_residual = solution.l2_primal_residual; - info.l2_dual_residual = solution.l2_dual_residual; - info.l2_relative_primal_residual = solution.l2_primal_residual / (1.0 + norm_user_objective); - info.l2_relative_dual_residual = solution.l2_dual_residual / (1.0 + norm_rhs); + info.solved_by_pdlp = false; + info.primal_objective = solution.user_objective; + info.dual_objective = solution.user_objective; + info.gap = 0.0; + info.relative_gap = 0.0; + info.solve_time = duration; + info.number_of_steps_taken = solution.iterations; + info.total_number_of_attempted_steps = solution.iterations; + info.l2_primal_residual = solution.l2_primal_residual; + info.l2_dual_residual = solution.l2_dual_residual; + info.l2_relative_primal_residual = solution.l2_primal_residual / (1.0 + norm_user_objective); + info.l2_relative_dual_residual = solution.l2_dual_residual / (1.0 + norm_rhs); + info.max_primal_ray_infeasibility = 0.0; + info.primal_ray_linear_objective = 0.0; + info.max_dual_ray_infeasibility = 0.0; + info.dual_ray_linear_objective = 0.0; pdlp_termination_status_t termination_status = to_termination_status(status); auto sol = optimization_problem_solution_t(final_primal_solution,