Enable Floating-Point exceptions on the host

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)

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<i_t, f_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<i_t, f_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));

cpp/src/dual_simplex/basis_updates.cpp

@@ -1,6 +1,6 @@
 /* clang-format off */
 /*
- * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
  * SPDX-License-Identifier: Apache-2.0
  */
 /* clang-format on */
@@ -1354,7 +1354,8 @@ i_t basis_update_mpf_t<i_t, f_t>::b_transpose_solve(const sparse_vector_t<i_t, f
     solution.from_dense(solution_dense);
   }
   UTsol = solution;
-  sum_U_transpose_ += static_cast<f_t>(solution.i.size()) / input_size;
+  // avoid div-by-zero FPE
+  if (input_size > 0.0) { sum_U_transpose_ += static_cast<f_t>(solution.i.size()) / input_size; }
 
 #ifdef CHECK_U_TRANSPOSE_SOLVE
   std::vector<f_t> UTsol_dense;
@@ -1383,7 +1384,8 @@ i_t basis_update_mpf_t<i_t, f_t>::b_transpose_solve(const sparse_vector_t<i_t, f
     l_transpose_solve(solution_dense);
     solution.from_dense(solution_dense);
   }
-  sum_L_transpose_ += static_cast<f_t>(solution.i.size()) / rhs_size;
+  // avoid div-by-zero FPE
+  if (rhs_size > 0.0) { sum_L_transpose_ += static_cast<f_t>(solution.i.size()) / rhs_size; }
 
 #ifdef CHECK_L_TRANSPOSE_SOLVE
   std::vector<f_t> solution_dense;
@@ -1656,7 +1658,7 @@ i_t basis_update_mpf_t<i_t, f_t>::b_solve(const sparse_vector_t<i_t, f_t>& rhs,
     solution.from_dense(solution_dense);
   }
   if (need_Lsol) { Lsol = solution; }
-  sum_L_ += static_cast<f_t>(solution.i.size()) / input_size;
+  if (input_size > 0.0) { sum_L_ += static_cast<f_t>(solution.i.size()) / input_size; }
 
 #ifdef CHECK_L_SOLVE
   std::vector<f_t> l_solve_dense;
