From 4c3fc360471d42accb88b3e20b4a3119c7cd6908 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Wed, 24 Sep 2025 15:41:18 +0000 Subject: [PATCH 01/15] initial testing --- benchmarks/linear_programming/cuopt/run_mip.cpp | 4 ++++ cpp/src/dual_simplex/CMakeLists.txt | 2 ++ cpp/src/dual_simplex/basis_updates.cpp | 12 +++++++----- cpp/src/dual_simplex/branch_and_bound.cpp | 1 + cpp/src/dual_simplex/phase2.cpp | 3 ++- cpp/src/dual_simplex/presolve.cpp | 2 ++ 6 files changed, 18 insertions(+), 6 deletions(-) diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index 6013dcaf5..1f03707bc 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -38,6 +38,8 @@ #include #include +#include + #include "initial_problem_check.hpp" void merge_result_files(const std::string& out_dir, @@ -291,6 +293,8 @@ void return_gpu_to_the_queue(std::unordered_map& pid_gpu_map, int main(int argc, char* argv[]) { + feenableexcept(FE_DIVBYZERO | FE_INVALID); + argparse::ArgumentParser program("solve_MIP"); // Define all arguments with appropriate defaults and help messages diff --git a/cpp/src/dual_simplex/CMakeLists.txt b/cpp/src/dual_simplex/CMakeLists.txt index e8a9b5dce..87a7076cb 100644 --- a/cpp/src/dual_simplex/CMakeLists.txt +++ b/cpp/src/dual_simplex/CMakeLists.txt @@ -36,5 +36,7 @@ set(DUAL_SIMPLEX_SRC_FILES # Uncomment to enable debug info #set_source_files_properties(${DUAL_SIMPLEX_SRC_FILES} DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g1") +set_source_files_properties(${DUAL_SIMPLEX_SRC_FILES} DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g;-Og") + set(CUOPT_SRC_FILES ${CUOPT_SRC_FILES} ${DUAL_SIMPLEX_SRC_FILES} PARENT_SCOPE) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 6b79f3c86..3f98d1132 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1620,11 +1620,12 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, sparse_vector_t& Lsol, bool need_Lsol) const { + // assert(rhs.n > 0); const i_t m = L0_.m; solution = rhs; solution.inverse_permute_vector(inverse_row_permutation_); -#ifdef CHECK_PERMUTATION +#if defined(CHECK_PERMUTATION) std::vector permuation_rhs; rhs.to_dense(permuation_rhs); std::vector finish_perm(m); @@ -1639,7 +1640,7 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, } #endif -#ifdef CHECK_L_SOLVE +#if defined(CHECK_L_SOLVE) std::vector l_solve_rhs; solution.to_dense(l_solve_rhs); #endif @@ -1658,7 +1659,7 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, if (need_Lsol) { Lsol = solution; } sum_L_ += static_cast(solution.i.size()) / input_size; -#ifdef CHECK_L_SOLVE +#if defined(CHECK_L_SOLVE) std::vector l_solve_dense; Lsol.to_dense(l_solve_dense); @@ -1671,12 +1672,13 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, if (max_err_l_solve > 1e-9) { printf("B solve L solve residual %e\n", max_err_l_solve); } #endif -#ifdef CHECK_U_SOLVE +#if defined(CHECK_U_SOLVE) std::vector rhs_dense; solution.to_dense(rhs_dense); #endif const f_t rhs_size = static_cast(solution.i.size()); + // assert(rhs_size > 0); estimate_solution_density(rhs_size, sum_U_, num_calls_U_, use_hypersparse); if (use_hypersparse) { u_solve(solution); @@ -1688,7 +1690,7 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, } sum_U_ += static_cast(solution.i.size()) / rhs_size; -#ifdef CHECK_U_SOLVE +#if defined(CHECK_U_SOLVE) std::vector solution_dense; solution.to_dense(solution_dense); diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 77acca8f7..53c757d46 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -159,6 +159,7 @@ f_t sgn(f_t x) template f_t relative_gap(f_t obj_value, f_t lower_bound) { + if (!std::isfinite(obj_value)) 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); diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 56298ef4d..9486d0c2e 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -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; diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index d247fbf67..bd016a264 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -13,6 +13,8 @@ #include #include +#include + #include #include #include From 3380d72a3e90a25acf6f5d41d0ea56cfafe2949a Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Mon, 6 Oct 2025 09:35:50 +0000 Subject: [PATCH 02/15] FJCPU inf timelimit fix --- cpp/src/mip/feasibility_jump/fj_cpu.cu | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/cpp/src/mip/feasibility_jump/fj_cpu.cu b/cpp/src/mip/feasibility_jump/fj_cpu.cu index 8d534dfff..e9f75022f 100644 --- a/cpp/src/mip/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip/feasibility_jump/fj_cpu.cu @@ -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, From 256836b5a2c046f6f0476507bab761971983550f Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Mon, 6 Oct 2025 12:58:05 +0000 Subject: [PATCH 03/15] catch FPEs in b_solve --- cpp/src/dual_simplex/basis_updates.cpp | 15 +++++++++++++++ cpp/src/dual_simplex/branch_and_bound.cpp | 15 +++++++++++---- 2 files changed, 26 insertions(+), 4 deletions(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 3f98d1132..247a96b7e 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -13,6 +13,19 @@ #include #include +#include + +class fpe_disable { + int old_mask; + + public: + explicit fpe_disable(int mask = FE_INVALID) : old_mask(fegetexcept()) { fedisableexcept(mask); } + ~fpe_disable() { feenableexcept(old_mask); } + + fpe_disable(const fpe_disable&) = delete; + fpe_disable& operator=(const fpe_disable&) = delete; +}; + namespace cuopt::linear_programming::dual_simplex { template @@ -1620,6 +1633,8 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, sparse_vector_t& Lsol, bool need_Lsol) const { + fpe_disable fpe_disabler; + // assert(rhs.n > 0); const i_t m = L0_.m; solution = rhs; diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 53c757d46..5cf60aacd 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -159,7 +159,10 @@ f_t sgn(f_t x) template f_t relative_gap(f_t obj_value, f_t lower_bound) { - if (!std::isfinite(obj_value)) return std::numeric_limits::infinity(); + // 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); @@ -172,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; } From 659fd5bf6c6000657b47d1d9191bac96ee32ce39 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Mon, 6 Oct 2025 14:57:04 +0000 Subject: [PATCH 04/15] more fixes --- .../mip/feasibility_jump/feasibility_jump_impl_common.cuh | 4 ++-- .../mip/local_search/feasibility_pump/feasibility_pump.cu | 7 ++++--- 2 files changed, 6 insertions(+), 5 deletions(-) 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..920fc8b48 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh +++ b/cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh @@ -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/local_search/feasibility_pump/feasibility_pump.cu b/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu index 76d930b45..debb79883 100644 --- a/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu +++ b/cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu @@ -85,9 +85,10 @@ void feasibility_pump_t::adjust_objective_with_original(solution_t Date: Fri, 12 Dec 2025 09:50:58 +0000 Subject: [PATCH 05/15] typos, relax FPE --- cpp/src/dual_simplex/basis_updates.cpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 247a96b7e..65517d51c 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1633,14 +1633,13 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, sparse_vector_t& Lsol, bool need_Lsol) const { - fpe_disable fpe_disabler; + // fpe_disable fpe_disabler; - // assert(rhs.n > 0); const i_t m = L0_.m; solution = rhs; solution.inverse_permute_vector(inverse_row_permutation_); -#if defined(CHECK_PERMUTATION) +#ifdef CHECK_PERMUTATION std::vector permuation_rhs; rhs.to_dense(permuation_rhs); std::vector finish_perm(m); @@ -1655,7 +1654,7 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, } #endif -#if defined(CHECK_L_SOLVE) +#ifdef CHECK_L_SOLVE std::vector l_solve_rhs; solution.to_dense(l_solve_rhs); #endif @@ -1674,7 +1673,7 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, if (need_Lsol) { Lsol = solution; } sum_L_ += static_cast(solution.i.size()) / input_size; -#if defined(CHECK_L_SOLVE) +#ifdef CHECK_L_SOLVE std::vector l_solve_dense; Lsol.to_dense(l_solve_dense); @@ -1687,13 +1686,12 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, if (max_err_l_solve > 1e-9) { printf("B solve L solve residual %e\n", max_err_l_solve); } #endif -#if defined(CHECK_U_SOLVE) +#ifdef CHECK_U_SOLVE std::vector rhs_dense; solution.to_dense(rhs_dense); #endif const f_t rhs_size = static_cast(solution.i.size()); - // assert(rhs_size > 0); estimate_solution_density(rhs_size, sum_U_, num_calls_U_, use_hypersparse); if (use_hypersparse) { u_solve(solution); @@ -1705,7 +1703,7 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, } sum_U_ += static_cast(solution.i.size()) / rhs_size; -#if defined(CHECK_U_SOLVE) +#ifdef CHECK_U_SOLVE std::vector solution_dense; solution.to_dense(solution_dense); From d0199d493d2862c35919e3800f0851c074ba1a10 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Thu, 8 Jan 2026 15:24:07 +0000 Subject: [PATCH 06/15] more FPE fixes --- .../linear_programming/cuopt/run_mip.cpp | 4 ---- cpp/CMakeLists.txt | 2 ++ cpp/src/dual_simplex/CMakeLists.txt | 6 ++---- cpp/src/dual_simplex/basis_updates.cpp | 15 --------------- cpp/src/dual_simplex/bounds_strengthening.cpp | 19 ++++++++++++++++++- cpp/src/dual_simplex/branch_and_bound.cpp | 6 ++++-- cpp/src/dual_simplex/phase2.cpp | 2 +- cpp/src/dual_simplex/presolve.cpp | 2 -- .../restart_strategy/pdlp_restart_strategy.cu | 17 +++++++++++------ cpp/src/linear_programming/solve.cu | 9 ++++++++- cpp/src/mip/solve.cu | 7 +++++++ 11 files changed, 53 insertions(+), 36 deletions(-) diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index 1f03707bc..6013dcaf5 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -38,8 +38,6 @@ #include #include -#include - #include "initial_problem_check.hpp" void merge_result_files(const std::string& out_dir, @@ -293,8 +291,6 @@ void return_gpu_to_the_queue(std::unordered_map& pid_gpu_map, int main(int argc, char* argv[]) { - feenableexcept(FE_DIVBYZERO | FE_INVALID); - argparse::ArgumentParser program("solve_MIP"); // Define all arguments with appropriate defaults and help messages 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/CMakeLists.txt b/cpp/src/dual_simplex/CMakeLists.txt index 87a7076cb..799b9b5a9 100644 --- a/cpp/src/dual_simplex/CMakeLists.txt +++ b/cpp/src/dual_simplex/CMakeLists.txt @@ -1,5 +1,5 @@ # cmake-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 # cmake-format: on @@ -34,9 +34,7 @@ set(DUAL_SIMPLEX_SRC_FILES ) # Uncomment to enable debug info -#set_source_files_properties(${DUAL_SIMPLEX_SRC_FILES} DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g1") - -set_source_files_properties(${DUAL_SIMPLEX_SRC_FILES} DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g;-Og") +#set_source_files_properties(${DUAL_SIMPLEX_SRC_FILES} DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g2") set(CUOPT_SRC_FILES ${CUOPT_SRC_FILES} ${DUAL_SIMPLEX_SRC_FILES} PARENT_SCOPE) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 65517d51c..6b79f3c86 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -13,19 +13,6 @@ #include #include -#include - -class fpe_disable { - int old_mask; - - public: - explicit fpe_disable(int mask = FE_INVALID) : old_mask(fegetexcept()) { fedisableexcept(mask); } - ~fpe_disable() { feenableexcept(old_mask); } - - fpe_disable(const fpe_disable&) = delete; - fpe_disable& operator=(const fpe_disable&) = delete; -}; - namespace cuopt::linear_programming::dual_simplex { template @@ -1633,8 +1620,6 @@ i_t basis_update_mpf_t::b_solve(const sparse_vector_t& rhs, sparse_vector_t& Lsol, bool need_Lsol) const { - // fpe_disable fpe_disabler; - const i_t m = L0_.m; solution = rhs; solution.inverse_permute_vector(inverse_row_permutation_); diff --git a/cpp/src/dual_simplex/bounds_strengthening.cpp b/cpp/src/dual_simplex/bounds_strengthening.cpp index f1bf52c1e..a6efd2461 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,19 @@ #include #include +#include + +class fpe_disable { + int old_mask; + + public: + explicit fpe_disable(int mask = FE_INVALID) : old_mask(fegetexcept()) { fedisableexcept(mask); } + ~fpe_disable() { 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 +108,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 a5ad5df96..f6eebe44d 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 */ @@ -471,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); diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 9486d0c2e..f6b22522e 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 */ diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index bd016a264..d247fbf67 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -13,8 +13,6 @@ #include #include -#include - #include #include #include 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..f15df6a6c 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_); @@ -817,7 +823,6 @@ void pdlp_restart_strategy_t::cupdlpx_restart( weighted_average_solution_.iterations_since_last_restart_ = 0; last_trial_fixed_point_error_ = std::numeric_limits::infinity(); } - template bool pdlp_restart_strategy_t::run_cupdlpx_restart( const convergence_information_t& current_convergence_information, diff --git a/cpp/src/linear_programming/solve.cu b/cpp/src/linear_programming/solve.cu index d038ade72..3c12598cd 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,6 +42,8 @@ #include // For std::thread +#include + namespace cuopt::linear_programming { // This serves as both a warm up but also a mandatory initial call to setup cuSparse and cuBLAS @@ -827,6 +829,11 @@ optimization_problem_solution_t solve_lp( // This needs to be called before pdlp is initialized init_handler(op_problem.get_handle_ptr()); +#ifndef NDEBUG + CUOPT_LOG_DEBUG("Enabling host FPEs"); + feenableexcept(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/solve.cu b/cpp/src/mip/solve.cu index e6a392d40..ed0c66f9f 100644 --- a/cpp/src/mip/solve.cu +++ b/cpp/src/mip/solve.cu @@ -38,6 +38,8 @@ #include +#include + namespace cuopt::linear_programming { // This serves as both a warm up but also a mandatory initial call to setup cuSparse and cuBLAS @@ -160,6 +162,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()); +#ifndef NDEBUG + CUOPT_LOG_DEBUG("Enabling host FPEs"); + feenableexcept(FE_DIVBYZERO | FE_INVALID); +#endif + print_version_info(); raft::common::nvtx::range fun_scope("Running solver"); From 395fd1ce32416b09614ca2657167592da70069be Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Thu, 8 Jan 2026 15:25:25 +0000 Subject: [PATCH 07/15] cleanup --- cpp/src/dual_simplex/CMakeLists.txt | 4 ++-- .../restart_strategy/pdlp_restart_strategy.cu | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/CMakeLists.txt b/cpp/src/dual_simplex/CMakeLists.txt index 799b9b5a9..e8a9b5dce 100644 --- a/cpp/src/dual_simplex/CMakeLists.txt +++ b/cpp/src/dual_simplex/CMakeLists.txt @@ -1,5 +1,5 @@ # cmake-format: off -# SPDX-FileCopyrightText: Copyright (c) 2024-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-FileCopyrightText: Copyright (c) 2024-2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. # SPDX-License-Identifier: Apache-2.0 # cmake-format: on @@ -34,7 +34,7 @@ set(DUAL_SIMPLEX_SRC_FILES ) # Uncomment to enable debug info -#set_source_files_properties(${DUAL_SIMPLEX_SRC_FILES} DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g2") +#set_source_files_properties(${DUAL_SIMPLEX_SRC_FILES} DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g1") set(CUOPT_SRC_FILES ${CUOPT_SRC_FILES} ${DUAL_SIMPLEX_SRC_FILES} PARENT_SCOPE) 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 f15df6a6c..4e7449289 100644 --- a/cpp/src/linear_programming/restart_strategy/pdlp_restart_strategy.cu +++ b/cpp/src/linear_programming/restart_strategy/pdlp_restart_strategy.cu @@ -823,6 +823,7 @@ void pdlp_restart_strategy_t::cupdlpx_restart( weighted_average_solution_.iterations_since_last_restart_ = 0; last_trial_fixed_point_error_ = std::numeric_limits::infinity(); } + template bool pdlp_restart_strategy_t::run_cupdlpx_restart( const convergence_information_t& current_convergence_information, From 9848ce25919dbe772dd31ab87205a4f25a2c8db8 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Fri, 9 Jan 2026 09:18:04 +0000 Subject: [PATCH 08/15] RAII guarding of FPE to prevent contamination across tests --- cpp/src/dual_simplex/bounds_strengthening.cpp | 6 +++++- cpp/src/linear_programming/solve.cu | 17 ++++++++++++++++- cpp/src/mip/solve.cu | 17 ++++++++++++++++- 3 files changed, 37 insertions(+), 3 deletions(-) diff --git a/cpp/src/dual_simplex/bounds_strengthening.cpp b/cpp/src/dual_simplex/bounds_strengthening.cpp index a6efd2461..826d901c0 100644 --- a/cpp/src/dual_simplex/bounds_strengthening.cpp +++ b/cpp/src/dual_simplex/bounds_strengthening.cpp @@ -17,7 +17,11 @@ class fpe_disable { public: explicit fpe_disable(int mask = FE_INVALID) : old_mask(fegetexcept()) { fedisableexcept(mask); } - ~fpe_disable() { feenableexcept(old_mask); } + ~fpe_disable() + { + fedisableexcept(FE_ALL_EXCEPT); + feenableexcept(old_mask); + } fpe_disable(const fpe_disable&) = delete; fpe_disable& operator=(const fpe_disable&) = delete; diff --git a/cpp/src/linear_programming/solve.cu b/cpp/src/linear_programming/solve.cu index 3c12598cd..9ec691d26 100644 --- a/cpp/src/linear_programming/solve.cu +++ b/cpp/src/linear_programming/solve.cu @@ -46,6 +46,21 @@ 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) { @@ -831,7 +846,7 @@ optimization_problem_solution_t solve_lp( #ifndef NDEBUG CUOPT_LOG_DEBUG("Enabling host FPEs"); - feenableexcept(FE_DIVBYZERO | FE_INVALID); + fpe_enable fpe_guard(FE_DIVBYZERO | FE_INVALID); #endif if (op_problem.has_quadratic_objective()) { diff --git a/cpp/src/mip/solve.cu b/cpp/src/mip/solve.cu index ed0c66f9f..78d8584e4 100644 --- a/cpp/src/mip/solve.cu +++ b/cpp/src/mip/solve.cu @@ -42,6 +42,21 @@ 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) { @@ -164,7 +179,7 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, #ifndef NDEBUG CUOPT_LOG_DEBUG("Enabling host FPEs"); - feenableexcept(FE_DIVBYZERO | FE_INVALID); + fpe_enable fpe_guard(FE_DIVBYZERO | FE_INVALID); #endif print_version_info(); From c950868869b2ca171c9741159dba9ab5f686eb9b Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Fri, 9 Jan 2026 09:30:04 +0000 Subject: [PATCH 09/15] fix copyrights --- cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh | 2 +- cpp/src/mip/feasibility_jump/fj_cpu.cu | 2 +- cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) 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 920fc8b48..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 */ diff --git a/cpp/src/mip/feasibility_jump/fj_cpu.cu b/cpp/src/mip/feasibility_jump/fj_cpu.cu index e9f75022f..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 */ 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 debb79883..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 */ From 16ecc396faede76c5ce3b4e1e9a9f90bd51d3f17 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Fri, 9 Jan 2026 17:06:18 +0000 Subject: [PATCH 10/15] debug signalign nan issue --- cpp/src/linear_programming/solve.cu | 4 +- cpp/tests/linear_programming/pdlp_test.cu | 1714 +++++++++++---------- 2 files changed, 874 insertions(+), 844 deletions(-) diff --git a/cpp/src/linear_programming/solve.cu b/cpp/src/linear_programming/solve.cu index 9ec691d26..a0e496aa9 100644 --- a/cpp/src/linear_programming/solve.cu +++ b/cpp/src/linear_programming/solve.cu @@ -572,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()); 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) From af77a56678cf8ca39cad9a798b3af890730205ab Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Fri, 9 Jan 2026 17:24:15 +0000 Subject: [PATCH 11/15] fpe always --- cpp/src/linear_programming/solve.cu | 2 +- cpp/src/mip/solve.cu | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/linear_programming/solve.cu b/cpp/src/linear_programming/solve.cu index a0e496aa9..a2c945218 100644 --- a/cpp/src/linear_programming/solve.cu +++ b/cpp/src/linear_programming/solve.cu @@ -846,7 +846,7 @@ optimization_problem_solution_t solve_lp( // This needs to be called before pdlp is initialized init_handler(op_problem.get_handle_ptr()); -#ifndef NDEBUG +#if 1 CUOPT_LOG_DEBUG("Enabling host FPEs"); fpe_enable fpe_guard(FE_DIVBYZERO | FE_INVALID); #endif diff --git a/cpp/src/mip/solve.cu b/cpp/src/mip/solve.cu index 78d8584e4..b3bba34a4 100644 --- a/cpp/src/mip/solve.cu +++ b/cpp/src/mip/solve.cu @@ -177,7 +177,7 @@ 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()); -#ifndef NDEBUG +#if 1 CUOPT_LOG_DEBUG("Enabling host FPEs"); fpe_enable fpe_guard(FE_DIVBYZERO | FE_INVALID); #endif From 976e0b66b39fb7c31eeec78b9890662f2593138d Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Mon, 12 Jan 2026 13:50:02 +0000 Subject: [PATCH 12/15] handle div-by-0 and 0/0 in b_solve --- cpp/src/dual_simplex/basis_updates.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 6b79f3c86..20eb145da 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 */ @@ -1656,7 +1656,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 +1677,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); From 5436ed333630a9165e2613d9f6fbdf4b17c12947 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 13 Jan 2026 16:35:50 +0000 Subject: [PATCH 13/15] basis update guards --- cpp/src/dual_simplex/basis_updates.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 20eb145da..c81820dd4 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -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; From ca3007fe295bfa0199c20e2129b5a0cb65588e28 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Thu, 8 Jan 2026 12:59:43 +0000 Subject: [PATCH 14/15] fix incorrect termination code when user termination is requested and incorrect PDLP objective problem space when passed to B&B from diversity_manager --- cpp/src/dual_simplex/basis_solves.cpp | 6 +-- cpp/src/dual_simplex/basis_updates.cpp | 47 +++++++++++++--------- cpp/src/dual_simplex/branch_and_bound.cpp | 5 +++ cpp/src/dual_simplex/crossover.cpp | 23 +++++++---- cpp/src/dual_simplex/phase2.cpp | 11 +++-- cpp/src/dual_simplex/primal.cpp | 9 +++-- cpp/src/dual_simplex/primal.hpp | 5 ++- cpp/src/mip/diversity/diversity_manager.cu | 17 ++++---- cpp/src/mip/problem/problem.cu | 6 +++ cpp/src/mip/problem/problem.cuh | 3 +- 10 files changed, 83 insertions(+), 49 deletions(-) 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 c81820dd4..7b4a4735c 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -2059,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); @@ -2089,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/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index f6eebe44d..47c13e1ba 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -1364,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 && + settings_.check_termination(exploration_stats_.start_time)) { + 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 f6b22522e..097ccc913 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -2242,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; } } @@ -2884,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; } @@ -2902,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/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/problem/problem.cu b/cpp/src/mip/problem/problem.cu index 815ef5aa0..635796cc9 100644 --- a/cpp/src/mip/problem/problem.cu +++ b/cpp/src/mip/problem/problem.cu @@ -1780,6 +1780,12 @@ f_t problem_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, From 66aae74f7915d7eb6201af761d7d7fa854424af9 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Wed, 14 Jan 2026 10:29:36 +0000 Subject: [PATCH 15/15] fix cherrypick --- cpp/src/dual_simplex/branch_and_bound.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 47c13e1ba..4ddaf4ac7 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -1365,7 +1365,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (root_status == lp_status_t::CONCURRENT_LIMIT && - settings_.check_termination(exploration_stats_.start_time)) { + toc(exploration_stats_.start_time) > settings_.time_limit) { solver_status_ = mip_exploration_status_t::TIME_LIMIT; return set_final_solution(solution, root_objective_); }