diff --git a/cpp/cuopt_cli.cpp b/cpp/cuopt_cli.cpp index 0cac35262f..ac568e07cf 100644 --- a/cpp/cuopt_cli.cpp +++ b/cpp/cuopt_cli.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include @@ -92,11 +93,15 @@ int run_single_file(const std::string& file_path, bool solve_relaxation, cuopt::linear_programming::solver_settings_t& settings) { + cuopt::init_logger_t log(settings.get_parameter(CUOPT_LOG_FILE), + settings.get_parameter(CUOPT_LOG_TO_CONSOLE)); + std::string base_filename = file_path.substr(file_path.find_last_of("/\\") + 1); constexpr bool input_mps_strict = false; cuopt::mps_parser::mps_data_model_t mps_data_model; bool parsing_failed = false; + auto timer = cuopt::timer_t(settings.get_parameter(CUOPT_TIME_LIMIT)); { CUOPT_LOG_INFO("Reading file %s", base_filename.c_str()); try { @@ -111,6 +116,7 @@ int run_single_file(const std::string& file_path, CUOPT_LOG_ERROR("Parsing MPS failed. Exiting!"); return -1; } + CUOPT_LOG_INFO("Read file %s in %.2f seconds", base_filename.c_str(), timer.elapsed_time()); // Determine memory backend and create problem using interface // Create handle only for GPU memory backend (avoid CUDA init on CPU-only hosts) diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index 075323744d..4da66abe77 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -2151,10 +2151,9 @@ void barrier_solver_t::gpu_compute_residual_norms(const rmm::device_uv std::max(device_vector_norm_inf(data.d_primal_residual_, stream_view_), device_vector_norm_inf(data.d_bound_residual_, stream_view_)); dual_residual_norm = device_vector_norm_inf(data.d_dual_residual_, stream_view_); - // TODO: CMM understand why rhs and not residual complementarity_residual_norm = - std::max(device_vector_norm_inf(data.d_complementarity_xz_rhs_, stream_view_), - device_vector_norm_inf(data.d_complementarity_wv_rhs_, stream_view_)); + std::max(device_vector_norm_inf(data.d_complementarity_xz_residual_, stream_view_), + device_vector_norm_inf(data.d_complementarity_wv_residual_, stream_view_)); } template @@ -3492,7 +3491,9 @@ lp_status_t barrier_solver_t::solve(f_t start_time, f_t relative_primal_residual = primal_residual_norm / (1.0 + norm_b); f_t relative_dual_residual = dual_residual_norm / (1.0 + norm_c); f_t relative_complementarity_residual = - complementarity_residual_norm / (1.0 + std::abs(primal_objective)); + complementarity_residual_norm / + (1.0 + std::min(std::abs(compute_user_objective(lp, primal_objective)), + std::abs(primal_objective))); dense_vector_t upper(lp.upper); data.gather_upper_bounds(upper, data.restrict_u_); @@ -3508,11 +3509,11 @@ lp_status_t barrier_solver_t::solve(f_t start_time, float64_t elapsed_time = toc(start_time); settings.log.printf("%3d %+.12e %+.12e %.2e %.2e %.2e %.1f\n", iter, - primal_objective, - dual_objective, - primal_residual_norm, - dual_residual_norm, - complementarity_residual_norm, + compute_user_objective(lp, primal_objective), + compute_user_objective(lp, dual_objective), + relative_primal_residual, + relative_dual_residual, + relative_complementarity_residual, elapsed_time); bool converged = primal_residual_norm < settings.barrier_relative_feasibility_tol && @@ -3654,7 +3655,9 @@ lp_status_t barrier_solver_t::solve(f_t start_time, relative_primal_residual = primal_residual_norm / (1.0 + norm_b); relative_dual_residual = dual_residual_norm / (1.0 + norm_c); relative_complementarity_residual = - complementarity_residual_norm / (1.0 + std::abs(primal_objective)); + complementarity_residual_norm / + (1.0 + std::min(std::abs(compute_user_objective(lp, primal_objective)), + std::abs(primal_objective))); if (relative_primal_residual < settings.barrier_relaxed_feasibility_tol && relative_dual_residual < settings.barrier_relaxed_optimality_tol && diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index d2a68d96de..c5ef847106 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -821,6 +821,168 @@ i_t presolve(const lp_problem_t& original, { problem = original; std::vector row_sense(problem.num_rows, '='); + // Check for free variables + i_t free_variables = 0; + for (i_t j = 0; j < problem.num_cols; j++) { + if (problem.lower[j] == -inf && problem.upper[j] == inf) { free_variables++; } + } + + if (settings.barrier_presolve && free_variables > 0) { + // Try to remove free variables + std::vector constraints_to_check; + std::vector current_free_variables; + std::vector row_marked(problem.num_rows, 0); + current_free_variables.reserve(problem.num_cols); + constraints_to_check.reserve(problem.num_rows); + for (i_t j = 0; j < problem.num_cols; j++) { + if (problem.lower[j] == -inf && problem.upper[j] == inf) { + current_free_variables.push_back(j); + const i_t col_start = problem.A.col_start[j]; + const i_t col_end = problem.A.col_start[j + 1]; + for (i_t p = col_start; p < col_end; p++) { + const i_t i = problem.A.i[p]; + if (row_marked[i] == 0) { + row_marked[i] = 1; + constraints_to_check.push_back(i); + } + } + } + } + + i_t removed_free_variables = 0; + + if (constraints_to_check.size() > 0) { + // Check if the constraints are feasible + csr_matrix_t Arow(0, 0, 0); + problem.A.to_compressed_row(Arow); + + // The constraints are in the form: + // sum_j a_j x_j = beta + for (i_t i : constraints_to_check) { + const i_t row_start = Arow.row_start[i]; + const i_t row_end = Arow.row_start[i + 1]; + f_t lower_activity_i = 0.0; + f_t upper_activity_i = 0.0; + i_t lower_inf_i = 0; + i_t upper_inf_i = 0; + i_t last_free_i = -1; + f_t last_free_coeff_i = 0.0; + for (i_t p = row_start; p < row_end; p++) { + const i_t j = Arow.j[p]; + const f_t aij = Arow.x[p]; + const f_t lower_j = problem.lower[j]; + const f_t upper_j = problem.upper[j]; + if (lower_j == -inf && upper_j == inf) { + last_free_i = j; + last_free_coeff_i = aij; + } + if (aij > 0) { + if (lower_j > -inf) { + lower_activity_i += aij * lower_j; + } else { + lower_inf_i++; + } + if (upper_j < inf) { + upper_activity_i += aij * upper_j; + } else { + upper_inf_i++; + } + } else { + if (upper_j < inf) { + lower_activity_i += aij * upper_j; + } else { + lower_inf_i++; + } + if (lower_j > -inf) { + upper_activity_i += aij * lower_j; + } else { + upper_inf_i++; + } + } + } + + if (last_free_i == -1) { continue; } + + // sum_j a_ij x_j == beta + + const f_t rhs = problem.rhs[i]; + // sum_{k != j} a_ik x_k + a_ij x_j == rhs + // Suppose that -inf < x_j < inf and all other variables x_k with k != j are bounded + // a_ij x_j == rhs - sum_{k != j} a_ik x_k + // So if a_ij > 0, we have + // x_j == 1/a_ij * (rhs - sum_{k != j} a_ik x_k) + // We can derive two bounds from this: + // x_j <= 1/a_ij * (rhs - lower_activity_i) and + // x_j >= 1/a_ij * (rhs - upper_activity_i) + + // If a_ij < 0, we have + // x_j == 1/a_ij * (rhs - sum_{k != j} a_ik x_k + // And we can derive two bounds from this: + // x_j >= 1/a_ij * (rhs - lower_activity_i) + // x_j <= 1/a_ij * (rhs - upper_activity_i) + const i_t j = last_free_i; + const f_t a_ij = last_free_coeff_i; + const f_t max_bound = 1e10; + bool bounded = false; + if (a_ij > 0) { + if (lower_inf_i == 1) { + const f_t new_upper = 1.0 / a_ij * (rhs - lower_activity_i); + if (new_upper < max_bound) { + problem.upper[j] = new_upper; + bounded = true; + } + } + if (upper_inf_i == 1) { + const f_t new_lower = 1.0 / a_ij * (rhs - upper_activity_i); + if (new_lower > -max_bound) { + problem.lower[j] = new_lower; + bounded = true; + } + } + } else if (a_ij < 0) { + if (lower_inf_i == 1) { + const f_t new_lower = 1.0 / a_ij * (rhs - lower_activity_i); + if (new_lower > -max_bound) { + problem.lower[j] = new_lower; + bounded = true; + } + } + if (upper_inf_i == 1) { + const f_t new_upper = 1.0 / a_ij * (rhs - upper_activity_i); + if (new_upper < max_bound) { + problem.upper[j] = new_upper; + bounded = true; + } + } + } + + if (bounded) { removed_free_variables++; } + } + } + + for (i_t j : current_free_variables) { + if (problem.lower[j] > -inf && problem.upper[j] < inf) { + // We don't need two bounds. Pick the smallest one. + if (std::abs(problem.lower[j]) < std::abs(problem.upper[j])) { + // Restore the inf in the upper bound. Barrier will not require an additional w variable + problem.upper[j] = inf; + } else { + // Restores the -inf in the lower bound. Barrier will require an additional w variable + problem.lower[j] = -inf; + } + } + } + + i_t new_free_variables = 0; + for (i_t j = 0; j < problem.num_cols; j++) { + if (problem.lower[j] == -inf && problem.upper[j] == inf) { new_free_variables++; } + } + if (removed_free_variables != 0) { + settings.log.printf("Bounded %d free variables\n", removed_free_variables); + } + assert(new_free_variables == free_variables - removed_free_variables); + free_variables = new_free_variables; + } // The original problem may have a variable without a lower bound // but a finite upper bound @@ -834,7 +996,43 @@ i_t presolve(const lp_problem_t& original, settings.log.printf("%d variables with no lower bound\n", no_lower_bound); } - // FIXME:: handle no lower bound case for barrier presolve + // Handle -inf < x_j <= u_j by substituting x'_j = -x_j, giving -u_j <= x'_j < inf + if (settings.barrier_presolve && no_lower_bound > 0) { + presolve_info.negated_variables.reserve(no_lower_bound); + for (i_t j = 0; j < problem.num_cols; j++) { + if (problem.lower[j] == -inf && problem.upper[j] < inf) { + presolve_info.negated_variables.push_back(j); + + problem.lower[j] = -problem.upper[j]; + problem.upper[j] = inf; + problem.objective[j] *= -1; + + const i_t col_start = problem.A.col_start[j]; + const i_t col_end = problem.A.col_start[j + 1]; + for (i_t p = col_start; p < col_end; p++) { + problem.A.x[p] *= -1.0; + } + } + } + + // (1/2) x^T Q x with x = D x' (D_ii = -1 for negated columns) is (1/2) x'^T D Q D x'. + // One pass: Q'_{ik} = D_{ii} D_{kk} Q_{ik} — flip iff exactly one of {i,k} is negated. + if (problem.Q.n > 0 && !presolve_info.negated_variables.empty()) { + std::vector is_negated(static_cast(problem.num_cols), false); + for (i_t const j : presolve_info.negated_variables) { + is_negated[static_cast(j)] = true; + } + for (i_t row = 0; row < problem.Q.m; ++row) { + const i_t q_start = problem.Q.row_start[row]; + const i_t q_end = problem.Q.row_start[row + 1]; + const bool is_negated_row = is_negated[static_cast(row)]; + for (i_t p = q_start; p < q_end; ++p) { + const i_t col = problem.Q.j[p]; + if (is_negated_row != is_negated[static_cast(col)]) { problem.Q.x[p] *= -1.0; } + } + } + } + } // The original problem may have nonzero lower bounds // 0 != l_j <= x_j <= u_j @@ -939,12 +1137,6 @@ i_t presolve(const lp_problem_t& original, remove_empty_cols(problem, num_empty_cols, presolve_info); } - // Check for free variables - i_t free_variables = 0; - for (i_t j = 0; j < problem.num_cols; j++) { - if (problem.lower[j] == -inf && problem.upper[j] == inf) { free_variables++; } - } - problem.Q.check_matrix("Before free variable expansion"); if (settings.barrier_presolve && free_variables > 0) { @@ -1511,6 +1703,14 @@ void uncrush_solution(const presolve_info_t& presolve_info, } settings.log.printf("Post-solve: Handling removed lower bounds %d\n", num_lower_bounds); } + + if (presolve_info.negated_variables.size() > 0) { + for (const i_t j : presolve_info.negated_variables) { + input_x[j] *= -1.0; + input_z[j] *= -1.0; + } + } + assert(uncrushed_x.size() == input_x.size()); assert(uncrushed_y.size() == input_y.size()); assert(uncrushed_z.size() == input_z.size()); diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index a068ed04ab..d570ea933e 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -202,6 +202,9 @@ struct presolve_info_t { std::vector removed_constraints; folding_info_t folding_info; + + // Variables that were negated to handle -inf < x_j <= u_j + std::vector negated_variables; }; template