@@ -1677,6 +1679,8 @@ i_t basis_update_mpf_t<i_t, f_t>::b_solve(const sparse_vector_t<i_t, f_t>& rhs,
 #endif
 
   const f_t rhs_size = static_cast<f_t>(solution.i.size());
+  // (nothing to solve, avoids 0/0 in statistics)
+  if (rhs_size == 0.0) { return 0; }
   estimate_solution_density(rhs_size, sum_U_, num_calls_U_, use_hypersparse);
   if (use_hypersparse) {
     u_solve(solution);
@@ -2055,16 +2059,21 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
 
   if (L0_.m != A.m) { resize(A.m); }
   std::vector<i_t> q;
-  if (factorize_basis(A,
-                      settings,
-                      basic_list,
-                      L0_,
-                      U0_,
-                      row_permutation_,
-                      inverse_row_permutation_,
-                      q,
-                      deficient,
-                      slacks_needed) == -1) {
+  i_t factorize_result = factorize_basis(A,
+                                         settings,
+                                         basic_list,
+                                         L0_,
+                                         U0_,
+                                         row_permutation_,
+                                         inverse_row_permutation_,
+                                         q,
+                                         deficient,
+                                         slacks_needed);
+  if (factorize_result == -2) {
+    // Concurrent halt requested, return early
+    return -2;
+  }
+  if (factorize_result == -1) {
     settings.log.debug("Initial factorization failed\n");
     basis_repair(A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
 
@@ -2085,16 +2094,18 @@ int basis_update_mpf_t<i_t, f_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

cpp/src/dual_simplex/bounds_strengthening.cpp

@@ -1,6 +1,6 @@
 /* clang-format off */
 /*
- * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
  * SPDX-License-Identifier: Apache-2.0
  */
 /* clang-format on */
@@ -10,6 +10,23 @@
 #include <algorithm>
 #include <cmath>
 
+#include <fenv.h>
+
+class fpe_disable {
+  int old_mask;
+
+ public:
+  explicit fpe_disable(int mask = FE_INVALID) : old_mask(fegetexcept()) { fedisableexcept(mask); }
+  ~fpe_disable()
+  {
+    fedisableexcept(FE_ALL_EXCEPT);
+    feenableexcept(old_mask);
+  }
+
+  fpe_disable(const fpe_disable&)            = delete;
+  fpe_disable& operator=(const fpe_disable&) = delete;
+};
+
 namespace cuopt::linear_programming::dual_simplex {
 
 template <typename f_t>
@@ -95,6 +112,10 @@ bool bounds_strengthening_t<i_t, f_t>::bounds_strengthening(
   std::vector<f_t>& upper_bounds,
   const simplex_solver_settings_t<i_t, f_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;
 

cpp/src/dual_simplex/branch_and_bound.cpp

@@ -1,6 +1,6 @@
 /* clang-format off */
 /*
- * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
  * SPDX-License-Identifier: Apache-2.0
  */
 /* clang-format on */
@@ -159,6 +159,10 @@ f_t sgn(f_t x)
 template <typename f_t>
 f_t relative_gap(f_t obj_value, f_t lower_bound)
 {
+  // avoid raising a FPE
+  if (!std::isfinite(obj_value) || !std::isfinite(lower_bound)) {
+    return std::numeric_limits<f_t>::infinity();
+  }
   f_t user_mip_gap = obj_value == 0.0
                        ? (lower_bound == 0.0 ? 0.0 : std::numeric_limits<f_t>::infinity())
                        : std::abs(obj_value - lower_bound) / std::abs(obj_value);
@@ -171,9 +175,13 @@ f_t user_relative_gap(const lp_problem_t<i_t, f_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<f_t>::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<f_t>::infinity();
+  }
+  f_t user_mip_gap = user_obj == 0.0
+                       ? (user_lower_bound == 0.0 ? 0.0 : std::numeric_limits<f_t>::infinity())
+                       : std::abs(user_obj - user_lower_bound) / std::abs(user_obj);
   if (std::isnan(user_mip_gap)) { return std::numeric_limits<f_t>::infinity(); }
   return user_mip_gap;
 }
@@ -463,7 +471,9 @@ mip_status_t branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t
   }
 
   f_t upper_bound      = get_upper_bound();
-  f_t gap              = upper_bound - lower_bound;
+  f_t gap              = (std::isfinite(upper_bound) && std::isfinite(lower_bound))
+                           ? upper_bound - lower_bound
+                           : std::numeric_limits<f_t>::infinity();
   f_t obj              = compute_user_objective(original_lp_, upper_bound);
   f_t user_bound       = compute_user_objective(original_lp_, lower_bound);
   f_t gap_rel          = user_relative_gap(original_lp_, upper_bound, lower_bound);
@@ -1354,6 +1364,11 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
     return mip_status_t::UNBOUNDED;
   }
 
+  if (root_status == lp_status_t::CONCURRENT_LIMIT &&
+      toc(exploration_stats_.start_time) > settings_.time_limit) {
+    solver_status_ = mip_exploration_status_t::TIME_LIMIT;
+    return set_final_solution(solution, root_objective_);
+  }
   if (root_status == lp_status_t::TIME_LIMIT) {
     solver_status_ = mip_exploration_status_t::TIME_LIMIT;
     return set_final_solution(solution, -inf);

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<i_t, f_t>& lp,
         std::vector<i_t> 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<i_t, f_t>& lp,
   std::vector<i_t> 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<i_t, f_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 {

cpp/src/dual_simplex/phase2.cpp

@@ -1,6 +1,6 @@
 /* clang-format off */
 /*
- * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
  * SPDX-License-Identifier: Apache-2.0
  */
 /* clang-format on */
@@ -767,7 +767,8 @@ i_t steepest_edge_pricing_with_infeasibilities(const lp_problem_t<i_t, f_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<f_t>::infinity()
+                                                          : squared_infeas / dy_steepest_edge[j];
     if (val > max_val || (val == max_val && j > leaving_index)) {
       max_val                = val;
       leaving_index          = j;
@@ -2241,9 +2242,9 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
     assert(superbasic_list.size() == 0);
     assert(nonbasic_list.size() == n - m);
 
-    if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
-      return dual::status_t::NUMERICAL;
-    }
+    i_t refactor_result = ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus);
+    if (refactor_result == -2) { return dual::status_t::CONCURRENT_LIMIT; }
+    if (refactor_result > 0) { return dual::status_t::NUMERICAL; }
 
     if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
   }
@@ -2883,7 +2884,9 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
 #endif
     if (should_refactor) {
       bool should_recompute_x = false;
-      if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
+      i_t refactor_result = ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus);
+      if (refactor_result == -2) { return dual::status_t::CONCURRENT_LIMIT; }
+      if (refactor_result > 0) {
         should_recompute_x = true;
         settings.log.printf("Failed to factorize basis. Iteration %d\n", iter);
         if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
@@ -2901,6 +2904,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
           count++;
           if (count > 10) { return dual::status_t::NUMERICAL; }
         }
+        if (deficient_size == -2) { return dual::status_t::CONCURRENT_LIMIT; }
 
         settings.log.printf("Successfully repaired basis. Iteration %d\n", iter);
       }

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<i_t> 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 {

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
 };
 }