Skip to content

Fix bug: dual simplex does not fill in additional info #33

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rapids-bot merged 4 commits into from
May 23, 2025
Merged

Fix bug: dual simplex does not fill in additional info #33

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
4c3fa98
Fix bug: dual simplex does not fill in additional info
chris-maes May 23, 2025
ee8b4c9
Add more missing fields
chris-maes May 23, 2025
90bb612
Merge branch 'branch-25.05-nvidia' into dual_simplex_term_info
chris-maes May 23, 2025
80350de
Merge branch 'branch-25.05-nvidia' into dual_simplex_term_info
chris-maes May 23, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 24 additions & 0 deletions cpp/src/dual_simplex/phase2.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,28 @@ namespace cuopt::linear_programming::dual_simplex {

namespace phase2 {

template <typename i_t, typename f_t>
f_t l2_dual_residual(const lp_problem_t<i_t, f_t>& lp, const lp_solution_t<i_t, f_t>& solution)
{
std::vector<f_t> 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<i_t, f_t>(dual_residual);
}

template <typename i_t, typename f_t>
f_t l2_primal_residual(const lp_problem_t<i_t, f_t>& lp, const lp_solution_t<i_t, f_t>& solution)
{
std::vector<f_t> primal_residual = lp.rhs;
matrix_vector_multiply(lp.A, 1.0, solution.x, -1.0, primal_residual);
return vector_norm2<i_t, f_t>(primal_residual);
}

template <typename i_t, typename f_t>
void compute_dual_solution_from_basis(const lp_problem_t<i_t, f_t>& lp,
basis_update_t<i_t, f_t>& ft,
Expand Down Expand Up @@ -976,6 +998,8 @@ void prepare_optimality(const lp_problem_t<i_t, f_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) {
Expand Down
13 changes: 11 additions & 2 deletions cpp/src/dual_simplex/solution.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,16 @@ namespace cuopt::linear_programming::dual_simplex {
template <typename i_t, typename f_t>
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<f_t>::quiet_NaN()),
user_objective(std::numeric_limits<f_t>::quiet_NaN()),
iterations(0),
l2_primal_residual(std::numeric_limits<f_t>::quiet_NaN()),
l2_dual_residual(std::numeric_limits<f_t>::quiet_NaN())
{
objective = std::numeric_limits<f_t>::quiet_NaN();
}

void resize(i_t m, i_t n)
Expand All @@ -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 <typename i_t, typename f_t>
Expand Down
20 changes: 12 additions & 8 deletions cpp/src/dual_simplex/solve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,10 +205,12 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original
unscale_solution<i_t, f_t>(column_scales, solution.x, solution.z, unscaled_x, unscaled_z);
uncrush_solution(
presolve_info, unscaled_x, unscaled_z, original_solution.x, original_solution.z);
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; }
Expand Down Expand Up @@ -251,10 +253,12 @@ lp_status_t solve_linear_program(const user_problem_t<i_t, f_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;
}

Expand Down
51 changes: 37 additions & 14 deletions cpp/src/linear_programming/solve.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -226,7 +226,9 @@ optimization_problem_solution_t<i_t, f_t> convert_dual_simplex_sol(
detail::problem_t<i_t, f_t>& problem,
const dual_simplex::lp_solution_t<i_t, f_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) {
Expand All @@ -253,10 +255,22 @@ optimization_problem_solution_t<i_t, f_t> convert_dual_simplex_sol(

// Should be filled with more information from dual simplex
typename optimization_problem_solution_t<i_t, f_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.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<i_t, f_t>(final_primal_solution,
Expand All @@ -279,12 +293,15 @@ optimization_problem_solution_t<i_t, f_t> convert_dual_simplex_sol(
}

template <typename i_t, typename f_t>
std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, dual_simplex::lp_status_t, f_t> run_dual_simplex(
dual_simplex::user_problem_t<i_t, f_t>& user_problem,
pdlp_solver_settings_t<i_t, f_t> const& settings)
std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, dual_simplex::lp_status_t, f_t, f_t, f_t>
run_dual_simplex(dual_simplex::user_problem_t<i_t, f_t>& user_problem,
pdlp_solver_settings_t<i_t, f_t> const& settings)
{
auto start_solver = std::chrono::high_resolution_clock::now();

f_t norm_user_objective = dual_simplex::vector_norm2<i_t, f_t>(user_problem.objective);
f_t norm_rhs = dual_simplex::vector_norm2<i_t, f_t>(user_problem.rhs);

dual_simplex::simplex_solver_settings_t<i_t, f_t> dual_simplex_settings;
dual_simplex_settings.time_limit = settings.time_limit;
dual_simplex_settings.iteration_limit = settings.iteration_limit;
Expand All @@ -310,7 +327,7 @@ std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, 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 <typename i_t, typename f_t>
Expand All @@ -324,7 +341,9 @@ optimization_problem_solution_t<i_t, f_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 <typename i_t, typename f_t>
Expand Down Expand Up @@ -421,11 +440,12 @@ void run_dual_simplex_thread(
dual_simplex::user_problem_t<i_t, f_t>& problem,
pdlp_solver_settings_t<i_t, f_t> const& settings,
std::unique_ptr<
std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, dual_simplex::lp_status_t, f_t>>& sol_ptr)
std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, 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_solution_t<i_t, f_t>, dual_simplex::lp_status_t, f_t>>(
std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, dual_simplex::lp_status_t, f_t, f_t, f_t>>(
run_dual_simplex(problem, settings));
}

Expand All @@ -452,7 +472,8 @@ optimization_problem_solution_t<i_t, f_t> run_concurrent(
dual_simplex::user_problem_t<i_t, f_t> dual_simplex_problem =
cuopt_problem_to_simplex_problem<i_t, f_t>(problem);
// Create a thread for dual simplex
std::unique_ptr<std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, dual_simplex::lp_status_t, f_t>>
std::unique_ptr<
std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, dual_simplex::lp_status_t, f_t, f_t, f_t>>
sol_dual_simplex_ptr;
std::thread dual_simplex_thread(run_dual_simplex_thread<i_t, f_t>,
std::ref(dual_simplex_problem),
Expand All @@ -469,7 +490,9 @@ optimization_problem_solution_t<i_t, f_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);
Expand Down