diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index 075323744d..14297f67df 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -1971,6 +1971,14 @@ int barrier_solver_t::initial_point(iteration_data_t& data) settings.log.printf("min v %e min z %e\n", data.v.minimum(), data.z.minimum()); #endif + if (x_pdhg.size() > 0) { + data.x = x_pdhg; + data.w = w_pdhg; + data.y = y_pdhg; + data.z = z_pdhg; + data.v = v_pdhg; + } + return 0; } diff --git a/cpp/src/barrier/barrier.hpp b/cpp/src/barrier/barrier.hpp index 014a46db35..0bdcb3e198 100644 --- a/cpp/src/barrier/barrier.hpp +++ b/cpp/src/barrier/barrier.hpp @@ -76,6 +76,12 @@ class barrier_solver_t { // To be able to directly pass lambdas to transform functions public: + + std::vector x_pdhg; + std::vector w_pdhg; + std::vector y_pdhg; + std::vector z_pdhg; + std::vector v_pdhg; void compute_next_iterate(iteration_data_t& data, f_t step_scale, f_t step_primal, diff --git a/cpp/src/dual_simplex/CMakeLists.txt b/cpp/src/dual_simplex/CMakeLists.txt index 7341c0de30..502bfc8577 100644 --- a/cpp/src/dual_simplex/CMakeLists.txt +++ b/cpp/src/dual_simplex/CMakeLists.txt @@ -7,7 +7,7 @@ set(DUAL_SIMPLEX_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/basis_solves.cpp ${CMAKE_CURRENT_SOURCE_DIR}/basis_updates.cpp ${CMAKE_CURRENT_SOURCE_DIR}/bound_flipping_ratio_test.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/crossover.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/crossover.cu ${CMAKE_CURRENT_SOURCE_DIR}/folding.cpp ${CMAKE_CURRENT_SOURCE_DIR}/initial_basis.cpp ${CMAKE_CURRENT_SOURCE_DIR}/phase1.cpp diff --git a/cpp/src/dual_simplex/crossover.cpp b/cpp/src/dual_simplex/crossover.cu similarity index 78% rename from cpp/src/dual_simplex/crossover.cpp rename to cpp/src/dual_simplex/crossover.cu index 988c9c50ad..fd7240819f 100644 --- a/cpp/src/dual_simplex/crossover.cpp +++ b/cpp/src/dual_simplex/crossover.cu @@ -11,10 +11,16 @@ #include #include #include +#include #include #include #include + +#include +#include +#include + #include #include @@ -411,7 +417,7 @@ i_t dual_push(const lp_problem_t& lp, // U^T*etilde = es. // We have that B^T = U^T*L^T so B^T delta_y = U^T*etilde = -delta_zs*es // So we need to divide by -delta_zs - for (i_t k = 0; k < UTsol_sparse.i.size(); ++k) { + for (size_t k = 0; k < UTsol_sparse.i.size(); ++k) { UTsol_sparse.x[k] /= -delta_zs; } @@ -675,7 +681,7 @@ f_t primal_ratio_test(const lp_problem_t& lp, const i_t n = lp.num_cols; f_t step_length = 1.0; constexpr f_t pivot_tol = 1e-9; - for (i_t k = 0; k < delta_xB.i.size(); ++k) { + for (size_t k = 0; k < delta_xB.i.size(); ++k) { const i_t j = basic_list[delta_xB.i[k]]; if (x[j] <= lp.upper[j] && delta_xB.x[k] > pivot_tol && lp.upper[j] < inf) { const f_t ratio = (lp.upper[j] - x[j]) / delta_xB.x[k]; @@ -782,7 +788,7 @@ i_t primal_push(const lp_problem_t& lp, // Compute delta_xB_trial = -delta_xs * w sparse_vector_t& delta_xB_trial = delta_xB_trials[push]; delta_xB_trial = w_sparse; - for (i_t k = 0; k < w_sparse.i.size(); ++k) { + for (size_t k = 0; k < w_sparse.i.size(); ++k) { delta_xB_trial.x[k] = -w_sparse.x[k] * delta_xs; } leaving_index_trials[push] = -1; @@ -823,7 +829,7 @@ i_t primal_push(const lp_problem_t& lp, #endif // xB <- xB + step_length * delta_xB - for (i_t k = 0; k < delta_xB.i.size(); ++k) { + for (size_t k = 0; k < delta_xB.i.size(); ++k) { const i_t j = basic_list[delta_xB.i[k]]; x[j] += step_length * delta_xB.x[k]; } @@ -1016,7 +1022,7 @@ i_t find_candidate_columns(const lp_problem_t& lp, const i_t n = lp.num_cols; const i_t m = lp.num_rows; f_t basis_threshold = 1e-10; - while (candidate_columns.size() < m && basis_threshold < 1.0) { + while (static_cast(candidate_columns.size()) < m && basis_threshold < 1.0) { basis_threshold *= 10.0; for (i_t j = 0; j < n; ++j) { const f_t lower_bound_slack = solution.x[j] - lp.lower[j]; @@ -1217,7 +1223,7 @@ crossover_status_t crossover(const lp_problem_t& lp, settings.log.debug("num basic %d from slacks\n", num_basic); } - for (i_t k = 0; k < candidate_columns.size(); k++) { + for (size_t k = 0; k < candidate_columns.size(); k++) { const i_t j = candidate_columns[k]; vstatus[j] = vstatus_for_candidates[k]; if (vstatus[j] == variable_status_t::BASIC) { num_basic++; } @@ -1572,6 +1578,444 @@ crossover_status_t crossover(const lp_problem_t& lp, return status; } +template +crossover_status_t central_path(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const lp_solution_t& initial_solution, + f_t start_time, + lp_solution_t& solution, + std::vector& vstatus) +{ + + printf("\n\n\nStarting central path solver\n"); + // First apply barrier presolve to the problem + simplex_solver_settings_t barrier_settings = settings; + barrier_settings.barrier_presolve = true; + barrier_settings.folding = false; + barrier_settings.barrier = true; + + std::vector original_dual_residual = initial_solution.z; + for (i_t j = 0; j < lp.num_cols; j++) { + original_dual_residual[j] -= lp.objective[j]; + } + matrix_transpose_vector_multiply(lp.A, 1.0, initial_solution.y, +1.0, original_dual_residual); + // printf("Original dual residual %e\n", vector_norm_inf(original_dual_residual)); + + lp_problem_t presolved_lp(lp.handle_ptr, 1, 1, 1); + presolve_info_t presolve_info; + presolve(lp, barrier_settings, presolved_lp, presolve_info); + + std::vector x; + std::vector y; + std::vector z; + crush_solution_to_presolve_space( + lp, presolve_info, initial_solution.x, initial_solution.y, initial_solution.z, x, y, z); + + // Compute primal infeasibility + std::vector residual = presolved_lp.rhs; + matrix_vector_multiply(presolved_lp.A, 1.0, x, -1.0, residual); + + // Compute dual infeasibility + std::vector dual_residual_presolved = z; + for (size_t j = 0; j < z.size(); j++) { + dual_residual_presolved[j] -= presolved_lp.objective[j]; + } + matrix_transpose_vector_multiply(presolved_lp.A, 1.0, y, +1.0, dual_residual_presolved); + + f_t primal_infeas = vector_norm_inf(residual); + f_t dual_infeas = vector_norm_inf(dual_residual_presolved); + + printf("Primal residual %e Dual residual %e\n", primal_infeas, dual_infeas); + + f_t primal_inf = 0.0; + i_t n = presolved_lp.num_cols; + for (i_t j = 0; j < n; ++j) { + if (x[j] < presolved_lp.lower[j]) { primal_inf += presolved_lp.lower[j] - x[j]; } + if (x[j] > presolved_lp.upper[j]) { primal_inf += x[j] - presolved_lp.upper[j]; } + } + + // printf("Primal infeasibility %e\n", primal_inf); + + // Currently we have + // minimize c^T x + // subject to A x = b + // 0 <= x <= u + + // We translate this into + // minimize c^T x + // subject to A x = b + // x + w = u + // x, w >= 0 + + lp_problem_t barrier_lp = presolved_lp; + + i_t n_upper_bounds = 0; + for (i_t j = 0; j < n; ++j) { + if (presolved_lp.upper[j] < inf) { n_upper_bounds++; } + } + + // Create the upper bounds vector + std::vector upper_bounds(n_upper_bounds); + n_upper_bounds = 0; + for (i_t j = 0; j < n; j++) { + if (presolved_lp.upper[j] < inf) { upper_bounds[n_upper_bounds++] = j; } + } + if (n_upper_bounds > 0) { + // printf("Upper bounds : %d\n", n_upper_bounds); + + csr_matrix_t Arow(1, 1, 1); + presolved_lp.A.to_compressed_row(Arow); + + Arow.row_start.resize(presolved_lp.num_rows + n_upper_bounds + 1); + i_t nz = presolved_lp.A.col_start[n]; + Arow.j.resize(nz + 2 * n_upper_bounds); + Arow.x.resize(nz + 2 * n_upper_bounds); + + for (i_t i = presolved_lp.num_rows; i < presolved_lp.num_rows + n_upper_bounds; i++) { + Arow.row_start[i] = nz; + const i_t k = i - presolved_lp.num_rows; + Arow.j[nz] = upper_bounds[k]; + Arow.x[nz] = 1.0; + nz++; + Arow.j[nz] = n + k; + Arow.x[nz] = 1.0; + nz++; + } + Arow.row_start[presolved_lp.num_rows + n_upper_bounds] = nz; + Arow.m = presolved_lp.num_rows + n_upper_bounds; + Arow.n = n + n_upper_bounds; + Arow.nz_max = nz; + + Arow.to_compressed_col(barrier_lp.A); + + barrier_lp.rhs.resize(presolved_lp.num_rows + n_upper_bounds); + for (i_t i = presolved_lp.num_rows; i < presolved_lp.num_rows + n_upper_bounds; i++) { + const i_t k = i - presolved_lp.num_rows; + barrier_lp.rhs[i] = presolved_lp.upper[upper_bounds[k]]; + } + + barrier_lp.objective.resize(n + n_upper_bounds, 0.0); + barrier_lp.lower.resize(n + n_upper_bounds, 0.0); + barrier_lp.upper.resize(n + n_upper_bounds, inf); + + // Originally we have that + // A^T y + z_l - z_u = c + + // After the transformation we want + // A^T y_1 + y_2 + z_x = c_x + // y_2 + z_w = c_w = 0 + // So z_x = z_l + // And y_2 = -z_u + + std::vector zl = z; + std::vector zu = z; + for (i_t j = 0; j < n; j++) { + if (presolved_lp.upper[j] < inf) { + if (z[j] >= 0.0) { + zl[j] = z[j]; + zu[j] = 0.0; + } else { + zl[j] = 0.0; + zu[j] = -z[j]; + } + } else { + zl[j] = z[j]; + zu[j] = 0.0; + } + } + + z.resize(n + n_upper_bounds); + y.resize(presolved_lp.num_rows + n_upper_bounds); + + for (i_t k = 0; k < n_upper_bounds; k++) { + const i_t j = upper_bounds[k]; + y[presolved_lp.num_rows + k] = -zu[j]; + z[n + k] = zu[j]; + } + for (i_t j = 0; j < n; j++) { + z[j] = zl[j]; + } + + x.resize(n + n_upper_bounds); + for (i_t j = n; j < n + n_upper_bounds; j++) { + const i_t k = j - n; + x[j] = presolved_lp.upper[upper_bounds[k]] - x[upper_bounds[k]]; + } + + barrier_lp.num_cols = n + n_upper_bounds; + barrier_lp.num_rows = presolved_lp.num_rows + n_upper_bounds; + } + + std::vector res1 = barrier_lp.rhs; + matrix_vector_multiply(barrier_lp.A, 1.0, x, -1.0, res1); + f_t primal_inf1 = vector_norm_inf(res1); + // printf("Barrier Primal residual %e\n", primal_inf1); + + std::vector res2 = z; + for (size_t j = 0; j < z.size(); j++) { + res2[j] -= barrier_lp.objective[j]; + } + matrix_transpose_vector_multiply(barrier_lp.A, 1.0, y, +1.0, res2); + f_t dual_inf2 = vector_norm_inf(res2); + // printf("Barrier Dual residual %e\n", dual_inf2); + + f_t barrier_primal_inf = 0.0; + for (size_t j = 0; j < x.size(); j++) { + if (x[j] < 0.0) { barrier_primal_inf += -x[j]; } + } + + // printf("Barrier Primal infeasibility %e\n", barrier_primal_inf); + + f_t barrier_dual_inf = 0.0; + for (size_t j = 0; j < z.size(); j++) { + if (z[j] < 0.0) { barrier_dual_inf += -z[j]; } + } + + // printf("Barrier Dual infeasibility %e\n", barrier_dual_inf); + + f_t mu_avg = 0.0; + for (size_t j = 0; j < z.size(); j++) { + mu_avg += x[j] * z[j]; + } + mu_avg /= z.size(); + printf("Mu avg %e\n", mu_avg); + + // Now implement Algorithm 1 + n = barrier_lp.num_cols; + i_t m = barrier_lp.num_rows; + + f_t mu_min = 1e-6; + f_t delta_max = 1e-4; + f_t mu_target = 1e-6; + f_t alpha_min = 1e-6; + + std::vector xprime = x; + std::vector zprime = z; + + for (i_t j = 0; j < n; j++) { + xprime[j] = std::max(x[j], alpha_min); + zprime[j] = std::max(z[j], alpha_min); + if (x[j] >= z[j]) { + xprime[j] = + std::min(std::max(mu_target / zprime[j], xprime[j] - delta_max), xprime[j] + delta_max); + zprime[j] = + std::min(std::max(mu_target / xprime[j], zprime[j] - delta_max), zprime[j] + delta_max); + } else { + zprime[j] = + std::min(std::max(mu_target / xprime[j], zprime[j] - delta_max), zprime[j] + delta_max); + xprime[j] = + std::min(std::max(mu_target / zprime[j], xprime[j] - delta_max), xprime[j] + delta_max); + } + } + + bool skip_pdhg = false; + bool converged = false; + if (!skip_pdhg) { + // Now do PDHG on the central path problem + i_t iter_max = 1000000; + + f_t norm_b = vector_norm2(barrier_lp.rhs); + f_t norm_c = vector_norm2(barrier_lp.objective); + f_t primal_weight = norm_b > 0.0 ? norm_c / norm_b : 1.0; + f_t eta = 0.998 / barrier_lp.A.norm2_estimate(1e-4); + f_t tau = eta / primal_weight; + f_t sigma = eta * primal_weight; + f_t mu = mu_target; + + printf("tau %e sigma %e mu %e\n", tau, sigma, mu); + + std::vector ysav = y; + std::vector xbar = xprime; + std::vector xold = x; + std::vector v = x; + std::vector w = x; + + std::vector z_term = x; + + std::vector primal_residual = barrier_lp.rhs; + std::vector dual_residual = barrier_lp.objective; + + cusparse_view_t cusparse_view(lp.handle_ptr, barrier_lp.A); + rmm::device_uvector d_delta_y(m, lp.handle_ptr->get_stream()); + rmm::device_uvector d_xbar(n, lp.handle_ptr->get_stream()); + rmm::device_uvector d_y(m, lp.handle_ptr->get_stream()); + rmm::device_uvector d_x(n, lp.handle_ptr->get_stream()); + rmm::device_uvector d_xold(n, lp.handle_ptr->get_stream()); + rmm::device_uvector d_v(n, lp.handle_ptr->get_stream()); + rmm::device_uvector d_w(n, lp.handle_ptr->get_stream()); + rmm::device_uvector d_c(n, lp.handle_ptr->get_stream()); + rmm::device_uvector d_b(m, lp.handle_ptr->get_stream()); + rmm::device_uvector d_primal_residual(m, lp.handle_ptr->get_stream()); + rmm::device_uvector d_dual_residual(n, lp.handle_ptr->get_stream()); + + raft::copy(d_b.data(), barrier_lp.rhs.data(), m, lp.handle_ptr->get_stream()); + raft::copy(d_c.data(), barrier_lp.objective.data(), n, lp.handle_ptr->get_stream()); + + raft::copy(d_xbar.data(), xbar.data(), n, lp.handle_ptr->get_stream()); + raft::copy(d_y.data(), y.data(), m, lp.handle_ptr->get_stream()); + raft::copy(d_x.data(), x.data(), n, lp.handle_ptr->get_stream()); + + printf("%8s %12s %12s %8s %8s %8s\n", + "Iter", + "Primal Obj.", + "Dual Obj.", + "P. Resid.", + "D. Resid.", + "Time"); + for (i_t iter = 0; iter < iter_max; iter++) { + // Dual update + // delta_y = -sigma * A * xbar + sigma * b = sigma * (b - A * xbar) + raft::copy(d_delta_y.data(), d_b.data(), m, lp.handle_ptr->get_stream()); + cusparse_view.spmv(-1, d_xbar, 1, d_delta_y); + cub::DeviceTransform::Transform( + cuda::std::make_tuple(d_delta_y.data(), d_y.data()), + d_y.data(), + d_delta_y.size(), + [sigma] HD(f_t delta_y_i, f_t dy_i) -> f_t { return dy_i + delta_y_i * sigma; }, + lp.handle_ptr->get_stream().value()); + RAFT_CHECK_CUDA(lp.handle_ptr->get_stream()); + + raft::copy(d_xold.data(), d_x.data(), n, lp.handle_ptr->get_stream()); + + // xold = x + // Primal gradient step + // v = x + tau * A' * y + + lp.handle_ptr->sync_stream(); + cusparse_view.transpose_spmv(1.0, d_y, 0.0, d_v); + cub::DeviceTransform::Transform( + cuda::std::make_tuple(d_x.data(), d_v.data()), + d_v.data(), + d_v.size(), + [tau] HD(f_t x_i, f_t v_i) -> f_t { return x_i + tau * v_i; }, + lp.handle_ptr->get_stream().value()); + RAFT_CHECK_CUDA(lp.handle_ptr->get_stream()); + lp.handle_ptr->sync_stream(); + + cub::DeviceTransform::Transform( + cuda::std::make_tuple(d_v.data(), d_c.data(), d_x.data(), d_xold.data()), + thrust::make_zip_iterator(d_x.data(), d_xbar.data()), + d_x.size(), + [tau, mu] HD(f_t v_j, f_t c_j, f_t x_j, f_t xold_j) -> thrust::tuple { + f_t w_j = v_j - tau * c_j; + f_t new_x_j = (w_j + std::sqrt(w_j * w_j + 4.0 * tau * mu)) / 2.0; + return {new_x_j, 2.0 * new_x_j - xold_j}; + }, + lp.handle_ptr->get_stream().value()); + RAFT_CHECK_CUDA(lp.handle_ptr->get_stream()); + lp.handle_ptr->sync_stream(); + + if (iter % 1000 == 0) { + raft::copy(x.data(), d_x.data(), n, lp.handle_ptr->get_stream()); + raft::copy(d_primal_residual.data(), d_b.data(), m, lp.handle_ptr->get_stream()); + cusparse_view.spmv(1.0, d_x, -1.0, d_primal_residual); + + cub::DeviceTransform::Transform( + cuda::std::make_tuple(d_xbar.data(), d_c.data()), + d_dual_residual.data(), + d_xbar.size(), + [mu_target] HD(f_t xbar_j, f_t c_j) -> f_t { return mu_target / xbar_j - c_j; }, + lp.handle_ptr->get_stream().value()); + RAFT_CHECK_CUDA(lp.handle_ptr->get_stream()); + lp.handle_ptr->sync_stream(); + + cusparse_view.transpose_spmv(1.0, d_y, 1.0, d_dual_residual); + + f_t vector_norm_2_primal = + std::sqrt(thrust::inner_product(rmm::exec_policy(lp.handle_ptr->get_stream()), + d_primal_residual.data(), + d_primal_residual.data() + m, + d_primal_residual.data(), + 0.0)); + f_t vector_norm_2_dual = + std::sqrt(thrust::inner_product(rmm::exec_policy(lp.handle_ptr->get_stream()), + d_dual_residual.data(), + d_dual_residual.data() + n, + d_dual_residual.data(), + 0.0)); + + f_t primal_relative = vector_norm_2_primal / (1.0 + norm_b); + f_t dual_relative = vector_norm_2_dual / (1.0 + norm_c); + + f_t primal_obj = thrust::inner_product(rmm::exec_policy(lp.handle_ptr->get_stream()), + d_c.data(), + d_c.data() + n, + d_xbar.data(), + 0.0); + f_t dual_obj = thrust::inner_product(rmm::exec_policy(lp.handle_ptr->get_stream()), + d_b.data(), + d_b.data() + m, + d_y.data(), + 0.0); + f_t time = toc(start_time); + printf("%8d %+12.6e %+12.6e %8.2e %8.2e %8.1fs\n", + iter, + primal_obj, + dual_obj, + primal_relative, + dual_relative, + time); + + if (primal_relative < 1e-5 && dual_relative < 1e-5) { + converged = true; + raft::copy(x.data(), d_x.data(), n, lp.handle_ptr->get_stream()); + raft::copy(y.data(), d_y.data(), m, lp.handle_ptr->get_stream()); + raft::copy(xbar.data(), d_xbar.data(), n, lp.handle_ptr->get_stream()); + break; + } + } + } + } else { + x = xprime; + z = zprime; + converged = true; + } + + if (converged) { + printf("PDHG found point on central path\n\n\n"); + // Solve using barrier + lp_solution_t barrier_solution(presolved_lp.num_rows, presolved_lp.num_cols); + + barrier_solver_t barrier_solver(presolved_lp, presolve_info, barrier_settings); + + // Convert the converged solution to (x, w, y, z, v) + std::vector x_pdhg(n - n_upper_bounds); + std::vector w_pdhg(n_upper_bounds); + std::vector y_pdhg(m - n_upper_bounds); + std::vector z_pdhg(n - n_upper_bounds); + std::vector v_pdhg(n_upper_bounds); + + for (i_t j = 0; j < n - n_upper_bounds; j++) { + x_pdhg[j] = x[j]; + z_pdhg[j] = mu_target / x[j]; + } + + for (i_t j = 0; j < n_upper_bounds; j++) { + w_pdhg[j] = x[n - n_upper_bounds + j]; + v_pdhg[j] = mu_target / x[n - n_upper_bounds + j]; + } + + for (i_t i = 0; i < m - n_upper_bounds; i++) { + y_pdhg[i] = y[i]; + } + + barrier_solver.x_pdhg = x_pdhg; + barrier_solver.w_pdhg = w_pdhg; + barrier_solver.y_pdhg = y_pdhg; + barrier_solver.z_pdhg = z_pdhg; + barrier_solver.v_pdhg = v_pdhg; + + barrier_solver_settings_t barrier_solver_settings; + lp_status_t barrier_status = + barrier_solver.solve(start_time, barrier_solver_settings, barrier_solution); + + return barrier_status == lp_status_t::OPTIMAL ? crossover_status_t::OPTIMAL + : crossover_status_t::NUMERICAL_ISSUES; + } else { + return crossover_status_t::NUMERICAL_ISSUES; + } +} + #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE template crossover_status_t crossover( @@ -1582,6 +2026,14 @@ template crossover_status_t crossover( lp_solution_t& solution, std::vector& vstatus); + template crossover_status_t central_path( + const lp_problem_t& problem, + const simplex_solver_settings_t& settings, + const lp_solution_t& initial_solution, + double start_time, + lp_solution_t& solution, + std::vector& vstatus); + #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/crossover.hpp b/cpp/src/dual_simplex/crossover.hpp index 0bb1ff6627..e30ebe2254 100644 --- a/cpp/src/dual_simplex/crossover.hpp +++ b/cpp/src/dual_simplex/crossover.hpp @@ -32,4 +32,12 @@ crossover_status_t crossover(const lp_problem_t& problem, lp_solution_t& solution, std::vector& vstatus); +template +crossover_status_t central_path(const lp_problem_t& problem, + const simplex_solver_settings_t& settings, + const lp_solution_t& initial_solution, + f_t start_time, + lp_solution_t& solution, + std::vector& vstatus); + } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index b9ee419517..71dc1dd3fe 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -835,6 +835,24 @@ i_t presolve(const lp_problem_t& original, } // FIXME:: handle no lower bound case for barrier presolve + 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; + } + } + } + } // The original problem may have nonzero lower bounds // 0 != l_j <= x_j <= u_j @@ -1134,6 +1152,58 @@ i_t presolve(const lp_problem_t& original, return 0; } +template +void crush_solution_to_presolve_space(const lp_problem_t& original_lp, + presolve_info_t& presolve_info, + const std::vector& original_x, + const std::vector& original_y, + const std::vector& original_z, + std::vector& crushed_x, + std::vector& crushed_y, + std::vector& crushed_z) +{ + crushed_x = original_x; + crushed_y = original_y; + crushed_z = original_z; + + if (presolve_info.negated_variables.size() > 0) { + for (i_t j = 0; j < presolve_info.negated_variables.size(); j++) { + crushed_x[presolve_info.negated_variables[j]] *= -1.0; + crushed_z[presolve_info.negated_variables[j]] *= -1.0; + } + } + + // Repeat all presolve steps to get the crushed solution + if (presolve_info.removed_lower_bounds.size() > 0) { + // We had l_j <= x_j <= u_j + // And we transformed it into + // 0 <= v_j = x_j - l_j + // So we need to subtract the removed lower bounds from the original solution + for (i_t j = 0; j < original_x.size(); j++) { + crushed_x[j] -= presolve_info.removed_lower_bounds[j]; + } + } + + // Handle empty rows + if (presolve_info.removed_constraints.size() > 0) { + printf("Can't handle removed constraints %ld\n", presolve_info.removed_constraints.size()); + exit(1); + } + + // Handle empty cols + if (presolve_info.removed_variables.size() > 0) { + printf("Can't handle removed variables %ld\n", presolve_info.removed_variables.size()); + exit(1); + } + + // Handle free variables + if (presolve_info.free_variable_pairs.size() > 0) { + printf("Can't handle free variable pairs %ld\n", presolve_info.free_variable_pairs.size()); + exit(1); + } + +} + template void convert_user_lp_with_guess(const user_problem_t& user_problem, const lp_solution_t& initial_solution, @@ -1569,6 +1639,15 @@ template void uncrush_solution(const presolve_info_t& std::vector& uncrushed_y, std::vector& uncrushed_z); +template void crush_solution_to_presolve_space( + const lp_problem_t& problem, + presolve_info_t& presolve_info, + const std::vector& original_x, + const std::vector& original_y, + const std::vector& original_z, + std::vector& x, + std::vector& y, + std::vector& z); #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 17e6176e3b..b6a621b626 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -87,6 +87,9 @@ struct presolve_info_t { // problem std::vector removed_constraints; + // Variables that were negated to handle no lower bound + std::vector negated_variables; + folding_info_t folding_info; }; @@ -183,4 +186,14 @@ void uncrush_solution(const presolve_info_t& presolve_info, std::vector& uncrushed_y, std::vector& uncrushed_z); +template +void crush_solution_to_presolve_space(const lp_problem_t& problem, + presolve_info_t& presolve_info, + const std::vector& original_x, + const std::vector& original_y, + const std::vector& original_z, + std::vector& x, + std::vector& y, + std::vector& z); + } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/random.hpp b/cpp/src/dual_simplex/random.hpp index 571a59857a..4898f022dd 100644 --- a/cpp/src/dual_simplex/random.hpp +++ b/cpp/src/dual_simplex/random.hpp @@ -30,6 +30,13 @@ class random_t { return distrib(gen); } + template + T random_value(T min, T max) + { + std::uniform_real_distribution<> distrib(min, max); + return distrib(gen); + } + private: std::mt19937 gen; }; diff --git a/cpp/src/dual_simplex/sparse_matrix.cpp b/cpp/src/dual_simplex/sparse_matrix.cpp index 8ccccd57cf..aa040f4455 100644 --- a/cpp/src/dual_simplex/sparse_matrix.cpp +++ b/cpp/src/dual_simplex/sparse_matrix.cpp @@ -8,7 +8,7 @@ // #include #include #include - +#include #include // #include @@ -845,6 +845,60 @@ f_t csc_matrix_t::norm1() const return norm; } +template +f_t csc_matrix_t::norm2_estimate(f_t tol) const +{ + i_t m = this->m; + i_t n = this->n; + f_t norm = 0.0; + + std::vector x(n); + std::vector Sx(m); + + for (i_t j = 0; j < n; ++j) { + const i_t col_start = this->col_start[j]; + const i_t col_end = this->col_start[j + 1]; + x[j] = 0.0; + for (i_t p = col_start; p < col_end; ++p) { + x[j] += std::abs(this->x[p]); + } + } + + f_t e = vector_norm2(x); + if (e == 0.0) { return 0.0; } + + for (i_t j = 0; j < n; ++j) { + x[j] /= e; + } + + f_t e0 = 0.0; + + i_t iter = 0; + const i_t max_iter = 100; + while (std::abs(e - e0) > tol * e) { + e0 = e; + matrix_vector_multiply(*this, 1.0, x, 0.0, Sx); + f_t Sx_norm = vector_norm2(Sx); + if (Sx_norm == 0.0) { + random_t rng(0); + for (i_t i = 0; i < m; ++i) { + Sx[i] = rng.random_value(0.0, 1.0); + } + Sx_norm = vector_norm2(Sx); + } + matrix_transpose_vector_multiply(*this, 1.0, Sx, 0.0, x); + f_t norm_x = vector_norm2(x); + e = norm_x / Sx_norm; + for (i_t j = 0; j < n; ++j) { + x[j] /= norm_x; + } + + iter++; + if (iter > max_iter) { break; } + } + return e; +} + template f_t vector_norm_inf(const std::vector& x) { diff --git a/cpp/src/dual_simplex/sparse_matrix.hpp b/cpp/src/dual_simplex/sparse_matrix.hpp index b6d3ee9aae..00b28f0aba 100644 --- a/cpp/src/dual_simplex/sparse_matrix.hpp +++ b/cpp/src/dual_simplex/sparse_matrix.hpp @@ -103,6 +103,8 @@ class csc_matrix_t { // Compute || A ||_1 = max_j (sum {i = 1 to m} | A(i, j) | ) f_t norm1() const; + f_t norm2_estimate(f_t tol) const; + // Compare two matrices void compare(csc_matrix_t const& B) const; diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 5e1e25bbee..42dd2e628a 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -620,7 +620,7 @@ optimization_problem_solution_t run_pdlp(detail::problem_t& dual_simplex_settings.concurrent_halt = settings.concurrent_halt; dual_simplex::lp_solution_t vertex_solution(lp.num_rows, lp.num_cols); std::vector vstatus(lp.num_cols); - dual_simplex::crossover_status_t crossover_status = dual_simplex::crossover( + dual_simplex::crossover_status_t crossover_status = dual_simplex::central_path( lp, dual_simplex_settings, initial_solution, timer.get_tic_start(), vertex_solution, vstatus); pdlp_termination_status_t termination_status = pdlp_termination_status_t::TimeLimit; auto to_termination_status = [](dual_simplex::crossover_status_t status) { diff --git a/cpp/src/pdlp/translate.hpp b/cpp/src/pdlp/translate.hpp index b8e0075733..440162eb48 100644 --- a/cpp/src/pdlp/translate.hpp +++ b/cpp/src/pdlp/translate.hpp @@ -11,6 +11,7 @@ #include #include +#include #include @@ -190,12 +191,18 @@ void translate_to_crossover_problem(const detail::problem_t& problem, CUOPT_LOG_DEBUG("Finished with x"); initial_solution.y = cuopt::host_copy(sol.get_dual_solution(), stream); - std::vector tmp_z = cuopt::host_copy(sol.get_reduced_cost(), stream); - stream.synchronize(); - std::copy(tmp_z.begin(), tmp_z.begin() + problem.n_variables, initial_solution.z.begin()); - for (i_t j = problem.n_variables; j < n; ++j) { - initial_solution.z[j] = initial_solution.y[j - problem.n_variables]; + // Compute z = c - A^T y (exact dual feasibility) +std::copy(lp.objective.begin(), lp.objective.end(), initial_solution.z.begin()); +dual_simplex::matrix_transpose_vector_multiply( + lp.A, -1.0, initial_solution.y, 1.0, initial_solution.z); + + std::vector dual_residual = initial_solution.z; + for (i_t j = 0; j < lp.num_cols; j++) { + dual_residual[j] -= lp.objective[j]; } + + dual_simplex::matrix_transpose_vector_multiply(lp.A, 1.0, initial_solution.y, +1.0, dual_residual); + CUOPT_LOG_DEBUG("Finished with z"); CUOPT_LOG_DEBUG("Finished translating");