diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 35f7b7e08..03e6d77ce 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -85,6 +85,8 @@ endif(BUILD_SANITIZER) if(DEFINE_ASSERT) add_definitions(-DASSERT_MODE) + list(APPEND CUOPT_CXX_FLAGS -ftrapping-math) + list(APPEND CUOPT_CUDA_FLAGS -Xcompiler=-ftrapping-math) endif(DEFINE_ASSERT) if(DEFINE_BENCHMARK) diff --git a/cpp/src/dual_simplex/basis_solves.cpp b/cpp/src/dual_simplex/basis_solves.cpp index db24f55a2..d6992bf84 100644 --- a/cpp/src/dual_simplex/basis_solves.cpp +++ b/cpp/src/dual_simplex/basis_solves.cpp @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -363,7 +363,7 @@ i_t factorize_basis(const csc_matrix_t& A, S_perm_inv); if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Concurrent halt\n"); - return -1; + return -2; // Use -2 to distinguish from rank deficiency (-1) } if (Srank != Sdim) { // Get the rank deficient columns @@ -582,7 +582,7 @@ i_t factorize_basis(const csc_matrix_t& A, } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Concurrent halt\n"); - return -1; + return -2; // Use -2 to distinguish from rank deficiency (-1) } if (verbose) { printf("Right Lnz+Unz %d t %.3f\n", L.col_start[m] + U.col_start[m], toc(fact_start)); diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 6b79f3c86..7b4a4735c 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -1354,7 +1354,8 @@ i_t basis_update_mpf_t::b_transpose_solve(const sparse_vector_t(solution.i.size()) / input_size; + // avoid div-by-zero FPE + if (input_size > 0.0) { sum_U_transpose_ += static_cast(solution.i.size()) / input_size; } #ifdef CHECK_U_TRANSPOSE_SOLVE std::vector UTsol_dense; @@ -1383,7 +1384,8 @@ i_t basis_update_mpf_t::b_transpose_solve(const sparse_vector_t(solution.i.size()) / rhs_size; + // avoid div-by-zero FPE + if (rhs_size > 0.0) { sum_L_transpose_ += static_cast(solution.i.size()) / rhs_size; } #ifdef CHECK_L_TRANSPOSE_SOLVE std::vector solution_dense; @@ -1656,7 +1658,7 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, solution.from_dense(solution_dense); } if (need_Lsol) { Lsol = solution; } - sum_L_ += static_cast(solution.i.size()) / input_size; + if (input_size > 0.0) { sum_L_ += static_cast(solution.i.size()) / input_size; } #ifdef CHECK_L_SOLVE std::vector l_solve_dense; @@ -1677,6 +1679,8 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, #endif const f_t rhs_size = static_cast(solution.i.size()); + // (nothing to solve, avoids 0/0 in statistics) + if (rhs_size == 0.0) { return 0; } estimate_solution_density(rhs_size, sum_U_, num_calls_U_, use_hypersparse); if (use_hypersparse) { u_solve(solution); @@ -2055,16 +2059,21 @@ int basis_update_mpf_t::refactor_basis( if (L0_.m != A.m) { resize(A.m); } std::vector q; - if (factorize_basis(A, - settings, - basic_list, - L0_, - U0_, - row_permutation_, - inverse_row_permutation_, - q, - deficient, - slacks_needed) == -1) { + i_t factorize_result = factorize_basis(A, + settings, + basic_list, + L0_, + U0_, + row_permutation_, + inverse_row_permutation_, + q, + deficient, + slacks_needed); + if (factorize_result == -2) { + // Concurrent halt requested, return early + return -2; + } + if (factorize_result == -1) { settings.log.debug("Initial factorization failed\n"); basis_repair(A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus); @@ -2085,16 +2094,18 @@ int basis_update_mpf_t::refactor_basis( } #endif - if (factorize_basis(A, - settings, - basic_list, - L0_, - U0_, - row_permutation_, - inverse_row_permutation_, - q, - deficient, - slacks_needed) == -1) { + factorize_result = factorize_basis(A, + settings, + basic_list, + L0_, + U0_, + row_permutation_, + inverse_row_permutation_, + q, + deficient, + slacks_needed); + if (factorize_result == -2) { return -2; } + if (factorize_result == -1) { #ifdef CHECK_L_FACTOR if (L0_.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); } #endif diff --git a/cpp/src/dual_simplex/bounds_strengthening.cpp b/cpp/src/dual_simplex/bounds_strengthening.cpp index f1bf52c1e..826d901c0 100644 --- a/cpp/src/dual_simplex/bounds_strengthening.cpp +++ b/cpp/src/dual_simplex/bounds_strengthening.cpp @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -10,6 +10,23 @@ #include #include +#include + +class fpe_disable { + int old_mask; + + public: + explicit fpe_disable(int mask = FE_INVALID) : old_mask(fegetexcept()) { fedisableexcept(mask); } + ~fpe_disable() + { + fedisableexcept(FE_ALL_EXCEPT); + feenableexcept(old_mask); + } + + fpe_disable(const fpe_disable&) = delete; + fpe_disable& operator=(const fpe_disable&) = delete; +}; + namespace cuopt::linear_programming::dual_simplex { template @@ -95,6 +112,10 @@ bool bounds_strengthening_t::bounds_strengthening( std::vector& upper_bounds, const simplex_solver_settings_t& settings) { + // Disable FPEs to avoid inf + (-inf) issues + // TODO: actually resolve the numerical issues + fpe_disable fpe_guard(FE_INVALID); + const i_t m = A.m; const i_t n = A.n; diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 6161f4d3f..4ddaf4ac7 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -159,6 +159,10 @@ f_t sgn(f_t x) template f_t relative_gap(f_t obj_value, f_t lower_bound) { + // avoid raising a FPE + if (!std::isfinite(obj_value) || !std::isfinite(lower_bound)) { + return std::numeric_limits::infinity(); + } f_t user_mip_gap = obj_value == 0.0 ? (lower_bound == 0.0 ? 0.0 : std::numeric_limits::infinity()) : std::abs(obj_value - lower_bound) / std::abs(obj_value); @@ -171,9 +175,13 @@ f_t user_relative_gap(const lp_problem_t& lp, f_t obj_value, f_t lower { f_t user_obj = compute_user_objective(lp, obj_value); f_t user_lower_bound = compute_user_objective(lp, lower_bound); - f_t user_mip_gap = user_obj == 0.0 - ? (user_lower_bound == 0.0 ? 0.0 : std::numeric_limits::infinity()) - : std::abs(user_obj - user_lower_bound) / std::abs(user_obj); + // avoid raising a FPE + if (!std::isfinite(user_obj) || !std::isfinite(user_lower_bound)) { + return std::numeric_limits::infinity(); + } + f_t user_mip_gap = user_obj == 0.0 + ? (user_lower_bound == 0.0 ? 0.0 : std::numeric_limits::infinity()) + : std::abs(user_obj - user_lower_bound) / std::abs(user_obj); if (std::isnan(user_mip_gap)) { return std::numeric_limits::infinity(); } return user_mip_gap; } @@ -463,7 +471,9 @@ mip_status_t branch_and_bound_t::set_final_solution(mip_solution_t::infinity(); f_t obj = compute_user_objective(original_lp_, upper_bound); f_t user_bound = compute_user_objective(original_lp_, lower_bound); f_t gap_rel = user_relative_gap(original_lp_, upper_bound, lower_bound); @@ -1354,6 +1364,11 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::UNBOUNDED; } + if (root_status == lp_status_t::CONCURRENT_LIMIT && + toc(exploration_stats_.start_time) > settings_.time_limit) { + solver_status_ = mip_exploration_status_t::TIME_LIMIT; + return set_final_solution(solution, root_objective_); + } if (root_status == lp_status_t::TIME_LIMIT) { solver_status_ = mip_exploration_status_t::TIME_LIMIT; return set_final_solution(solution, -inf); diff --git a/cpp/src/dual_simplex/crossover.cpp b/cpp/src/dual_simplex/crossover.cpp index 23d9a0e8e..d2b4466ba 100644 --- a/cpp/src/dual_simplex/crossover.cpp +++ b/cpp/src/dual_simplex/crossover.cpp @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -783,12 +783,15 @@ i_t primal_push(const lp_problem_t& lp, std::vector slacks_needed; i_t rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); + if (rank == -2) { return -2; } // Concurrent halt if (rank != m) { settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m); basis_repair( lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus); - if (factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == -1) { + i_t repair_rank = + factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); + if (repair_rank == -2) { return -2; } // Concurrent halt + if (repair_rank == -1) { settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m); return -1; } else { @@ -1130,11 +1133,14 @@ crossover_status_t crossover(const lp_problem_t& lp, std::vector slacks_needed; rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); + if (rank == -2) { return crossover_status_t::CONCURRENT_LIMIT; } // Concurrent halt if (rank != m) { settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m); basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus); - if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == - -1) { + i_t repair_rank = + factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); + if (repair_rank == -2) { return crossover_status_t::CONCURRENT_LIMIT; } // Concurrent halt + if (repair_rank == -1) { settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m); return crossover_status_t::NUMERICAL_ISSUES; } else { @@ -1321,11 +1327,14 @@ crossover_status_t crossover(const lp_problem_t& lp, get_basis_from_vstatus(m, vstatus, basic_list, nonbasic_list, superbasic_list); rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); + if (rank == -2) { return crossover_status_t::CONCURRENT_LIMIT; } // Concurrent halt if (rank != m) { settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m); basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus); - if (factorize_basis( - lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == -1) { + i_t repair_rank = + factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); + if (repair_rank == -2) { return crossover_status_t::CONCURRENT_LIMIT; } // Concurrent halt + if (repair_rank == -1) { settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m); return crossover_status_t::NUMERICAL_ISSUES; } else { diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 56298ef4d..097ccc913 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -767,7 +767,8 @@ i_t steepest_edge_pricing_with_infeasibilities(const lp_problem_t& lp, for (i_t k = 0; k < nz; ++k) { const i_t j = infeasibility_indices[k]; const f_t squared_infeas = squared_infeasibilities[j]; - const f_t val = squared_infeas / dy_steepest_edge[j]; + const f_t val = dy_steepest_edge[j] == 0.0 ? std::numeric_limits::infinity() + : squared_infeas / dy_steepest_edge[j]; if (val > max_val || (val == max_val && j > leaving_index)) { max_val = val; leaving_index = j; @@ -2241,9 +2242,9 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, assert(superbasic_list.size() == 0); assert(nonbasic_list.size() == n - m); - if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) { - return dual::status_t::NUMERICAL; - } + i_t refactor_result = ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus); + if (refactor_result == -2) { return dual::status_t::CONCURRENT_LIMIT; } + if (refactor_result > 0) { return dual::status_t::NUMERICAL; } if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; } } @@ -2883,7 +2884,9 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, #endif if (should_refactor) { bool should_recompute_x = false; - if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) { + i_t refactor_result = ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus); + if (refactor_result == -2) { return dual::status_t::CONCURRENT_LIMIT; } + if (refactor_result > 0) { should_recompute_x = true; settings.log.printf("Failed to factorize basis. Iteration %d\n", iter); if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; } @@ -2901,6 +2904,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, count++; if (count > 10) { return dual::status_t::NUMERICAL; } } + if (deficient_size == -2) { return dual::status_t::CONCURRENT_LIMIT; } settings.log.printf("Successfully repaired basis. Iteration %d\n", iter); } diff --git a/cpp/src/dual_simplex/primal.cpp b/cpp/src/dual_simplex/primal.cpp index 80406dcf0..f9ca2c5e5 100644 --- a/cpp/src/dual_simplex/primal.cpp +++ b/cpp/src/dual_simplex/primal.cpp @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -296,11 +296,14 @@ primal::status_t primal_phase2(i_t phase, std::vector slacks_needed; i_t rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); + if (rank == -2) { return primal::status_t::CONCURRENT_HALT; } if (rank != m) { settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m); basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus); - if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == - -1) { + i_t repair_rank = + factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); + if (repair_rank == -2) { return primal::status_t::CONCURRENT_HALT; } + if (repair_rank == -1) { settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m); return primal::status_t::NUMERICAL; } else { diff --git a/cpp/src/dual_simplex/primal.hpp b/cpp/src/dual_simplex/primal.hpp index a5d356fdb..cd4b0cce4 100644 --- a/cpp/src/dual_simplex/primal.hpp +++ b/cpp/src/dual_simplex/primal.hpp @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -23,7 +23,8 @@ enum class status_t { PRIMAL_UNBOUNDED = 1, NUMERICAL = 2, NOT_LOADED = 3, - ITERATION_LIMIT = 4 + ITERATION_LIMIT = 4, + CONCURRENT_HALT = 5 }; } diff --git a/cpp/src/linear_programming/restart_strategy/pdlp_restart_strategy.cu b/cpp/src/linear_programming/restart_strategy/pdlp_restart_strategy.cu index 565377682..4e7449289 100644 --- a/cpp/src/linear_programming/restart_strategy/pdlp_restart_strategy.cu +++ b/cpp/src/linear_programming/restart_strategy/pdlp_restart_strategy.cu @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2022-2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -729,9 +729,13 @@ void pdlp_restart_strategy_t::cupdlpx_restart( current_convergence_information.get_relative_l2_dual_residual_value(); const f_t relative_l2_primal_residual_value = current_convergence_information.get_relative_l2_primal_residual_value(); - const f_t ratio_infeas = (relative_l2_primal_residual_value == f_t(0.0)) - ? std::numeric_limits::infinity() - : relative_l2_dual_residual_value / relative_l2_primal_residual_value; + + f_t ratio_infeas; + if (relative_l2_primal_residual_value < 1e-16) { + ratio_infeas = std::numeric_limits::infinity(); + } else { + ratio_infeas = relative_l2_dual_residual_value / relative_l2_primal_residual_value; + } if (l2_primal_distance > f_t(1e-16) && l2_dual_distance > f_t(1e-16) && l2_primal_distance < f_t(1e12) && l2_dual_distance < f_t(1e12) && ratio_infeas > f_t(1e-8) && @@ -774,7 +778,9 @@ void pdlp_restart_strategy_t::cupdlpx_restart( // TODO possible error if relative_l2_primal_residual_value == 0 const f_t primal_dual_residual_gap = - std::abs(std::log10(relative_l2_dual_residual_value / relative_l2_primal_residual_value)); + relative_l2_primal_residual_value > f_t(1e-16) + ? std::abs(std::log10(relative_l2_dual_residual_value / relative_l2_primal_residual_value)) + : std::numeric_limits::infinity(); if (primal_dual_residual_gap < best_primal_dual_residual_gap_) { best_primal_dual_residual_gap_ = primal_dual_residual_gap; const f_t new_primal_weight = primal_weight.value(stream_view_); diff --git a/cpp/src/linear_programming/solve.cu b/cpp/src/linear_programming/solve.cu index d038ade72..a2c945218 100644 --- a/cpp/src/linear_programming/solve.cu +++ b/cpp/src/linear_programming/solve.cu @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2022-2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -42,8 +42,25 @@ #include // For std::thread +#include + namespace cuopt::linear_programming { +class fpe_enable { + int old_mask; + + public: + explicit fpe_enable(int mask = FE_INVALID) : old_mask(fegetexcept()) { feenableexcept(mask); } + ~fpe_enable() + { + fedisableexcept(FE_ALL_EXCEPT); + feenableexcept(old_mask); + } + + fpe_enable(const fpe_enable&) = delete; + fpe_enable& operator=(const fpe_enable&) = delete; +}; + // This serves as both a warm up but also a mandatory initial call to setup cuSparse and cuBLAS static void init_handler(const raft::handle_t* handle_ptr) { @@ -555,7 +572,9 @@ optimization_problem_solution_t run_pdlp(detail::problem_t& if (sol.get_termination_status() != pdlp_termination_status_t::ConcurrentLimit) { CUOPT_LOG_INFO("Status: %s Objective: %.8e Iterations: %d Time: %.3fs, Total time %.3fs", sol.get_termination_status_string().c_str(), - sol.get_objective_value(), + // printf doesn't like signaling NaNs for some reason + std::isnan(sol.get_objective_value()) ? std::numeric_limits::quiet_NaN() + : sol.get_objective_value(), sol.get_additional_termination_information().number_of_steps_taken, pdlp_solve_time, sol.get_solve_time()); @@ -827,6 +846,11 @@ optimization_problem_solution_t solve_lp( // This needs to be called before pdlp is initialized init_handler(op_problem.get_handle_ptr()); +#if 1 + CUOPT_LOG_DEBUG("Enabling host FPEs"); + fpe_enable fpe_guard(FE_DIVBYZERO | FE_INVALID); +#endif + if (op_problem.has_quadratic_objective()) { CUOPT_LOG_INFO("Problem has a quadratic objective. Using Barrier."); settings.method = method_t::Barrier; diff --git a/cpp/src/mip/diversity/diversity_manager.cu b/cpp/src/mip/diversity/diversity_manager.cu index 483ffeb68..d53977ad5 100644 --- a/cpp/src/mip/diversity/diversity_manager.cu +++ b/cpp/src/mip/diversity/diversity_manager.cu @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2024-2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2024-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -432,15 +432,14 @@ solution_t diversity_manager_t::run_solver() problem_ptr->handle_ptr->get_stream()); problem_ptr->handle_ptr->sync_stream(); - auto user_obj = problem_ptr->get_user_obj_from_solver_obj(lp_result.get_objective_value()); + // PDLP returns user-space objective (it applies objective_scaling_factor internally) + auto user_obj = lp_result.get_objective_value(); + auto solver_obj = problem_ptr->get_solver_obj_from_user_obj(user_obj); auto iterations = lp_result.get_additional_termination_information().number_of_steps_taken; - // Set for the B&B - problem_ptr->set_root_relaxation_solution_callback(host_primal, - host_dual, - host_reduced_costs, - lp_result.get_objective_value(), - user_obj, - iterations); + + // Set for the B&B (param4 expects solver space, param5 expects user space) + problem_ptr->set_root_relaxation_solution_callback( + host_primal, host_dual, host_reduced_costs, solver_obj, user_obj, iterations); } // in case the pdlp returned var boudns that are out of bounds diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh b/cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh index fbc5a7b39..6d76a31dd 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh +++ b/cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -166,12 +166,12 @@ HDI std::pair feas_score_constraint( // simple improvement else if (!old_sat && !new_sat && old_lhs > new_lhs) { cuopt_assert(old_viol && new_viol, ""); - base_feas += (i_t)(cstr_weight * fj.settings->parameters.excess_improvement_weight); + base_feas += round(cstr_weight * fj.settings->parameters.excess_improvement_weight); } // simple worsening else if (!old_sat && !new_sat && old_lhs <= new_lhs) { cuopt_assert(old_viol && new_viol, ""); - base_feas -= (i_t)(cstr_weight * fj.settings->parameters.excess_improvement_weight); + base_feas -= round(cstr_weight * fj.settings->parameters.excess_improvement_weight); } // robustness score bonus if this would leave some strick slack diff --git a/cpp/src/mip/feasibility_jump/fj_cpu.cu b/cpp/src/mip/feasibility_jump/fj_cpu.cu index 8d534dfff..013a380a7 100644 --- a/cpp/src/mip/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip/feasibility_jump/fj_cpu.cu @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -938,15 +938,17 @@ bool fj_t::cpu_solve(fj_cpu_climber_t& fj_cpu, f_t in_time_l { raft::common::nvtx::range scope("fj_cpu"); + auto time_limit_ms = + std::isfinite(in_time_limit) ? (int)(in_time_limit * 1000) : std::numeric_limits::max(); + i_t local_mins = 0; auto loop_start = std::chrono::high_resolution_clock::now(); - auto time_limit = std::chrono::milliseconds((int)(in_time_limit * 1000)); + auto time_limit = std::chrono::milliseconds(time_limit_ms); auto loop_time_start = std::chrono::high_resolution_clock::now(); while (!fj_cpu.halted && !fj_cpu.preemption_flag.load()) { // Check if 5 seconds have passed auto now = std::chrono::high_resolution_clock::now(); - if (in_time_limit < std::numeric_limits::infinity() && - now - loop_time_start > time_limit) { + if (time_limit_ms != std::numeric_limits::max() && now - loop_time_start > time_limit) { CUOPT_LOG_TRACE("%sTime limit of %.4f seconds reached, breaking loop at iteration %d\n", fj_cpu.log_prefix.c_str(), time_limit.count() / 1000.f, diff --git a/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu b/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu index 76d930b45..35b5c3dcf 100644 --- a/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu +++ b/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2024-2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2024-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -85,9 +85,10 @@ void feasibility_pump_t::adjust_objective_with_original(solution_t::get_user_obj_from_solver_obj(f_t solver_obj) const return presolve_data.objective_scaling_factor * (solver_obj + presolve_data.objective_offset); } +template +f_t problem_t::get_solver_obj_from_user_obj(f_t user_obj) const +{ + return (user_obj / presolve_data.objective_scaling_factor) - presolve_data.objective_offset; +} + template void problem_t::compute_vars_with_objective_coeffs() { diff --git a/cpp/src/mip/problem/problem.cuh b/cpp/src/mip/problem/problem.cuh index ed0adb971..3e81fc610 100644 --- a/cpp/src/mip/problem/problem.cuh +++ b/cpp/src/mip/problem/problem.cuh @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2024-2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2024-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -91,6 +91,7 @@ class problem_t { void post_process_solution(solution_t& solution); void compute_transpose_of_problem(); f_t get_user_obj_from_solver_obj(f_t solver_obj) const; + f_t get_solver_obj_from_user_obj(f_t user_obj) const; bool is_objective_integral() const { return objective_is_integral; } void compute_integer_fixed_problem(); void fill_integer_fixed_problem(rmm::device_uvector& assignment, diff --git a/cpp/src/mip/solve.cu b/cpp/src/mip/solve.cu index e6a392d40..b3bba34a4 100644 --- a/cpp/src/mip/solve.cu +++ b/cpp/src/mip/solve.cu @@ -38,8 +38,25 @@ #include +#include + namespace cuopt::linear_programming { +class fpe_enable { + int old_mask; + + public: + explicit fpe_enable(int mask = FE_INVALID) : old_mask(fegetexcept()) { feenableexcept(mask); } + ~fpe_enable() + { + fedisableexcept(FE_ALL_EXCEPT); + feenableexcept(old_mask); + } + + fpe_enable(const fpe_enable&) = delete; + fpe_enable& operator=(const fpe_enable&) = delete; +}; + // This serves as both a warm up but also a mandatory initial call to setup cuSparse and cuBLAS static void init_handler(const raft::handle_t* handle_ptr) { @@ -160,6 +177,11 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, // This needs to be called before pdlp is initialized init_handler(op_problem.get_handle_ptr()); +#if 1 + CUOPT_LOG_DEBUG("Enabling host FPEs"); + fpe_enable fpe_guard(FE_DIVBYZERO | FE_INVALID); +#endif + print_version_info(); raft::common::nvtx::range fun_scope("Running solver"); diff --git a/cpp/tests/linear_programming/pdlp_test.cu b/cpp/tests/linear_programming/pdlp_test.cu index 01e10ca61..bd298feee 100644 --- a/cpp/tests/linear_programming/pdlp_test.cu +++ b/cpp/tests/linear_programming/pdlp_test.cu @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2022-2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -55,856 +55,884 @@ static bool is_incorrect_objective(double reference, double objective) return std::abs((reference - objective) / reference) > 0.01; } -TEST(pdlp_class, run_double) -{ - const raft::handle_t handle_{}; - - auto path = make_path_absolute("linear_programming/afiro_original.mps"); - cuopt::mps_parser::mps_data_model_t op_problem = - cuopt::mps_parser::parse_mps(path, true); - - auto solver_settings = pdlp_solver_settings_t{}; - solver_settings.method = cuopt::linear_programming::method_t::PDLP; - - optimization_problem_solution_t solution = - solve_lp(&handle_, op_problem, solver_settings); - EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); - EXPECT_FALSE(is_incorrect_objective( - afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); -} - -TEST(pdlp_class, run_double_very_low_accuracy) -{ - const raft::handle_t handle_{}; - - auto path = make_path_absolute("linear_programming/afiro_original.mps"); - cuopt::mps_parser::mps_data_model_t op_problem = - cuopt::mps_parser::parse_mps(path, true); - - cuopt::linear_programming::pdlp_solver_settings_t settings = - cuopt::linear_programming::pdlp_solver_settings_t{}; - // With all 0 afiro with return an error - // Setting absolute tolerance to the minimal value of 1e-12 will make it work - settings.tolerances.absolute_dual_tolerance = settings.minimal_absolute_tolerance; - settings.tolerances.relative_dual_tolerance = 0.0; - settings.tolerances.absolute_primal_tolerance = settings.minimal_absolute_tolerance; - settings.tolerances.relative_primal_tolerance = 0.0; - settings.tolerances.absolute_gap_tolerance = settings.minimal_absolute_tolerance; - settings.tolerances.relative_gap_tolerance = 0.0; - settings.method = cuopt::linear_programming::method_t::PDLP; - - optimization_problem_solution_t solution = solve_lp(&handle_, op_problem, settings); - EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); - EXPECT_FALSE(is_incorrect_objective( - afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); -} - -TEST(pdlp_class, run_double_initial_solution) -{ - const raft::handle_t handle_{}; - - auto path = make_path_absolute("linear_programming/afiro_original.mps"); - cuopt::mps_parser::mps_data_model_t op_problem = - cuopt::mps_parser::parse_mps(path, true); - - std::vector inital_primal_sol(op_problem.get_n_variables()); - std::fill(inital_primal_sol.begin(), inital_primal_sol.end(), 1.0); - op_problem.set_initial_primal_solution(inital_primal_sol.data(), inital_primal_sol.size()); - - auto solver_settings = pdlp_solver_settings_t{}; - solver_settings.method = cuopt::linear_programming::method_t::PDLP; - - optimization_problem_solution_t solution = - solve_lp(&handle_, op_problem, solver_settings); - EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); - EXPECT_FALSE(is_incorrect_objective( - afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); -} - -TEST(pdlp_class, run_iteration_limit) -{ - const raft::handle_t handle_{}; - - auto path = make_path_absolute("linear_programming/afiro_original.mps"); - cuopt::mps_parser::mps_data_model_t op_problem = - cuopt::mps_parser::parse_mps(path, true); - - cuopt::linear_programming::pdlp_solver_settings_t settings = - cuopt::linear_programming::pdlp_solver_settings_t{}; - - settings.iteration_limit = 10; - // To make sure it doesn't return before the iteration limit - settings.set_optimality_tolerance(0); - settings.method = cuopt::linear_programming::method_t::PDLP; - - optimization_problem_solution_t solution = solve_lp(&handle_, op_problem, settings); - EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_ITERATION_LIMIT); - // By default we would return all 0, we now return what we currently have so not all 0 - EXPECT_FALSE(thrust::all_of(handle_.get_thrust_policy(), - solution.get_primal_solution().begin(), - solution.get_primal_solution().end(), - thrust::placeholders::_1 == 0.0)); -} - -TEST(pdlp_class, run_time_limit) -{ - const raft::handle_t handle_{}; - auto path = make_path_absolute("linear_programming/savsched1/savsched1.mps"); - cuopt::mps_parser::mps_data_model_t op_problem = - cuopt::mps_parser::parse_mps(path); - - cuopt::linear_programming::pdlp_solver_settings_t settings = - cuopt::linear_programming::pdlp_solver_settings_t{}; - - constexpr double time_limit_seconds = 2; - settings.time_limit = time_limit_seconds; - // To make sure it doesn't return before the time limit - settings.set_optimality_tolerance(0); - settings.method = cuopt::linear_programming::method_t::PDLP; - - optimization_problem_solution_t solution = solve_lp(&handle_, op_problem, settings); - - EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_TIME_LIMIT); - // By default we would return all 0, we now return what we currently have so not all 0 - EXPECT_FALSE(thrust::all_of(handle_.get_thrust_policy(), - solution.get_primal_solution().begin(), - solution.get_primal_solution().end(), - thrust::placeholders::_1 == 0.0)); - // Check that indeed it didn't run for more than x time - EXPECT_TRUE(solution.get_additional_termination_information().solve_time < - (time_limit_seconds * 5) * 1000); -} - -TEST(pdlp_class, run_sub_mittleman) -{ - std::vector> // Expected objective value - instances{{"graph40-40", -300.0}, - {"ex10", 100.0003411893773}, - {"datt256_lp", 255.9992298290425}, - {"woodlands09", 0.0}, - {"savsched1", 217.4054085795689}, - // {"nug08-3rd", 214.0141488989151}, // TODO: Fix this instance - {"qap15", 1040.999546647414}, - {"scpm1", 413.7787723060584}, - // {"neos3", 27773.54059633068}, // TODO: Fix this instance - {"a2864", -282.9962521965164}}; - - for (const auto& entry : instances) { - const auto& name = entry.first; - const auto expected_objective_value = entry.second; - - auto path = make_path_absolute("linear_programming/" + name + "/" + name + ".mps"); - cuopt::mps_parser::mps_data_model_t op_problem = - cuopt::mps_parser::parse_mps(path); - - // Testing for each solver_mode is ok as it's parsing that is the bottleneck here, not - // solving - auto solver_mode_list = { - cuopt::linear_programming::pdlp_solver_mode_t::Stable3, - cuopt::linear_programming::pdlp_solver_mode_t::Stable2, - cuopt::linear_programming::pdlp_solver_mode_t::Stable1, - cuopt::linear_programming::pdlp_solver_mode_t::Methodical1, - cuopt::linear_programming::pdlp_solver_mode_t::Fast1, - }; - for (auto solver_mode : solver_mode_list) { - auto settings = pdlp_solver_settings_t{}; - settings.pdlp_solver_mode = solver_mode; - settings.dual_postsolve = false; - for (auto [presolve, epsilon] : {std::pair{true, 1e-1}, std::pair{false, 1e-6}}) { - settings.presolve = presolve; - settings.method = cuopt::linear_programming::method_t::PDLP; - const raft::handle_t handle_{}; - optimization_problem_solution_t solution = - solve_lp(&handle_, op_problem, settings); - printf( - "running %s mode %d presolve? %d\n", name.c_str(), (int)solver_mode, settings.presolve); - EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); - EXPECT_FALSE(is_incorrect_objective( - expected_objective_value, - solution.get_additional_termination_information().primal_objective)); - test_objective_sanity(op_problem, - solution.get_primal_solution(), - solution.get_additional_termination_information().primal_objective, - epsilon); - test_constraint_sanity(op_problem, solution, epsilon, presolve); - } - } - } -} +// TEST(pdlp_class, run_double) +// { +// const raft::handle_t handle_{}; + +// auto path = make_path_absolute("linear_programming/afiro_original.mps"); +// cuopt::mps_parser::mps_data_model_t op_problem = +// cuopt::mps_parser::parse_mps(path, true); + +// auto solver_settings = pdlp_solver_settings_t{}; +// solver_settings.method = cuopt::linear_programming::method_t::PDLP; + +// optimization_problem_solution_t solution = +// solve_lp(&handle_, op_problem, solver_settings); +// EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); +// EXPECT_FALSE(is_incorrect_objective( +// afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); +// } + +// TEST(pdlp_class, run_double_very_low_accuracy) +// { +// const raft::handle_t handle_{}; + +// auto path = make_path_absolute("linear_programming/afiro_original.mps"); +// cuopt::mps_parser::mps_data_model_t op_problem = +// cuopt::mps_parser::parse_mps(path, true); + +// cuopt::linear_programming::pdlp_solver_settings_t settings = +// cuopt::linear_programming::pdlp_solver_settings_t{}; +// // With all 0 afiro with return an error +// // Setting absolute tolerance to the minimal value of 1e-12 will make it work +// settings.tolerances.absolute_dual_tolerance = settings.minimal_absolute_tolerance; +// settings.tolerances.relative_dual_tolerance = 0.0; +// settings.tolerances.absolute_primal_tolerance = settings.minimal_absolute_tolerance; +// settings.tolerances.relative_primal_tolerance = 0.0; +// settings.tolerances.absolute_gap_tolerance = settings.minimal_absolute_tolerance; +// settings.tolerances.relative_gap_tolerance = 0.0; +// settings.method = cuopt::linear_programming::method_t::PDLP; + +// optimization_problem_solution_t solution = solve_lp(&handle_, op_problem, +// settings); EXPECT_EQ((int)solution.get_termination_status(), +// CUOPT_TERIMINATION_STATUS_OPTIMAL); EXPECT_FALSE(is_incorrect_objective( +// afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); +// } + +// TEST(pdlp_class, run_double_initial_solution) +// { +// const raft::handle_t handle_{}; + +// auto path = make_path_absolute("linear_programming/afiro_original.mps"); +// cuopt::mps_parser::mps_data_model_t op_problem = +// cuopt::mps_parser::parse_mps(path, true); + +// std::vector inital_primal_sol(op_problem.get_n_variables()); +// std::fill(inital_primal_sol.begin(), inital_primal_sol.end(), 1.0); +// op_problem.set_initial_primal_solution(inital_primal_sol.data(), inital_primal_sol.size()); + +// auto solver_settings = pdlp_solver_settings_t{}; +// solver_settings.method = cuopt::linear_programming::method_t::PDLP; + +// optimization_problem_solution_t solution = +// solve_lp(&handle_, op_problem, solver_settings); +// EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); +// EXPECT_FALSE(is_incorrect_objective( +// afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); +// } + +// TEST(pdlp_class, run_iteration_limit) +// { +// const raft::handle_t handle_{}; + +// auto path = make_path_absolute("linear_programming/afiro_original.mps"); +// cuopt::mps_parser::mps_data_model_t op_problem = +// cuopt::mps_parser::parse_mps(path, true); + +// cuopt::linear_programming::pdlp_solver_settings_t settings = +// cuopt::linear_programming::pdlp_solver_settings_t{}; + +// settings.iteration_limit = 10; +// // To make sure it doesn't return before the iteration limit +// settings.set_optimality_tolerance(0); +// settings.method = cuopt::linear_programming::method_t::PDLP; + +// optimization_problem_solution_t solution = solve_lp(&handle_, op_problem, +// settings); EXPECT_EQ((int)solution.get_termination_status(), +// CUOPT_TERIMINATION_STATUS_ITERATION_LIMIT); +// // By default we would return all 0, we now return what we currently have so not all 0 +// EXPECT_FALSE(thrust::all_of(handle_.get_thrust_policy(), +// solution.get_primal_solution().begin(), +// solution.get_primal_solution().end(), +// thrust::placeholders::_1 == 0.0)); +// } + +// TEST(pdlp_class, run_time_limit) +// { +// const raft::handle_t handle_{}; +// auto path = make_path_absolute("linear_programming/savsched1/savsched1.mps"); +// cuopt::mps_parser::mps_data_model_t op_problem = +// cuopt::mps_parser::parse_mps(path); + +// cuopt::linear_programming::pdlp_solver_settings_t settings = +// cuopt::linear_programming::pdlp_solver_settings_t{}; + +// constexpr double time_limit_seconds = 2; +// settings.time_limit = time_limit_seconds; +// // To make sure it doesn't return before the time limit +// settings.set_optimality_tolerance(0); +// settings.method = cuopt::linear_programming::method_t::PDLP; + +// optimization_problem_solution_t solution = solve_lp(&handle_, op_problem, +// settings); + +// EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_TIME_LIMIT); +// // By default we would return all 0, we now return what we currently have so not all 0 +// EXPECT_FALSE(thrust::all_of(handle_.get_thrust_policy(), +// solution.get_primal_solution().begin(), +// solution.get_primal_solution().end(), +// thrust::placeholders::_1 == 0.0)); +// // Check that indeed it didn't run for more than x time +// EXPECT_TRUE(solution.get_additional_termination_information().solve_time < +// (time_limit_seconds * 5) * 1000); +// } + +// TEST(pdlp_class, run_sub_mittleman) +// { +// std::vector> // Expected objective value +// instances{{"graph40-40", -300.0}, +// {"ex10", 100.0003411893773}, +// {"datt256_lp", 255.9992298290425}, +// {"woodlands09", 0.0}, +// {"savsched1", 217.4054085795689}, +// // {"nug08-3rd", 214.0141488989151}, // TODO: Fix this instance +// {"qap15", 1040.999546647414}, +// {"scpm1", 413.7787723060584}, +// // {"neos3", 27773.54059633068}, // TODO: Fix this instance +// {"a2864", -282.9962521965164}}; + +// for (const auto& entry : instances) { +// const auto& name = entry.first; +// const auto expected_objective_value = entry.second; + +// auto path = make_path_absolute("linear_programming/" + name + "/" + name + ".mps"); +// cuopt::mps_parser::mps_data_model_t op_problem = +// cuopt::mps_parser::parse_mps(path); + +// // Testing for each solver_mode is ok as it's parsing that is the bottleneck here, not +// // solving +// auto solver_mode_list = { +// cuopt::linear_programming::pdlp_solver_mode_t::Stable3, +// cuopt::linear_programming::pdlp_solver_mode_t::Stable2, +// cuopt::linear_programming::pdlp_solver_mode_t::Stable1, +// cuopt::linear_programming::pdlp_solver_mode_t::Methodical1, +// cuopt::linear_programming::pdlp_solver_mode_t::Fast1, +// }; +// for (auto solver_mode : solver_mode_list) { +// auto settings = pdlp_solver_settings_t{}; +// settings.pdlp_solver_mode = solver_mode; +// settings.dual_postsolve = false; +// for (auto [presolve, epsilon] : {std::pair{true, 1e-1}, std::pair{false, 1e-6}}) { +// settings.presolve = presolve; +// settings.method = cuopt::linear_programming::method_t::PDLP; +// const raft::handle_t handle_{}; +// optimization_problem_solution_t solution = +// solve_lp(&handle_, op_problem, settings); +// printf( +// "running %s mode %d presolve? %d\n", name.c_str(), (int)solver_mode, +// settings.presolve); +// EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); +// EXPECT_FALSE(is_incorrect_objective( +// expected_objective_value, +// solution.get_additional_termination_information().primal_objective)); +// test_objective_sanity(op_problem, +// solution.get_primal_solution(), +// solution.get_additional_termination_information().primal_objective, +// epsilon); +// test_constraint_sanity(op_problem, solution, epsilon, presolve); +// } +// } +// } +// } constexpr double initial_step_size_afiro = 1.4893; constexpr double initial_primal_weight_afiro = 0.0141652; constexpr double factor_tolerance = 1e-4f; -// Should be added to google test +// // Should be added to google test #define EXPECT_NOT_NEAR(val1, val2, abs_error) \ EXPECT_FALSE((std::abs((val1) - (val2)) <= (abs_error))) -TEST(pdlp_class, initial_solution_test) -{ - const raft::handle_t handle_{}; - - auto path = make_path_absolute("linear_programming/afiro_original.mps"); - cuopt::mps_parser::mps_data_model_t mps_data_model = - cuopt::mps_parser::parse_mps(path); - - auto op_problem = cuopt::linear_programming::mps_data_model_to_optimization_problem( - &handle_, mps_data_model); - cuopt::linear_programming::detail::problem_t problem(op_problem); - - auto solver_settings = pdlp_solver_settings_t{}; - // We are just testing initial scaling on initial solution scheme so we don't care about solver - solver_settings.iteration_limit = 0; - solver_settings.method = cuopt::linear_programming::method_t::PDLP; - // Empty call solve to set the parameters and init the handler since calling pdlp object directly - // doesn't - solver_settings.pdlp_solver_mode = cuopt::linear_programming::pdlp_solver_mode_t::Methodical1; - solve_lp(op_problem, solver_settings); - EXPECT_EQ(cuopt::linear_programming::pdlp_hyper_params::initial_step_size_scaling, 1); - EXPECT_EQ(cuopt::linear_programming::pdlp_hyper_params::default_l_inf_ruiz_iterations, 5); - EXPECT_TRUE(cuopt::linear_programming::pdlp_hyper_params::do_pock_chambolle_scaling); - EXPECT_TRUE(cuopt::linear_programming::pdlp_hyper_params::do_ruiz_scaling); - EXPECT_EQ(cuopt::linear_programming::pdlp_hyper_params::default_alpha_pock_chambolle_rescaling, - 1.0); - - EXPECT_FALSE(cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution); - EXPECT_FALSE( - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution); - - { - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - solver.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - } - - // First add an initial primal then dual, then both, which shouldn't influence the values as the - // scale on initial option is not toggled - { - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 1); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - solver.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - } - { - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_dual(op_problem.get_n_constraints(), 1); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - } - { - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 1); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - std::vector initial_dual(op_problem.get_n_constraints(), 1); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - } - - // Toggle the scale on initial solution while not providing should yield the same - { - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; - solver.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; - } - { - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = true; - solver.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = false; - } - { - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = true; - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; - solver.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = false; - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; - } - - // Asking for initial scaling on step size with initial solution being only primal or only dual - // should not break but not modify the step size - { - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 1); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - solver.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; - } - { - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_dual(op_problem.get_n_constraints(), 1); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; - } - - // Asking for initial scaling on primal weight with initial solution being only primal or only - // dual should *not* break but the primal weight should not change - { - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 1); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - solver.run_solver(pdlp_timer); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = false; - } - { - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_dual(op_problem.get_n_constraints(), 1); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = false; - } - - // All 0 solution when given an initial primal and dual with scale on the step size should not - // break but not change primal weight and step size - { - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 0); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - std::vector initial_dual(op_problem.get_n_constraints(), 0); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; - } - - // All 0 solution when given an initial primal and/or dual with scale on the primal weight is - // *not* an error but should not change primal weight and step size - { - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 0); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - solver.run_solver(pdlp_timer); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = false; - } - { - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_dual(op_problem.get_n_constraints(), 0); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = false; - } - { - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 0); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - std::vector initial_dual(op_problem.get_n_constraints(), 0); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = false; - } - - // A non-all-0 vector for both initial primal and dual set should trigger a modification in primal - // weight and step size - { - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 1); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - std::vector initial_dual(op_problem.get_n_constraints(), 1); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NOT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = false; - } - { - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 1); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - std::vector initial_dual(op_problem.get_n_constraints(), 1); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - EXPECT_NOT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; - } - { - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = true; - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 1); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - std::vector initial_dual(op_problem.get_n_constraints(), 1); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - EXPECT_NOT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); - EXPECT_NOT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = false; - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; - } -} - -TEST(pdlp_class, initial_primal_weight_step_size_test) -{ - const raft::handle_t handle_{}; - - auto path = make_path_absolute("linear_programming/afiro_original.mps"); - cuopt::mps_parser::mps_data_model_t mps_data_model = - cuopt::mps_parser::parse_mps(path); - - auto op_problem = cuopt::linear_programming::mps_data_model_to_optimization_problem( - &handle_, mps_data_model); - cuopt::linear_programming::detail::problem_t problem(op_problem); - - auto solver_settings = pdlp_solver_settings_t{}; - // We are just testing initial scaling on initial solution scheme so we don't care about solver - solver_settings.iteration_limit = 0; - solver_settings.method = cuopt::linear_programming::method_t::PDLP; - // Select the default/legacy solver with no action upon the initial scaling on initial solution - solver_settings.pdlp_solver_mode = cuopt::linear_programming::pdlp_solver_mode_t::Methodical1; - EXPECT_FALSE(cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution); - EXPECT_FALSE( - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution); - - // Check setting an initial primal weight and step size - { - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - constexpr double test_initial_step_size = 1.0; - constexpr double test_initial_primal_weight = 2.0; - solver.set_initial_primal_weight(test_initial_primal_weight); - solver.set_initial_step_size(test_initial_step_size); - solver.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_EQ(test_initial_step_size, solver.get_step_size_h()); - EXPECT_EQ(test_initial_primal_weight, solver.get_primal_weight_h()); - } - - // Check that after setting an initial step size and primal weight, the computed one when adding - // an initial primal / dual is indeed different - { - // Launching without an inital step size / primal weight and query the value - cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = true; - cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - auto pdlp_timer = timer_t(solver_settings.time_limit); - std::vector initial_primal(op_problem.get_n_variables(), 1); - auto d_initial_primal = device_copy(initial_primal, handle_.get_stream()); - solver.set_initial_primal_solution(d_initial_primal); - std::vector initial_dual(op_problem.get_n_constraints(), 1); - auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); - solver.set_initial_dual_solution(d_initial_dual); - solver.run_solver(pdlp_timer); - const double previous_step_size = solver.get_step_size_h(); - const double previous_primal_weight = solver.get_primal_weight_h(); - - // Start again but with an initial and check the impact - cuopt::linear_programming::detail::pdlp_solver_t solver2(problem, solver_settings); - pdlp_timer = timer_t(solver_settings.time_limit); - constexpr double test_initial_step_size = 1.0; - constexpr double test_initial_primal_weight = 2.0; - solver2.set_initial_primal_weight(test_initial_primal_weight); - solver2.set_initial_step_size(test_initial_step_size); - solver2.set_initial_primal_solution(d_initial_primal); - solver2.set_initial_dual_solution(d_initial_dual); - solver2.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - const double sovler2_step_size = solver2.get_step_size_h(); - const double sovler2_primal_weight = solver2.get_primal_weight_h(); - EXPECT_NOT_NEAR(previous_step_size, sovler2_step_size, factor_tolerance); - EXPECT_NOT_NEAR(previous_primal_weight, sovler2_primal_weight, factor_tolerance); - - // Again but with an initial k which should change the step size only, not the primal weight - cuopt::linear_programming::detail::pdlp_solver_t solver3(problem, solver_settings); - pdlp_timer = timer_t(solver_settings.time_limit); - solver3.set_initial_primal_weight(test_initial_primal_weight); - solver3.set_initial_step_size(test_initial_step_size); - solver3.set_initial_primal_solution(d_initial_primal); - solver3.set_initial_k(10000); - solver3.set_initial_dual_solution(d_initial_dual); - solver3.set_initial_dual_solution(d_initial_dual); - solver3.run_solver(pdlp_timer); - RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); - EXPECT_NOT_NEAR(sovler2_step_size, solver3.get_step_size_h(), factor_tolerance); - EXPECT_NEAR(sovler2_primal_weight, solver3.get_primal_weight_h(), factor_tolerance); - } -} - -TEST(pdlp_class, initial_rhs_and_c) -{ - const raft::handle_t handle_{}; - - auto path = make_path_absolute("linear_programming/afiro_original.mps"); - cuopt::mps_parser::mps_data_model_t mps_data_model = - cuopt::mps_parser::parse_mps(path); - - auto op_problem = cuopt::linear_programming::mps_data_model_to_optimization_problem( - &handle_, mps_data_model); - cuopt::linear_programming::detail::problem_t problem(op_problem); - - cuopt::linear_programming::detail::pdlp_solver_t solver(problem); - constexpr double test_initial_primal_factor = 1.0; - constexpr double test_initial_dual_factor = 2.0; - solver.set_relative_dual_tolerance_factor(test_initial_dual_factor); - solver.set_relative_primal_tolerance_factor(test_initial_primal_factor); - - EXPECT_EQ(solver.get_relative_dual_tolerance_factor(), test_initial_dual_factor); - EXPECT_EQ(solver.get_relative_primal_tolerance_factor(), test_initial_primal_factor); -} - -TEST(pdlp_class, per_constraint_test) -{ - /* - * Define the following LP: - * x1=0.01 <= 0 - * x2=0.01 <= 0 - * x3=0.1 <= 0 - * - * With a tol of 0.1 per constraint will pass but the L2 version will not as L2 of primal residual - * will be 0.1009 - */ - raft::handle_t handle; - auto op_problem = optimization_problem_t(&handle); - - std::vector A_host = {1.0, 1.0, 1.0}; - std::vector indices_host = {0, 1, 2}; - std::vector offset_host = {0, 1, 2, 3}; - std::vector b_host = {0.0, 0.0, 0.0}; - std::vector h_initial_primal = {0.02, 0.03, 0.1}; - rmm::device_uvector d_initial_primal(3, handle.get_stream()); - raft::copy( - d_initial_primal.data(), h_initial_primal.data(), h_initial_primal.size(), handle.get_stream()); - - op_problem.set_csr_constraint_matrix(A_host.data(), - A_host.size(), - indices_host.data(), - indices_host.size(), - offset_host.data(), - offset_host.size()); - op_problem.set_constraint_lower_bounds(b_host.data(), b_host.size()); - op_problem.set_constraint_upper_bounds(b_host.data(), b_host.size()); - op_problem.set_objective_coefficients(b_host.data(), b_host.size()); - - auto problem = cuopt::linear_programming::detail::problem_t(op_problem); - - pdlp_solver_settings_t solver_settings; - solver_settings.tolerances.relative_primal_tolerance = 0; // Shouldn't matter - solver_settings.tolerances.absolute_primal_tolerance = 0.1; - solver_settings.tolerances.relative_dual_tolerance = 0; // Shoudln't matter - solver_settings.tolerances.absolute_dual_tolerance = 0.1; - solver_settings.method = cuopt::linear_programming::method_t::PDLP; - - // First solve without the per constraint and it should break - { - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - - raft::copy(solver.pdhg_solver_.get_primal_solution().data(), - d_initial_primal.data(), - d_initial_primal.size(), - handle.get_stream()); - - auto& current_termination_strategy = solver.get_current_termination_strategy(); - pdlp_termination_status_t termination_average = - current_termination_strategy.evaluate_termination_criteria( - solver.pdhg_solver_, - d_initial_primal, - d_initial_primal, - solver.pdhg_solver_.get_dual_slack(), - problem.combined_bounds, - problem.objective_coefficients); - - EXPECT_TRUE(termination_average != pdlp_termination_status_t::Optimal); - } - { - solver_settings.per_constraint_residual = true; - cuopt::linear_programming::detail::pdlp_solver_t solver(problem, solver_settings); - - raft::copy(solver.pdhg_solver_.get_primal_solution().data(), - d_initial_primal.data(), - d_initial_primal.size(), - handle.get_stream()); - - auto& current_termination_strategy = solver.get_current_termination_strategy(); - pdlp_termination_status_t termination_average = - current_termination_strategy.evaluate_termination_criteria( - solver.pdhg_solver_, - d_initial_primal, - d_initial_primal, - solver.pdhg_solver_.get_dual_slack(), - problem.combined_bounds, - problem.objective_coefficients); - EXPECT_EQ(current_termination_strategy.get_convergence_information() - .get_relative_linf_primal_residual() - .value(handle.get_stream()), - 0.1); - } -} - -TEST(pdlp_class, best_primal_so_far_iteration) -{ - GTEST_SKIP() << "Skipping test: best_primal_so_far_iteration. Enable when ready to run."; - const raft::handle_t handle1{}; - const raft::handle_t handle2{}; - - auto path = make_path_absolute("linear_programming/ns1687037/ns1687037.mps"); - auto solver_settings = pdlp_solver_settings_t{}; - solver_settings.iteration_limit = 3000; - solver_settings.per_constraint_residual = true; - solver_settings.method = cuopt::linear_programming::method_t::PDLP; - - cuopt::mps_parser::mps_data_model_t op_problem1 = - cuopt::mps_parser::parse_mps(path); - cuopt::mps_parser::mps_data_model_t op_problem2 = - cuopt::mps_parser::parse_mps(path); - - optimization_problem_solution_t solution1 = - solve_lp(&handle1, op_problem1, solver_settings); - RAFT_CUDA_TRY(cudaDeviceSynchronize()); - solver_settings.save_best_primal_so_far = true; - optimization_problem_solution_t solution2 = - solve_lp(&handle2, op_problem2, solver_settings); - RAFT_CUDA_TRY(cudaDeviceSynchronize()); - - EXPECT_TRUE(solution2.get_additional_termination_information().l2_primal_residual < - solution1.get_additional_termination_information().l2_primal_residual); -} - -TEST(pdlp_class, best_primal_so_far_time) -{ - GTEST_SKIP() << "Skipping test: best_primal_so_far_time. Enable when ready to run."; - const raft::handle_t handle1{}; - const raft::handle_t handle2{}; - - auto path = make_path_absolute("linear_programming/ns1687037/ns1687037.mps"); - auto solver_settings = pdlp_solver_settings_t{}; - solver_settings.time_limit = 2; - solver_settings.per_constraint_residual = true; - solver_settings.pdlp_solver_mode = cuopt::linear_programming::pdlp_solver_mode_t::Stable1; - solver_settings.method = cuopt::linear_programming::method_t::PDLP; - - cuopt::mps_parser::mps_data_model_t op_problem1 = - cuopt::mps_parser::parse_mps(path); - cuopt::mps_parser::mps_data_model_t op_problem2 = - cuopt::mps_parser::parse_mps(path); - - optimization_problem_solution_t solution1 = - solve_lp(&handle1, op_problem1, solver_settings); - RAFT_CUDA_TRY(cudaDeviceSynchronize()); - solver_settings.save_best_primal_so_far = true; - optimization_problem_solution_t solution2 = - solve_lp(&handle2, op_problem2, solver_settings); - RAFT_CUDA_TRY(cudaDeviceSynchronize()); - - EXPECT_TRUE(solution2.get_additional_termination_information().l2_primal_residual < - solution1.get_additional_termination_information().l2_primal_residual); -} - -TEST(pdlp_class, first_primal_feasible) -{ - GTEST_SKIP() << "Skipping test: first_primal_feasible. Enable when ready to run."; - const raft::handle_t handle1{}; - const raft::handle_t handle2{}; - - auto path = make_path_absolute("linear_programming/ns1687037/ns1687037.mps"); - auto solver_settings = pdlp_solver_settings_t{}; - solver_settings.iteration_limit = 1000; - solver_settings.per_constraint_residual = true; - solver_settings.set_optimality_tolerance(1e-2); - solver_settings.method = cuopt::linear_programming::method_t::PDLP; - - cuopt::mps_parser::mps_data_model_t op_problem1 = - cuopt::mps_parser::parse_mps(path); - cuopt::mps_parser::mps_data_model_t op_problem2 = - cuopt::mps_parser::parse_mps(path); - - optimization_problem_solution_t solution1 = - solve_lp(&handle1, op_problem1, solver_settings); - RAFT_CUDA_TRY(cudaDeviceSynchronize()); - solver_settings.first_primal_feasible = true; - optimization_problem_solution_t solution2 = - solve_lp(&handle2, op_problem2, solver_settings); - RAFT_CUDA_TRY(cudaDeviceSynchronize()); - - EXPECT_EQ(solution1.get_termination_status(), pdlp_termination_status_t::IterationLimit); - EXPECT_EQ(solution2.get_termination_status(), pdlp_termination_status_t::PrimalFeasible); -} - -TEST(pdlp_class, warm_start) -{ - std::vector instance_names{"graph40-40", - "ex10", - "datt256_lp", - "woodlands09", - "savsched1", - // "nug08-3rd", // TODO: Fix this instance - "qap15", - "scpm1", - // "neos3", // TODO: Fix this instance - "a2864"}; - for (auto instance_name : instance_names) { - const raft::handle_t handle{}; - - auto path = - make_path_absolute("linear_programming/" + instance_name + "/" + instance_name + ".mps"); - auto solver_settings = pdlp_solver_settings_t{}; - solver_settings.pdlp_solver_mode = cuopt::linear_programming::pdlp_solver_mode_t::Stable2; - solver_settings.set_optimality_tolerance(1e-2); - solver_settings.detect_infeasibility = false; - solver_settings.method = cuopt::linear_programming::method_t::PDLP; - - cuopt::mps_parser::mps_data_model_t mps_data_model = - cuopt::mps_parser::parse_mps(path); - auto op_problem1 = - cuopt::linear_programming::mps_data_model_to_optimization_problem( - &handle, mps_data_model); - - // Solving from scratch until 1e-2 - optimization_problem_solution_t solution1 = solve_lp(op_problem1, solver_settings); - - // Solving until 1e-1 to use the result as a warm start - solver_settings.set_optimality_tolerance(1e-1); - auto op_problem2 = - cuopt::linear_programming::mps_data_model_to_optimization_problem( - &handle, mps_data_model); - optimization_problem_solution_t solution2 = solve_lp(op_problem2, solver_settings); - - // Solving until 1e-2 using the previous state as a warm start - solver_settings.set_optimality_tolerance(1e-2); - auto op_problem3 = - cuopt::linear_programming::mps_data_model_to_optimization_problem( - &handle, mps_data_model); - solver_settings.set_pdlp_warm_start_data(solution2.get_pdlp_warm_start_data()); - optimization_problem_solution_t solution3 = solve_lp(op_problem3, solver_settings); - - EXPECT_EQ(solution1.get_additional_termination_information().number_of_steps_taken, - solution3.get_additional_termination_information().number_of_steps_taken + - solution2.get_additional_termination_information().number_of_steps_taken); - } -} - -TEST(pdlp_class, dual_postsolve_size) -{ - const raft::handle_t handle_{}; - - auto path = make_path_absolute("linear_programming/afiro_original.mps"); - cuopt::mps_parser::mps_data_model_t op_problem = - cuopt::mps_parser::parse_mps(path, true); - - auto solver_settings = pdlp_solver_settings_t{}; - solver_settings.method = cuopt::linear_programming::method_t::PDLP; - solver_settings.presolve = true; - - { - solver_settings.dual_postsolve = true; - optimization_problem_solution_t solution = - solve_lp(&handle_, op_problem, solver_settings); - EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); - EXPECT_EQ(solution.get_dual_solution().size(), op_problem.get_n_constraints()); - } - - { - solver_settings.dual_postsolve = false; - optimization_problem_solution_t solution = - solve_lp(&handle_, op_problem, solver_settings); - EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); - EXPECT_EQ(solution.get_dual_solution().size(), 0); - } -} - -TEST(dual_simplex, afiro) -{ - cuopt::linear_programming::pdlp_solver_settings_t settings = - cuopt::linear_programming::pdlp_solver_settings_t{}; - settings.method = cuopt::linear_programming::method_t::DualSimplex; - - const raft::handle_t handle_{}; - - auto path = make_path_absolute("linear_programming/afiro_original.mps"); - cuopt::mps_parser::mps_data_model_t op_problem = - cuopt::mps_parser::parse_mps(path, true); - - optimization_problem_solution_t solution = solve_lp(&handle_, op_problem, settings); - EXPECT_EQ(solution.get_termination_status(), pdlp_termination_status_t::Optimal); - EXPECT_FALSE(is_incorrect_objective( - afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); -} +// TEST(pdlp_class, initial_solution_test) +// { +// const raft::handle_t handle_{}; + +// auto path = make_path_absolute("linear_programming/afiro_original.mps"); +// cuopt::mps_parser::mps_data_model_t mps_data_model = +// cuopt::mps_parser::parse_mps(path); + +// auto op_problem = cuopt::linear_programming::mps_data_model_to_optimization_problem( +// &handle_, mps_data_model); +// cuopt::linear_programming::detail::problem_t problem(op_problem); + +// auto solver_settings = pdlp_solver_settings_t{}; +// // We are just testing initial scaling on initial solution scheme so we don't care about solver +// solver_settings.iteration_limit = 0; +// solver_settings.method = cuopt::linear_programming::method_t::PDLP; +// // Empty call solve to set the parameters and init the handler since calling pdlp object +// directly +// // doesn't +// solver_settings.pdlp_solver_mode = cuopt::linear_programming::pdlp_solver_mode_t::Methodical1; +// solve_lp(op_problem, solver_settings); +// EXPECT_EQ(cuopt::linear_programming::pdlp_hyper_params::initial_step_size_scaling, 1); +// EXPECT_EQ(cuopt::linear_programming::pdlp_hyper_params::default_l_inf_ruiz_iterations, 5); +// EXPECT_TRUE(cuopt::linear_programming::pdlp_hyper_params::do_pock_chambolle_scaling); +// EXPECT_TRUE(cuopt::linear_programming::pdlp_hyper_params::do_ruiz_scaling); +// EXPECT_EQ(cuopt::linear_programming::pdlp_hyper_params::default_alpha_pock_chambolle_rescaling, +// 1.0); + +// EXPECT_FALSE(cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution); +// EXPECT_FALSE( +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution); + +// { +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); +// solver.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// } + +// // First add an initial primal then dual, then both, which shouldn't influence the values as +// the +// // scale on initial option is not toggled +// { +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 1); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// solver.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// } +// { +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_dual(op_problem.get_n_constraints(), 1); auto d_initial_dual = +// device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// } +// { +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 1); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// std::vector initial_dual(op_problem.get_n_constraints(), 1); +// auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// } + +// // Toggle the scale on initial solution while not providing should yield the same +// { +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); +// cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; +// solver.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; +// } +// { +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// true; solver.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// false; +// } +// { +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// true; cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = +// true; solver.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// false; cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = +// false; +// } + +// // Asking for initial scaling on step size with initial solution being only primal or only dual +// // should not break but not modify the step size +// { +// cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 1); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// solver.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; +// } +// { +// cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_dual(op_problem.get_n_constraints(), 1); auto d_initial_dual = +// device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; +// } + +// // Asking for initial scaling on primal weight with initial solution being only primal or only +// // dual should *not* break but the primal weight should not change +// { +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// true; cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 1); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// solver.run_solver(pdlp_timer); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// false; +// } +// { +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// true; cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_dual(op_problem.get_n_constraints(), 1); auto d_initial_dual = +// device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// false; +// } + +// // All 0 solution when given an initial primal and dual with scale on the step size should not +// // break but not change primal weight and step size +// { +// cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 0); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// std::vector initial_dual(op_problem.get_n_constraints(), 0); +// auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; +// } + +// // All 0 solution when given an initial primal and/or dual with scale on the primal weight is +// // *not* an error but should not change primal weight and step size +// { +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// true; cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 0); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// solver.run_solver(pdlp_timer); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// false; +// } +// { +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// true; cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_dual(op_problem.get_n_constraints(), 0); auto d_initial_dual = +// device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// false; +// } +// { +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// true; cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 0); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// std::vector initial_dual(op_problem.get_n_constraints(), 0); +// auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// false; +// } + +// // A non-all-0 vector for both initial primal and dual set should trigger a modification in +// primal +// // weight and step size +// { +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// true; cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 1); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// std::vector initial_dual(op_problem.get_n_constraints(), 1); +// auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// EXPECT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NOT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// false; +// } +// { +// cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = true; +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 1); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// std::vector initial_dual(op_problem.get_n_constraints(), 1); +// auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// EXPECT_NOT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = false; +// } +// { +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// true; cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = +// true; cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 1); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// std::vector initial_dual(op_problem.get_n_constraints(), 1); +// auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// EXPECT_NOT_NEAR(initial_step_size_afiro, solver.get_step_size_h(), factor_tolerance); +// EXPECT_NOT_NEAR(initial_primal_weight_afiro, solver.get_primal_weight_h(), factor_tolerance); +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// false; cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = +// false; +// } +// } + +// TEST(pdlp_class, initial_primal_weight_step_size_test) +// { +// const raft::handle_t handle_{}; + +// auto path = make_path_absolute("linear_programming/afiro_original.mps"); +// cuopt::mps_parser::mps_data_model_t mps_data_model = +// cuopt::mps_parser::parse_mps(path); + +// auto op_problem = cuopt::linear_programming::mps_data_model_to_optimization_problem( +// &handle_, mps_data_model); +// cuopt::linear_programming::detail::problem_t problem(op_problem); + +// auto solver_settings = pdlp_solver_settings_t{}; +// // We are just testing initial scaling on initial solution scheme so we don't care about solver +// solver_settings.iteration_limit = 0; +// solver_settings.method = cuopt::linear_programming::method_t::PDLP; +// // Select the default/legacy solver with no action upon the initial scaling on initial solution +// solver_settings.pdlp_solver_mode = cuopt::linear_programming::pdlp_solver_mode_t::Methodical1; +// EXPECT_FALSE(cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution); +// EXPECT_FALSE( +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution); + +// // Check setting an initial primal weight and step size +// { +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = +// timer_t(solver_settings.time_limit); constexpr double test_initial_step_size = 1.0; +// constexpr double test_initial_primal_weight = 2.0; +// solver.set_initial_primal_weight(test_initial_primal_weight); +// solver.set_initial_step_size(test_initial_step_size); +// solver.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_EQ(test_initial_step_size, solver.get_step_size_h()); +// EXPECT_EQ(test_initial_primal_weight, solver.get_primal_weight_h()); +// } + +// // Check that after setting an initial step size and primal weight, the computed one when +// adding +// // an initial primal / dual is indeed different +// { +// // Launching without an inital step size / primal weight and query the value +// cuopt::linear_programming::pdlp_hyper_params::update_primal_weight_on_initial_solution = +// true; cuopt::linear_programming::pdlp_hyper_params::update_step_size_on_initial_solution = +// true; cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); auto pdlp_timer = timer_t(solver_settings.time_limit); std::vector +// initial_primal(op_problem.get_n_variables(), 1); auto d_initial_primal = +// device_copy(initial_primal, handle_.get_stream()); +// solver.set_initial_primal_solution(d_initial_primal); +// std::vector initial_dual(op_problem.get_n_constraints(), 1); +// auto d_initial_dual = device_copy(initial_dual, handle_.get_stream()); +// solver.set_initial_dual_solution(d_initial_dual); +// solver.run_solver(pdlp_timer); +// const double previous_step_size = solver.get_step_size_h(); +// const double previous_primal_weight = solver.get_primal_weight_h(); + +// // Start again but with an initial and check the impact +// cuopt::linear_programming::detail::pdlp_solver_t solver2(problem, +// solver_settings); pdlp_timer = +// timer_t(solver_settings.time_limit); constexpr double test_initial_step_size = 1.0; +// constexpr double test_initial_primal_weight = 2.0; +// solver2.set_initial_primal_weight(test_initial_primal_weight); +// solver2.set_initial_step_size(test_initial_step_size); +// solver2.set_initial_primal_solution(d_initial_primal); +// solver2.set_initial_dual_solution(d_initial_dual); +// solver2.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// const double sovler2_step_size = solver2.get_step_size_h(); +// const double sovler2_primal_weight = solver2.get_primal_weight_h(); +// EXPECT_NOT_NEAR(previous_step_size, sovler2_step_size, factor_tolerance); +// EXPECT_NOT_NEAR(previous_primal_weight, sovler2_primal_weight, factor_tolerance); + +// // Again but with an initial k which should change the step size only, not the primal weight +// cuopt::linear_programming::detail::pdlp_solver_t solver3(problem, +// solver_settings); pdlp_timer = timer_t(solver_settings.time_limit); +// solver3.set_initial_primal_weight(test_initial_primal_weight); +// solver3.set_initial_step_size(test_initial_step_size); +// solver3.set_initial_primal_solution(d_initial_primal); +// solver3.set_initial_k(10000); +// solver3.set_initial_dual_solution(d_initial_dual); +// solver3.set_initial_dual_solution(d_initial_dual); +// solver3.run_solver(pdlp_timer); +// RAFT_CUDA_TRY(cudaStreamSynchronize(handle_.get_stream())); +// EXPECT_NOT_NEAR(sovler2_step_size, solver3.get_step_size_h(), factor_tolerance); +// EXPECT_NEAR(sovler2_primal_weight, solver3.get_primal_weight_h(), factor_tolerance); +// } +// } + +// TEST(pdlp_class, initial_rhs_and_c) +// { +// const raft::handle_t handle_{}; + +// auto path = make_path_absolute("linear_programming/afiro_original.mps"); +// cuopt::mps_parser::mps_data_model_t mps_data_model = +// cuopt::mps_parser::parse_mps(path); + +// auto op_problem = cuopt::linear_programming::mps_data_model_to_optimization_problem( +// &handle_, mps_data_model); +// cuopt::linear_programming::detail::problem_t problem(op_problem); + +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem); +// constexpr double test_initial_primal_factor = 1.0; +// constexpr double test_initial_dual_factor = 2.0; +// solver.set_relative_dual_tolerance_factor(test_initial_dual_factor); +// solver.set_relative_primal_tolerance_factor(test_initial_primal_factor); + +// EXPECT_EQ(solver.get_relative_dual_tolerance_factor(), test_initial_dual_factor); +// EXPECT_EQ(solver.get_relative_primal_tolerance_factor(), test_initial_primal_factor); +// } + +// TEST(pdlp_class, per_constraint_test) +// { +// /* +// * Define the following LP: +// * x1=0.01 <= 0 +// * x2=0.01 <= 0 +// * x3=0.1 <= 0 +// * +// * With a tol of 0.1 per constraint will pass but the L2 version will not as L2 of primal +// residual +// * will be 0.1009 +// */ +// raft::handle_t handle; +// auto op_problem = optimization_problem_t(&handle); + +// std::vector A_host = {1.0, 1.0, 1.0}; +// std::vector indices_host = {0, 1, 2}; +// std::vector offset_host = {0, 1, 2, 3}; +// std::vector b_host = {0.0, 0.0, 0.0}; +// std::vector h_initial_primal = {0.02, 0.03, 0.1}; +// rmm::device_uvector d_initial_primal(3, handle.get_stream()); +// raft::copy( +// d_initial_primal.data(), h_initial_primal.data(), h_initial_primal.size(), +// handle.get_stream()); + +// op_problem.set_csr_constraint_matrix(A_host.data(), +// A_host.size(), +// indices_host.data(), +// indices_host.size(), +// offset_host.data(), +// offset_host.size()); +// op_problem.set_constraint_lower_bounds(b_host.data(), b_host.size()); +// op_problem.set_constraint_upper_bounds(b_host.data(), b_host.size()); +// op_problem.set_objective_coefficients(b_host.data(), b_host.size()); + +// auto problem = cuopt::linear_programming::detail::problem_t(op_problem); + +// pdlp_solver_settings_t solver_settings; +// solver_settings.tolerances.relative_primal_tolerance = 0; // Shouldn't matter +// solver_settings.tolerances.absolute_primal_tolerance = 0.1; +// solver_settings.tolerances.relative_dual_tolerance = 0; // Shoudln't matter +// solver_settings.tolerances.absolute_dual_tolerance = 0.1; +// solver_settings.method = +// cuopt::linear_programming::method_t::PDLP; + +// // First solve without the per constraint and it should break +// { +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); + +// raft::copy(solver.pdhg_solver_.get_primal_solution().data(), +// d_initial_primal.data(), +// d_initial_primal.size(), +// handle.get_stream()); + +// auto& current_termination_strategy = solver.get_current_termination_strategy(); +// pdlp_termination_status_t termination_average = +// current_termination_strategy.evaluate_termination_criteria( +// solver.pdhg_solver_, +// d_initial_primal, +// d_initial_primal, +// solver.pdhg_solver_.get_dual_slack(), +// problem.combined_bounds, +// problem.objective_coefficients); + +// EXPECT_TRUE(termination_average != pdlp_termination_status_t::Optimal); +// } +// { +// solver_settings.per_constraint_residual = true; +// cuopt::linear_programming::detail::pdlp_solver_t solver(problem, +// solver_settings); + +// raft::copy(solver.pdhg_solver_.get_primal_solution().data(), +// d_initial_primal.data(), +// d_initial_primal.size(), +// handle.get_stream()); + +// auto& current_termination_strategy = solver.get_current_termination_strategy(); +// pdlp_termination_status_t termination_average = +// current_termination_strategy.evaluate_termination_criteria( +// solver.pdhg_solver_, +// d_initial_primal, +// d_initial_primal, +// solver.pdhg_solver_.get_dual_slack(), +// problem.combined_bounds, +// problem.objective_coefficients); +// EXPECT_EQ(current_termination_strategy.get_convergence_information() +// .get_relative_linf_primal_residual() +// .value(handle.get_stream()), +// 0.1); +// } +// } + +// TEST(pdlp_class, best_primal_so_far_iteration) +// { +// GTEST_SKIP() << "Skipping test: best_primal_so_far_iteration. Enable when ready to run."; +// const raft::handle_t handle1{}; +// const raft::handle_t handle2{}; + +// auto path = make_path_absolute("linear_programming/ns1687037/ns1687037.mps"); +// auto solver_settings = pdlp_solver_settings_t{}; +// solver_settings.iteration_limit = 3000; +// solver_settings.per_constraint_residual = true; +// solver_settings.method = cuopt::linear_programming::method_t::PDLP; + +// cuopt::mps_parser::mps_data_model_t op_problem1 = +// cuopt::mps_parser::parse_mps(path); +// cuopt::mps_parser::mps_data_model_t op_problem2 = +// cuopt::mps_parser::parse_mps(path); + +// optimization_problem_solution_t solution1 = +// solve_lp(&handle1, op_problem1, solver_settings); +// RAFT_CUDA_TRY(cudaDeviceSynchronize()); +// solver_settings.save_best_primal_so_far = true; +// optimization_problem_solution_t solution2 = +// solve_lp(&handle2, op_problem2, solver_settings); +// RAFT_CUDA_TRY(cudaDeviceSynchronize()); + +// EXPECT_TRUE(solution2.get_additional_termination_information().l2_primal_residual < +// solution1.get_additional_termination_information().l2_primal_residual); +// } + +// TEST(pdlp_class, best_primal_so_far_time) +// { +// GTEST_SKIP() << "Skipping test: best_primal_so_far_time. Enable when ready to run."; +// const raft::handle_t handle1{}; +// const raft::handle_t handle2{}; + +// auto path = make_path_absolute("linear_programming/ns1687037/ns1687037.mps"); +// auto solver_settings = pdlp_solver_settings_t{}; +// solver_settings.time_limit = 2; +// solver_settings.per_constraint_residual = true; +// solver_settings.pdlp_solver_mode = +// cuopt::linear_programming::pdlp_solver_mode_t::Stable1; solver_settings.method = +// cuopt::linear_programming::method_t::PDLP; + +// cuopt::mps_parser::mps_data_model_t op_problem1 = +// cuopt::mps_parser::parse_mps(path); +// cuopt::mps_parser::mps_data_model_t op_problem2 = +// cuopt::mps_parser::parse_mps(path); + +// optimization_problem_solution_t solution1 = +// solve_lp(&handle1, op_problem1, solver_settings); +// RAFT_CUDA_TRY(cudaDeviceSynchronize()); +// solver_settings.save_best_primal_so_far = true; +// optimization_problem_solution_t solution2 = +// solve_lp(&handle2, op_problem2, solver_settings); +// RAFT_CUDA_TRY(cudaDeviceSynchronize()); + +// EXPECT_TRUE(solution2.get_additional_termination_information().l2_primal_residual < +// solution1.get_additional_termination_information().l2_primal_residual); +// } + +// TEST(pdlp_class, first_primal_feasible) +// { +// GTEST_SKIP() << "Skipping test: first_primal_feasible. Enable when ready to run."; +// const raft::handle_t handle1{}; +// const raft::handle_t handle2{}; + +// auto path = make_path_absolute("linear_programming/ns1687037/ns1687037.mps"); +// auto solver_settings = pdlp_solver_settings_t{}; +// solver_settings.iteration_limit = 1000; +// solver_settings.per_constraint_residual = true; +// solver_settings.set_optimality_tolerance(1e-2); +// solver_settings.method = cuopt::linear_programming::method_t::PDLP; + +// cuopt::mps_parser::mps_data_model_t op_problem1 = +// cuopt::mps_parser::parse_mps(path); +// cuopt::mps_parser::mps_data_model_t op_problem2 = +// cuopt::mps_parser::parse_mps(path); + +// optimization_problem_solution_t solution1 = +// solve_lp(&handle1, op_problem1, solver_settings); +// RAFT_CUDA_TRY(cudaDeviceSynchronize()); +// solver_settings.first_primal_feasible = true; +// optimization_problem_solution_t solution2 = +// solve_lp(&handle2, op_problem2, solver_settings); +// RAFT_CUDA_TRY(cudaDeviceSynchronize()); + +// EXPECT_EQ(solution1.get_termination_status(), pdlp_termination_status_t::IterationLimit); +// EXPECT_EQ(solution2.get_termination_status(), pdlp_termination_status_t::PrimalFeasible); +// } + +// TEST(pdlp_class, warm_start) +// { +// std::vector instance_names{"graph40-40", +// "ex10", +// "datt256_lp", +// "woodlands09", +// "savsched1", +// // "nug08-3rd", // TODO: Fix this instance +// "qap15", +// "scpm1", +// // "neos3", // TODO: Fix this instance +// "a2864"}; +// for (auto instance_name : instance_names) { +// const raft::handle_t handle{}; + +// auto path = +// make_path_absolute("linear_programming/" + instance_name + "/" + instance_name + ".mps"); +// auto solver_settings = pdlp_solver_settings_t{}; +// solver_settings.pdlp_solver_mode = cuopt::linear_programming::pdlp_solver_mode_t::Stable2; +// solver_settings.set_optimality_tolerance(1e-2); +// solver_settings.detect_infeasibility = false; +// solver_settings.method = cuopt::linear_programming::method_t::PDLP; + +// cuopt::mps_parser::mps_data_model_t mps_data_model = +// cuopt::mps_parser::parse_mps(path); +// auto op_problem1 = +// cuopt::linear_programming::mps_data_model_to_optimization_problem( +// &handle, mps_data_model); + +// // Solving from scratch until 1e-2 +// optimization_problem_solution_t solution1 = solve_lp(op_problem1, +// solver_settings); + +// // Solving until 1e-1 to use the result as a warm start +// solver_settings.set_optimality_tolerance(1e-1); +// auto op_problem2 = +// cuopt::linear_programming::mps_data_model_to_optimization_problem( +// &handle, mps_data_model); +// optimization_problem_solution_t solution2 = solve_lp(op_problem2, +// solver_settings); + +// // Solving until 1e-2 using the previous state as a warm start +// solver_settings.set_optimality_tolerance(1e-2); +// auto op_problem3 = +// cuopt::linear_programming::mps_data_model_to_optimization_problem( +// &handle, mps_data_model); +// solver_settings.set_pdlp_warm_start_data(solution2.get_pdlp_warm_start_data()); +// optimization_problem_solution_t solution3 = solve_lp(op_problem3, +// solver_settings); + +// EXPECT_EQ(solution1.get_additional_termination_information().number_of_steps_taken, +// solution3.get_additional_termination_information().number_of_steps_taken + +// solution2.get_additional_termination_information().number_of_steps_taken); +// } +// } + +// TEST(pdlp_class, dual_postsolve_size) +// { +// const raft::handle_t handle_{}; + +// auto path = make_path_absolute("linear_programming/afiro_original.mps"); +// cuopt::mps_parser::mps_data_model_t op_problem = +// cuopt::mps_parser::parse_mps(path, true); + +// auto solver_settings = pdlp_solver_settings_t{}; +// solver_settings.method = cuopt::linear_programming::method_t::PDLP; +// solver_settings.presolve = true; + +// { +// solver_settings.dual_postsolve = true; +// optimization_problem_solution_t solution = +// solve_lp(&handle_, op_problem, solver_settings); +// EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); +// EXPECT_EQ(solution.get_dual_solution().size(), op_problem.get_n_constraints()); +// } + +// { +// solver_settings.dual_postsolve = false; +// optimization_problem_solution_t solution = +// solve_lp(&handle_, op_problem, solver_settings); +// EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); +// EXPECT_EQ(solution.get_dual_solution().size(), 0); +// } +// } + +// TEST(dual_simplex, afiro) +// { +// cuopt::linear_programming::pdlp_solver_settings_t settings = +// cuopt::linear_programming::pdlp_solver_settings_t{}; +// settings.method = cuopt::linear_programming::method_t::DualSimplex; + +// const raft::handle_t handle_{}; + +// auto path = make_path_absolute("linear_programming/afiro_original.mps"); +// cuopt::mps_parser::mps_data_model_t op_problem = +// cuopt::mps_parser::parse_mps(path, true); + +// optimization_problem_solution_t solution = solve_lp(&handle_, op_problem, +// settings); EXPECT_EQ(solution.get_termination_status(), pdlp_termination_status_t::Optimal); +// EXPECT_FALSE(is_incorrect_objective( +// afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); +// } // Should return a numerical error TEST(pdlp_class, run_empty_matrix_pdlp)