diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index e2790858e..0e1dfc100 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -503,8 +503,6 @@ i_t convert_less_than_to_equal(const user_problem_t& user_problem, // We must convert rows in the form: a_i^T x <= beta // into: a_i^T x + s_i = beta, s_i >= 0 - csr_matrix_t Arow(0, 0, 0); - problem.A.to_compressed_row(Arow); i_t num_cols = problem.num_cols + less_rows; i_t nnz = problem.A.col_start[problem.num_cols] + less_rows; problem.A.col_start.resize(num_cols + 1); diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index dadb60cc6..d3d816a74 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -388,18 +388,62 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us settings.log.printf("Primal objective: %e\n", dot(dualize_info.primal_problem.objective, primal_solution.x)); - std::vector primal_residual = dualize_info.primal_problem.rhs; - matrix_vector_multiply( - dualize_info.primal_problem.A, 1.0, primal_solution.x, -1.0, primal_residual); std::vector inequality_rows(dualize_info.primal_problem.num_rows, 1); for (i_t i : dualize_info.equality_rows) { inequality_rows[i] = 0; } + i_t less_rows = 0; for (i_t i = 0; i < dualize_info.primal_problem.num_rows; ++i) { if (inequality_rows[i] == 1) { - primal_residual[i] = std::max(primal_residual[i], 0.0); // a_i^T x - b_i <= 0 + less_rows++; } } + // Add slack variables to the primal problem + if (less_rows > 0) + { + std::vector slack_info = dualize_info.primal_problem.rhs; + matrix_vector_multiply(dualize_info.primal_problem.A, -1.0, primal_solution.x, 1.0, slack_info); + + + lp_problem_t& problem = dualize_info.primal_problem; + i_t num_cols = problem.num_cols + less_rows; + i_t nnz = problem.A.col_start[problem.num_cols] + less_rows; + problem.A.col_start.resize(num_cols + 1); + problem.A.i.resize(nnz); + problem.A.x.resize(nnz); + problem.lower.resize(num_cols); + problem.upper.resize(num_cols); + problem.objective.resize(num_cols); + primal_solution.x.resize(num_cols); + primal_solution.z.resize(num_cols); + + i_t p = problem.A.col_start[problem.num_cols]; + i_t j = problem.num_cols; + for (i_t i = 0; i < problem.num_rows; i++) { + if (inequality_rows[i] == 1) { + problem.lower[j] = 0.0; + problem.upper[j] = INFINITY; + problem.objective[j] = 0.0; + problem.A.i[p] = i; + problem.A.x[p] = 1.0; + primal_solution.x[j] = slack_info[i]; + primal_solution.z[j] = -primal_solution.y[i]; + problem.A.col_start[j++] = p++; + inequality_rows[i] = 0; + less_rows--; + } + } + problem.A.col_start[num_cols] = p; + assert(less_rows == 0); + assert(p == nnz); + problem.A.n = num_cols; + problem.num_cols = num_cols; + } + + std::vector primal_residual = dualize_info.primal_problem.rhs; + matrix_vector_multiply( + dualize_info.primal_problem.A, 1.0, primal_solution.x, -1.0, primal_residual); + f_t primal_residual_norm = vector_norm_inf(primal_residual); const f_t norm_b = vector_norm_inf(dualize_info.primal_problem.rhs); f_t primal_relative_residual = primal_residual_norm / (1.0 + norm_b); @@ -436,6 +480,13 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us if (!settings.crossover) { return barrier_status; } if (settings.crossover && barrier_status == lp_status_t::OPTIMAL) { + + { + std::vector rhs = original_lp.rhs; + matrix_vector_multiply(original_lp.A, 1.0, lp_solution.x, -1.0, rhs); + f_t primal_residual = vector_norm_inf(rhs); + settings.log.printf("Primal residual before adding artificial variables: %e\n", primal_residual); + } // Check to see if we need to add artifical variables std::vector artificial_variables; artificial_variables.reserve(original_lp.num_rows); @@ -481,6 +532,12 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t& us lp_solution.x.size(), lp_solution.z.size()); #endif + + std::vector rhs = original_lp.rhs; + matrix_vector_multiply(original_lp.A, 1.0, lp_solution.x, -1.0, rhs); + f_t primal_residual = vector_norm_inf(rhs); + settings.log.printf("Primal residual after adding artificial variables: %e\n", primal_residual); + } // Run crossover diff --git a/cpp/src/mip/diversity/weights.cuh b/cpp/src/mip/diversity/weights.cuh index 9f53b8847..d8563ea73 100644 --- a/cpp/src/mip/diversity/weights.cuh +++ b/cpp/src/mip/diversity/weights.cuh @@ -17,6 +17,7 @@ #pragma once +#include #include #include #include diff --git a/cpp/src/mip/problem/presolve_data.cuh b/cpp/src/mip/problem/presolve_data.cuh index 20ae5fe93..e9ba5f3b3 100644 --- a/cpp/src/mip/problem/presolve_data.cuh +++ b/cpp/src/mip/problem/presolve_data.cuh @@ -20,6 +20,7 @@ #include #include +#include #include namespace cuopt { diff --git a/cpp/src/mip/relaxed_lp/lp_state.cuh b/cpp/src/mip/relaxed_lp/lp_state.cuh index 3bfa00955..8662df754 100644 --- a/cpp/src/mip/relaxed_lp/lp_state.cuh +++ b/cpp/src/mip/relaxed_lp/lp_state.cuh @@ -17,6 +17,7 @@ #pragma once +#include #include #include