diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp
index ba2b63983e..41ed408379 100644
--- a/cpp/src/branch_and_bound/branch_and_bound.cpp
+++ b/cpp/src/branch_and_bound/branch_and_bound.cpp
@@ -6,8 +6,10 @@
 /* clang-format on */
 
 #include <branch_and_bound/branch_and_bound.hpp>
+#include <branch_and_bound/diving_heuristics.hpp>
 #include <branch_and_bound/mip_node.hpp>
 #include <branch_and_bound/pseudo_costs.hpp>
+#include <branch_and_bound/reduced_cost_fixing.hpp>
 
 #include <cuts/cuts.hpp>
 #include <mip_heuristics/presolve/conflict_graph/clique_table.cuh>
@@ -35,15 +37,12 @@
 #include <deque>
 #include <future>
 #include <limits>
-#include <map>
 #include <optional>
 #include <string>
 #include <thread>
-#include <unordered_map>
 #include <vector>
 
 namespace cuopt::linear_programming::dual_simplex {
-
 namespace {
 
 template <typename f_t>
@@ -196,6 +195,16 @@ f_t user_relative_gap(const lp_problem_t<i_t, f_t>& lp, f_t obj_value, f_t lower
   return user_mip_gap;
 }
 
+template <typename i_t, typename f_t>
+f_t user_absolute_gap(const lp_problem_t<i_t, f_t>& lp, f_t obj_value, f_t lower_bound)
+{
+  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     = std::abs(user_obj - user_lower_bound);
+  if (std::isnan(user_mip_gap)) { return std::numeric_limits<f_t>::infinity(); }
+  return user_mip_gap;
+}
+
 template <typename f_t>
 std::string user_mip_gap(f_t obj_value, f_t lower_bound)
 {
@@ -380,66 +389,6 @@ void branch_and_bound_t<i_t, f_t>::report(
   }
 }
 
-template <typename i_t, typename f_t>
-i_t branch_and_bound_t<i_t, f_t>::find_reduced_cost_fixings(f_t upper_bound,
-                                                            std::vector<f_t>& lower_bounds,
-                                                            std::vector<f_t>& upper_bounds)
-{
-  std::vector<f_t> reduced_costs = root_relax_soln_.z;
-  lower_bounds                   = original_lp_.lower;
-  upper_bounds                   = original_lp_.upper;
-  std::vector<bool> bounds_changed(original_lp_.num_cols, false);
-  const f_t root_obj    = compute_objective(original_lp_, root_relax_soln_.x);
-  const f_t threshold   = 100.0 * settings_.integer_tol;
-  const f_t weaken      = settings_.integer_tol;
-  const f_t fixed_tol   = settings_.fixed_tol;
-  i_t num_improved      = 0;
-  i_t num_fixed         = 0;
-  i_t num_cols_to_check = reduced_costs.size();  // Reduced costs will be smaller than the original
-                                                 // problem because we have added slacks for cuts
-  for (i_t j = 0; j < num_cols_to_check; j++) {
-    if (std::isfinite(reduced_costs[j]) && std::abs(reduced_costs[j]) > threshold) {
-      const f_t lower_j            = original_lp_.lower[j];
-      const f_t upper_j            = original_lp_.upper[j];
-      const f_t abs_gap            = upper_bound - root_obj;
-      f_t reduced_cost_upper_bound = upper_j;
-      f_t reduced_cost_lower_bound = lower_j;
-      if (lower_j > -inf && reduced_costs[j] > 0) {
-        const f_t new_upper_bound = lower_j + abs_gap / reduced_costs[j];
-        reduced_cost_upper_bound  = var_types_[j] == variable_type_t::INTEGER
-                                      ? std::floor(new_upper_bound + weaken)
-                                      : new_upper_bound;
-        if (reduced_cost_upper_bound < upper_j && var_types_[j] == variable_type_t::INTEGER) {
-          num_improved++;
-          upper_bounds[j]   = reduced_cost_upper_bound;
-          bounds_changed[j] = true;
-        }
-      }
-      if (upper_j < inf && reduced_costs[j] < 0) {
-        const f_t new_lower_bound = upper_j + abs_gap / reduced_costs[j];
-        reduced_cost_lower_bound  = var_types_[j] == variable_type_t::INTEGER
-                                      ? std::ceil(new_lower_bound - weaken)
-                                      : new_lower_bound;
-        if (reduced_cost_lower_bound > lower_j && var_types_[j] == variable_type_t::INTEGER) {
-          num_improved++;
-          lower_bounds[j]   = reduced_cost_lower_bound;
-          bounds_changed[j] = true;
-        }
-      }
-      if (var_types_[j] == variable_type_t::INTEGER &&
-          reduced_cost_upper_bound <= reduced_cost_lower_bound + fixed_tol) {
-        num_fixed++;
-      }
-    }
-  }
-
-  if (num_fixed > 0 || num_improved > 0) {
-    settings_.log.printf(
-      "Reduced costs: Found %d improved bounds and %d fixed variables\n", num_improved, num_fixed);
-  }
-  return num_fixed;
-}
-
 template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::update_user_bound(f_t lower_bound)
 {
@@ -449,23 +398,64 @@ void branch_and_bound_t<i_t, f_t>::update_user_bound(f_t lower_bound)
 }
 
 template <typename i_t, typename f_t>
-void branch_and_bound_t<i_t, f_t>::set_new_solution(const std::vector<f_t>& solution)
+void branch_and_bound_t<i_t, f_t>::root_reduced_cost_fixing(f_t upper_bound)
 {
+  if (settings_.reduced_cost_strengthening < 3) { return; }
+
   mutex_original_lp_.lock();
+  std::vector<f_t> lower = original_lp_.lower;
+  std::vector<f_t> upper = original_lp_.upper;
+  std::vector<bool> bounds_changed;
+  mutex_original_lp_.unlock();
+
+  auto [num_fixed, num_improved] = reduced_cost_fixing(root_relax_soln_.z,
+                                                       var_types_,
+                                                       settings_,
+                                                       root_objective_,
+                                                       upper_bound,
+                                                       lower,
+                                                       upper,
+                                                       bounds_changed);
+  if (num_fixed > 0 || num_improved > 0) {
+    std::vector<char> row_sense;
+    bounds_strengthening_t<i_t, f_t> node_presolve(original_lp_, Arow_, row_sense, var_types_);
+
+    mutex_original_lp_.lock();
+    original_lp_.lower = lower;
+    original_lp_.upper = upper;
+    bool feasible      = node_presolve.bounds_strengthening(
+      settings_, bounds_changed, original_lp_.lower, original_lp_.upper);
+    mutex_original_lp_.unlock();
+
+    worker_pool_.broadcast_root_bounds_change();
+
+    if (!feasible) {
+      settings_.log.printf(
+        "Bound strengthening failed when updating the bounds at the root node!\n");
+      solver_status_ = mip_status_t::NUMERICAL;
+    }
+  }
+}
+
+template <typename i_t, typename f_t>
+void branch_and_bound_t<i_t, f_t>::set_new_solution(const std::vector<f_t>& solution)
+{
   if (solution.size() != original_problem_.num_cols) {
     settings_.log.printf(
       "Solution size mismatch %ld %d\n", solution.size(), original_problem_.num_cols);
   }
+
+  mutex_original_lp_.lock();
   std::vector<f_t> crushed_solution;
   crush_primal_solution<i_t, f_t>(
     original_problem_, original_lp_, solution, new_slacks_, crushed_solution);
   f_t obj = compute_objective(original_lp_, crushed_solution);
   mutex_original_lp_.unlock();
-  bool is_feasible    = false;
-  bool attempt_repair = false;
-  mutex_upper_.lock();
+
+  bool is_feasible        = false;
+  bool attempt_repair     = false;
   f_t current_upper_bound = upper_bound_;
-  mutex_upper_.unlock();
+
   if (obj < current_upper_bound) {
     f_t primal_err;
     f_t bound_err;
@@ -484,6 +474,7 @@ void branch_and_bound_t<i_t, f_t>::set_new_solution(const std::vector<f_t>& solu
     if (is_feasible && obj < upper_bound_) {
       upper_bound_ = obj;
       incumbent_.set_incumbent_solution(obj, crushed_solution);
+      root_reduced_cost_fixing(obj);
     } else {
       attempt_repair         = true;
       constexpr bool verbose = false;
@@ -632,10 +623,15 @@ void branch_and_bound_t<i_t, f_t>::repair_heuristic_solutions()
     settings_.log.debug("Attempting to repair %ld injected solutions\n", to_repair.size());
     for (const std::vector<f_t>& uncrushed_solution : to_repair) {
       std::vector<f_t> crushed_solution;
+
+      mutex_original_lp_.lock();
       crush_primal_solution<i_t, f_t>(
         original_problem_, original_lp_, uncrushed_solution, new_slacks_, crushed_solution);
+      mutex_original_lp_.unlock();
+
       std::vector<f_t> repaired_solution;
       f_t repaired_obj;
+
       bool is_feasible =
         repair_solution(edge_norms_, crushed_solution, repaired_obj, repaired_solution);
       if (is_feasible) {
@@ -651,8 +647,9 @@ void branch_and_bound_t<i_t, f_t>::repair_heuristic_solutions()
             uncrush_primal_solution(original_problem_, original_lp_, repaired_solution, original_x);
             settings_.solution_callback(original_x, repaired_obj);
           }
-        }
 
+          root_reduced_cost_fixing(repaired_obj);
+        }
         mutex_upper_.unlock();
       }
     }
@@ -710,9 +707,9 @@ void branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t, f_t>&
     settings_.heuristic_preemption_callback();
   }
 
-  f_t gap              = upper_bound_ - lower_bound;
   f_t obj              = compute_user_objective(original_lp_, upper_bound_.load());
   f_t user_bound       = compute_user_objective(original_lp_, lower_bound);
+  f_t gap              = std::abs(obj - user_bound);
   f_t gap_rel          = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
   bool is_maximization = original_lp_.obj_scale < 0.0;
 
@@ -788,6 +785,7 @@ void branch_and_bound_t<i_t, f_t>::add_feasible_solution(f_t leaf_objective,
     incumbent_.set_incumbent_solution(leaf_objective, leaf_solution);
     upper_bound_ = leaf_objective;
     report(feasible_solution_symbol(thread_type), leaf_objective, get_lower_bound(), leaf_depth, 0);
+    root_reduced_cost_fixing(leaf_objective);
     send_solution = true;
   }
 
@@ -857,7 +855,7 @@ branch_variable_t<i_t> branch_and_bound_t<i_t, f_t>::variable_selection(
 
     case search_strategy_t::COEFFICIENT_DIVING:
       return coefficient_diving(
-        original_lp_, fractional, solution, var_up_locks_, var_down_locks_, log);
+        worker->leaf_problem, fractional, solution, var_up_locks_, var_down_locks_, log);
 
     case search_strategy_t::LINE_SEARCH_DIVING:
       return line_search_diving(fractional, solution, root_relax_soln_.x, log);
@@ -1126,7 +1124,7 @@ struct deterministic_diving_policy_t
       case search_strategy_t::COEFFICIENT_DIVING: {
         logger_t log;
         log.log = false;
-        return coefficient_diving<i_t, f_t>(this->bnb.original_lp_,
+        return coefficient_diving<i_t, f_t>(this->worker.leaf_problem,
                                             fractional,
                                             x,
                                             this->bnb.var_up_locks_,
@@ -1175,7 +1173,6 @@ std::pair<node_status_t, rounding_direction_t> branch_and_bound_t<i_t, f_t>::upd
   dual::status_t lp_status,
   Policy& policy)
 {
-  constexpr f_t inf                      = std::numeric_limits<f_t>::infinity();
   const f_t abs_fathom_tol               = settings_.absolute_mip_gap_tol / 10;
   lp_problem_t<i_t, f_t>& leaf_problem   = worker->leaf_problem;
   lp_solution_t<i_t, f_t>& leaf_solution = worker->leaf_solution;
@@ -1183,6 +1180,9 @@ std::pair<node_status_t, rounding_direction_t> branch_and_bound_t<i_t, f_t>::upd
   node_status_t status                   = node_status_t::PENDING;
   rounding_direction_t round_dir         = rounding_direction_t::NONE;
 
+  worker->recompute_basis  = true;
+  worker->recompute_bounds = true;
+
   if (lp_status == dual::status_t::DUAL_UNBOUNDED) {
     node_ptr->lower_bound = inf;
     policy.graphviz(search_tree, node_ptr, "infeasible", 0.0);
@@ -1242,6 +1242,8 @@ std::pair<node_status_t, rounding_direction_t> branch_and_bound_t<i_t, f_t>::upd
       assert(dir != rounding_direction_t::NONE);
 
       policy.update_objective_estimate(node_ptr, leaf_fractional, leaf_solution.x);
+      worker->recompute_basis  = false;
+      worker->recompute_bounds = false;
 
       logger_t log;
       log.log = false;
@@ -1376,6 +1378,31 @@ dual::status_t branch_and_bound_t<i_t, f_t>::solve_node_lp(
     node_ptr->vstatus[node_ptr->branch_var]);
 #endif
 
+  if (worker->start_bounds_updated) {
+    worker->start_bounds_updated = false;
+    worker->recompute_bounds     = true;
+
+    if (worker->search_strategy == BEST_FIRST) {
+      mutex_original_lp_.lock();
+      worker->start_lower = original_lp_.lower;
+      worker->start_upper = original_lp_.upper;
+      mutex_original_lp_.unlock();
+    } else {
+      // When diving, we are working on a separated subtree so we no longer can
+      // propagate the root bound changes from the root to the current node.
+      // Instead, we apply the reduced cost fixing directly
+      // to the starting bounds.
+      reduced_cost_fixing(root_relax_soln_.z,
+                          var_types_,
+                          settings_,
+                          root_objective_,
+                          upper_bound_.load(),
+                          worker->start_lower,
+                          worker->start_upper,
+                          worker->bounds_changed);
+    }
+  }
+
   bool feasible            = worker->set_lp_variable_bounds(node_ptr, settings_);
   dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED;
   worker->leaf_edge_norms  = edge_norms_;
@@ -1437,13 +1464,12 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(branch_and_bound_worker_t<i_t, f_
 
     f_t lower_bound = node_ptr->lower_bound;
     f_t upper_bound = upper_bound_;
-    f_t rel_gap     = user_relative_gap(original_lp_, upper_bound, lower_bound);
 
     // This is based on three assumptions:
     // - The stack only contains sibling nodes, i.e., the current node and it sibling (if
     // applicable)
-    // - The current node and its siblings uses the lower bound of the parent before solving the LP
-    // relaxation
+    // - The current node and its siblings uses the lower bound of the parent before solving the
+    // LP relaxation
     // - The lower bound of the parent is lower or equal to its children
     worker->lower_bound = lower_bound;
 
@@ -1481,10 +1507,6 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(branch_and_bound_worker_t<i_t, f_
 
     auto [node_status, round_dir] =
       update_tree(node_ptr, search_tree_, worker, lp_status, settings_.log);
-
-    worker->recompute_basis  = node_status != node_status_t::HAS_CHILDREN;
-    worker->recompute_bounds = node_status != node_status_t::HAS_CHILDREN;
-
     if (node_status == node_status_t::HAS_CHILDREN) {
       // The stack should only contain the children of the current parent.
       // If the stack size is greater than 0,
@@ -1554,7 +1576,6 @@ void branch_and_bound_t<i_t, f_t>::dive_with(branch_and_bound_worker_t<i_t, f_t>
 
     f_t lower_bound     = node_ptr->lower_bound;
     f_t upper_bound     = upper_bound_;
-    f_t rel_gap         = user_relative_gap(original_lp_, upper_bound, lower_bound);
     worker->lower_bound = lower_bound;
 
     if (node_ptr->lower_bound > upper_bound) {
@@ -1579,9 +1600,6 @@ void branch_and_bound_t<i_t, f_t>::dive_with(branch_and_bound_worker_t<i_t, f_t>
 
     auto [node_status, round_dir] = update_tree(node_ptr, dive_tree, worker, lp_status, log);
 
-    worker->recompute_basis  = node_status != node_status_t::HAS_CHILDREN;
-    worker->recompute_bounds = node_status != node_status_t::HAS_CHILDREN;
-
     if (node_status == node_status_t::HAS_CHILDREN) {
       if (round_dir == rounding_direction_t::UP) {
         stack.push_front(node_ptr->get_down_child());
@@ -1614,7 +1632,9 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
   std::array<i_t, num_search_strategies> max_num_workers_per_type =
     get_max_workers(num_workers, strategies);
 
+  mutex_original_lp_.lock();
   worker_pool_.init(num_workers, original_lp_, Arow_, var_types_, settings_);
+  mutex_original_lp_.unlock();
   active_workers_per_strategy_.fill(0);
 
 #ifdef CUOPT_LOG_DEBUG
@@ -1627,7 +1647,7 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
 #endif
 
   f_t lower_bound     = get_lower_bound();
-  f_t abs_gap         = upper_bound_ - lower_bound;
+  f_t abs_gap         = user_absolute_gap(original_lp_, upper_bound_.load(), lower_bound);
   f_t rel_gap         = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
   i_t last_node_depth = 0;
   i_t last_int_infeas = 0;
@@ -1637,7 +1657,7 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
          (active_workers_per_strategy_[0] > 0 || node_queue_.best_first_queue_size() > 0)) {
     bool launched_any_task = false;
     lower_bound            = get_lower_bound();
-    abs_gap                = upper_bound_ - lower_bound;
+    abs_gap                = user_absolute_gap(original_lp_, upper_bound_.load(), lower_bound);
     rel_gap                = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
 
     repair_heuristic_solutions();
@@ -1705,9 +1725,21 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
           continue;
         }
 
+        mutex_original_lp_.lock();
+        bool feasible = worker->init_best_first(start_node.value(), original_lp_);
+        mutex_original_lp_.unlock();
+
+        if (!feasible) {
+          // This node was put on the heap earlier but its variables bounds now violates the
+          // bounds at the root node
+          search_tree_.graphviz_node(
+            settings_.log, start_node.value(), "cutoff", start_node.value()->lower_bound);
+          search_tree_.update(start_node.value(), node_status_t::FATHOMED);
+          continue;
+        }
+
         // Remove the worker from the idle list.
         worker_pool_.pop_idle_worker();
-        worker->init_best_first(start_node.value(), original_lp_);
         last_node_depth = start_node.value()->depth;
         last_int_infeas = start_node.value()->integer_infeasible;
         active_workers_per_strategy_[strategy]++;
@@ -1725,8 +1757,10 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
           continue;
         }
 
+        mutex_original_lp_.lock();
         bool is_feasible =
           worker->init_diving(start_node.value(), strategy, original_lp_, settings_);
+        mutex_original_lp_.unlock();
         if (!is_feasible) { continue; }
 
         // Remove the worker from the idle list.
@@ -1752,14 +1786,14 @@ void branch_and_bound_t<i_t, f_t>::single_threaded_solve()
   branch_and_bound_worker_t<i_t, f_t> worker(0, original_lp_, Arow_, var_types_, settings_);
 
   f_t lower_bound = get_lower_bound();
-  f_t abs_gap     = upper_bound_ - lower_bound;
+  f_t abs_gap     = user_absolute_gap(original_lp_, upper_bound_.load(), lower_bound);
   f_t rel_gap     = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
 
   while (solver_status_ == mip_status_t::UNSET && abs_gap > settings_.absolute_mip_gap_tol &&
          rel_gap > settings_.relative_mip_gap_tol && node_queue_.best_first_queue_size() > 0) {
     bool launched_any_task = false;
     lower_bound            = get_lower_bound();
-    abs_gap                = upper_bound_ - lower_bound;
+    abs_gap                = user_absolute_gap(original_lp_, upper_bound_.load(), lower_bound);
     rel_gap                = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
 
     repair_heuristic_solutions();
@@ -1797,7 +1831,19 @@ void branch_and_bound_t<i_t, f_t>::single_threaded_solve()
       continue;
     }
 
-    worker.init_best_first(start_node.value(), original_lp_);
+    mutex_original_lp_.lock();
+    bool feasible = worker.init_best_first(start_node.value(), original_lp_);
+    mutex_original_lp_.unlock();
+
+    if (!feasible) {
+      // This node was put on the heap earlier but its variables bounds now violates the
+      // bounds at the root node
+      search_tree_.graphviz_node(
+        settings_.log, start_node.value(), "cutoff", start_node.value()->lower_bound);
+      search_tree_.update(start_node.value(), node_status_t::FATHOMED);
+      continue;
+    }
+
     plunge_with(&worker);
   }
 }
@@ -2275,32 +2321,46 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
         return mip_status_t::NUMERICAL;
       }
 
+      std::vector<bool> bounds_changed;
+      std::vector<f_t> new_lower;
+      std::vector<f_t> new_upper;
+
       if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) {
+        mutex_original_lp_.lock();
+        new_lower = original_lp_.lower;
+        new_upper = original_lp_.upper;
+        mutex_original_lp_.unlock();
         mutex_upper_.lock();
         last_upper_bound = upper_bound_.load();
-        std::vector<f_t> lower_bounds;
-        std::vector<f_t> upper_bounds;
-        find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds);
+        reduced_cost_fixing(root_relax_soln_.z,
+                            var_types_,
+                            settings_,
+                            root_objective_,
+                            upper_bound_.load(),
+                            new_lower,
+                            new_upper,
+                            bounds_changed);
         mutex_upper_.unlock();
         mutex_original_lp_.lock();
-        original_lp_.lower = lower_bounds;
-        original_lp_.upper = upper_bounds;
+        original_lp_.lower = new_lower;
+        original_lp_.upper = new_upper;
         mutex_original_lp_.unlock();
       }
 
-      // Try to do bound strengthening
-      std::vector<bool> bounds_changed(original_lp_.num_cols, true);
-      std::vector<char> row_sense;
 #ifdef CHECK_MATRICES
       settings_.log.printf("Before A check\n");
       original_lp_.A.check_matrix();
 #endif
       original_lp_.A.to_compressed_row(Arow_);
 
+      // Try to do bound strengthening
+      std::vector<char> row_sense;
       f_t node_presolve_start_time = tic();
       bounds_strengthening_t<i_t, f_t> node_presolve(original_lp_, Arow_, row_sense, var_types_);
-      std::vector<f_t> new_lower = original_lp_.lower;
-      std::vector<f_t> new_upper = original_lp_.upper;
+      bounds_changed.assign(original_lp_.num_cols, true);
+      new_lower = original_lp_.lower;
+      new_upper = original_lp_.upper;
+
       bool feasible =
         node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper);
       mutex_original_lp_.lock();
@@ -2410,7 +2470,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
       report(' ', obj, root_objective_, 0, num_fractional);
 
       f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_);
-      f_t abs_gap = upper_bound_.load() - root_objective_;
+      f_t abs_gap = user_absolute_gap(original_lp_, upper_bound_.load(), root_objective_);
       if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) {
         set_solution_at_root(solution, cut_info);
         set_final_solution(solution, root_objective_);
@@ -2468,13 +2528,26 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   }
 
   if (settings_.reduced_cost_strengthening >= 2 && upper_bound_.load() < last_upper_bound) {
-    std::vector<f_t> lower_bounds;
-    std::vector<f_t> upper_bounds;
-    i_t num_fixed = find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds);
-    if (num_fixed > 0) {
-      std::vector<bool> bounds_changed(original_lp_.num_cols, true);
-      std::vector<char> row_sense;
+    std::vector<bool> bounds_changed;
 
+    mutex_original_lp_.lock();
+    std::vector<f_t> lower_bounds = original_lp_.lower;
+    std::vector<f_t> upper_bounds = original_lp_.upper;
+    mutex_original_lp_.unlock();
+
+    mutex_upper_.lock();
+    auto [num_fixed, num_improved] = reduced_cost_fixing(root_relax_soln_.z,
+                                                         var_types_,
+                                                         settings_,
+                                                         root_objective_,
+                                                         upper_bound_.load(),
+                                                         lower_bounds,
+                                                         upper_bounds,
+                                                         bounds_changed);
+    mutex_upper_.unlock();
+
+    if (num_fixed > 0 || num_improved > 0) {
+      std::vector<char> row_sense;
       bounds_strengthening_t<i_t, f_t> node_presolve(original_lp_, Arow_, row_sense, var_types_);
 
       mutex_original_lp_.lock();
@@ -2483,6 +2556,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
       bool feasible      = node_presolve.bounds_strengthening(
         settings_, bounds_changed, original_lp_.lower, original_lp_.upper);
       mutex_original_lp_.unlock();
+
       if (!feasible) {
         settings_.log.printf("Bound strengthening failed\n");
         return mip_status_t::NUMERICAL;  // We had a feasible integer solution, but bound
@@ -2601,6 +2675,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
 //  Deterministic implementation
 // ============================================================================
 
+// TODO: reduced cost fixing is not enabled for the deterministic mode
 // The deterministic BSP model is based on letting independent workers execute during virtual time
 // intervals, and exchange data during serialized interval sync points.
 /*
@@ -3011,7 +3086,7 @@ void branch_and_bound_t<i_t, f_t>::deterministic_sync_callback()
 
   f_t lower_bound = deterministic_compute_lower_bound();
   f_t upper_bound = upper_bound_.load();
-  f_t abs_gap     = upper_bound - lower_bound;
+  f_t abs_gap     = user_absolute_gap(original_lp_, upper_bound, lower_bound);
   f_t rel_gap     = user_relative_gap(original_lp_, upper_bound, lower_bound);
 
   if (abs_gap <= settings_.absolute_mip_gap_tol || rel_gap <= settings_.relative_mip_gap_tol) {
@@ -3083,17 +3158,30 @@ node_status_t branch_and_bound_t<i_t, f_t>::solve_node_deterministic(
   raft::common::nvtx::range scope("BB::solve_node_deterministic");
 
   double work_units_at_start = worker.work_context.global_work_units_elapsed;
+  bool feasible              = true;
 
   std::fill(worker.bounds_changed.begin(), worker.bounds_changed.end(), false);
 
+  // FIXME: Add mutex protection for original_lp_ bounds access in solve_node_deterministic.
+  // The code reads original_lp_.lower and original_lp_.upper without mutex protection at lines
+  // 3171-3172 and 3177-3178. However, update_root_bounds (called from set_new_solution during
+  // solution processing) modifies these same bounds while holding mutex_original_lp_. This creates
+  // a data race: workers can read bounds while the main thread updates them via heuristic
+  // solutions. Acquire mutex_original_lp_ before accessing the bounds, or take a snapshot at the
+  // start of solve_node_deterministic to avoid repeated locking.
+  // (This fix is needed for deterministic mode + reduced cost strengthening)
   if (worker.recompute_bounds_and_basis) {
-    worker.leaf_problem.lower = original_lp_.lower;
-    worker.leaf_problem.upper = original_lp_.upper;
-    node_ptr->get_variable_bounds(
-      worker.leaf_problem.lower, worker.leaf_problem.upper, worker.bounds_changed);
+    feasible = node_ptr->get_variable_bounds(original_lp_.lower,
+                                             original_lp_.upper,
+                                             worker.leaf_problem.lower,
+                                             worker.leaf_problem.upper,
+                                             worker.bounds_changed);
   } else {
-    node_ptr->update_branched_variable_bounds(
-      worker.leaf_problem.lower, worker.leaf_problem.upper, worker.bounds_changed);
+    feasible = node_ptr->update_branched_variable_bounds(original_lp_.lower,
+                                                         original_lp_.upper,
+                                                         worker.leaf_problem.lower,
+                                                         worker.leaf_problem.upper,
+                                                         worker.bounds_changed);
   }
 
   double remaining_time = settings_.time_limit - toc(exploration_stats_.start_time);
@@ -3107,11 +3195,12 @@ node_status_t branch_and_bound_t<i_t, f_t>::solve_node_deterministic(
   lp_settings.time_limit    = remaining_time;
   lp_settings.scale_columns = false;
 
-  bool feasible = true;
 #ifndef DETERMINISM_DISABLE_BOUNDS_STRENGTHENING
   raft::common::nvtx::range scope_bs("BB::bound_strengthening");
-  feasible = worker.node_presolver.bounds_strengthening(
-    lp_settings, worker.bounds_changed, worker.leaf_problem.lower, worker.leaf_problem.upper);
+  if (feasible) {
+    feasible = worker.node_presolver.bounds_strengthening(
+      lp_settings, worker.bounds_changed, worker.leaf_problem.lower, worker.leaf_problem.upper);
+  }
 
   if (settings_.deterministic) {
     // TEMP APPROXIMATION;
@@ -3696,15 +3785,20 @@ void branch_and_bound_t<i_t, f_t>::deterministic_dive(
 
     // Setup bounds for this node
     std::fill(worker.bounds_changed.begin(), worker.bounds_changed.end(), false);
+    bool feasible = true;
 
     if (worker.recompute_bounds_and_basis) {
-      worker.leaf_problem.lower = worker.dive_lower;
-      worker.leaf_problem.upper = worker.dive_upper;
-      node_ptr->get_variable_bounds(
-        worker.leaf_problem.lower, worker.leaf_problem.upper, worker.bounds_changed);
+      feasible = node_ptr->get_variable_bounds(worker.dive_lower,
+                                               worker.dive_upper,
+                                               worker.leaf_problem.lower,
+                                               worker.leaf_problem.upper,
+                                               worker.bounds_changed);
     } else {
-      node_ptr->update_branched_variable_bounds(
-        worker.leaf_problem.lower, worker.leaf_problem.upper, worker.bounds_changed);
+      feasible = node_ptr->update_branched_variable_bounds(worker.dive_lower,
+                                                           worker.dive_upper,
+                                                           worker.leaf_problem.lower,
+                                                           worker.leaf_problem.upper,
+                                                           worker.bounds_changed);
     }
 
     double remaining_time = settings_.time_limit - toc(exploration_stats_.start_time);
@@ -3719,8 +3813,10 @@ void branch_and_bound_t<i_t, f_t>::deterministic_dive(
     lp_settings.scale_columns = false;
 
 #ifndef DETERMINISM_DISABLE_BOUNDS_STRENGTHENING
-    bool feasible = worker.node_presolver.bounds_strengthening(
-      lp_settings, worker.bounds_changed, worker.leaf_problem.lower, worker.leaf_problem.upper);
+    if (feasible) {
+      feasible = worker.node_presolver.bounds_strengthening(
+        lp_settings, worker.bounds_changed, worker.leaf_problem.lower, worker.leaf_problem.upper);
+    }
 
     if (settings_.deterministic) {
       // TEMP APPROXIMATION;
diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp
index 98674b7f9e..da4ae53840 100644
--- a/cpp/src/branch_and_bound/branch_and_bound.hpp
+++ b/cpp/src/branch_and_bound/branch_and_bound.hpp
@@ -8,12 +8,12 @@
 #pragma once
 
 #include <branch_and_bound/bb_event.hpp>
-#include <branch_and_bound/branch_and_bound_worker.hpp>
 #include <branch_and_bound/deterministic_workers.hpp>
-#include <branch_and_bound/diving_heuristics.hpp>
 #include <branch_and_bound/mip_node.hpp>
 #include <branch_and_bound/node_queue.hpp>
 #include <branch_and_bound/pseudo_costs.hpp>
+#include <branch_and_bound/worker.hpp>
+#include <branch_and_bound/worker_pool.hpp>
 
 #include <cuts/cuts.hpp>
 
@@ -30,8 +30,6 @@
 #include <utilities/work_limit_context.hpp>
 #include <utilities/work_unit_scheduler.hpp>
 
-#include <omp.h>
-
 #include <atomic>
 #include <functional>
 #include <future>
@@ -102,7 +100,7 @@ class branch_and_bound_t {
     }
   }
 
-  // Set a solution based on the user problem during the course of the solve
+  // Set a solution based on the user problem during solve time
   void set_new_solution(const std::vector<f_t>& solution);
 
   // This queues the solution to be processed at the correct work unit timestamp
@@ -133,10 +131,6 @@ class branch_and_bound_t {
                                     std::vector<i_t>& nonbasic_list,
                                     std::vector<f_t>& edge_norms);
 
-  i_t find_reduced_cost_fixings(f_t upper_bound,
-                                std::vector<f_t>& lower_bounds,
-                                std::vector<f_t>& upper_bounds);
-
   // The main entry routine. Returns the solver status and populates solution with the incumbent.
   mip_status_t solve(mip_solution_t<i_t, f_t>& solution);
 
@@ -307,6 +301,8 @@ class branch_and_bound_t {
     dual::status_t lp_status,
     logger_t& log);
 
+  void root_reduced_cost_fixing(f_t upper_bound);
+
   // ============================================================================
   // Deterministic BSP (Bulk Synchronous Parallel) methods for deterministic parallel B&B
   // ============================================================================
diff --git a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp
deleted file mode 100644
index 4de2b43cae..0000000000
--- a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp
+++ /dev/null
@@ -1,281 +0,0 @@
-/* clang-format off */
-/*
- * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
- * SPDX-License-Identifier: Apache-2.0
- */
-/* clang-format on */
-
-#pragma once
-
-#include <branch_and_bound/mip_node.hpp>
-
-#include <dual_simplex/basis_updates.hpp>
-#include <dual_simplex/bounds_strengthening.hpp>
-#include <dual_simplex/phase2.hpp>
-
-#include <utilities/pcgenerator.hpp>
-
-#include <array>
-#include <deque>
-#include <mutex>
-#include <vector>
-
-namespace cuopt::linear_programming::dual_simplex {
-
-constexpr int num_search_strategies = 5;
-
-// Indicate the search and variable selection algorithms used by each thread
-// in B&B (See [1]).
-//
-// [1] T. Achterberg, “Constraint Integer Programming,” PhD, Technischen Universität Berlin,
-// Berlin, 2007. doi: 10.14279/depositonce-1634.
-enum search_strategy_t : int {
-  BEST_FIRST         = 0,  // Best-First + Plunging.
-  PSEUDOCOST_DIVING  = 1,  // Pseudocost diving (9.2.5)
-  LINE_SEARCH_DIVING = 2,  // Line search diving (9.2.4)
-  GUIDED_DIVING      = 3,  // Guided diving (9.2.3).
-  COEFFICIENT_DIVING = 4   // Coefficient diving (9.2.1)
-};
-
-template <typename i_t, typename f_t>
-struct branch_and_bound_stats_t {
-  f_t start_time                         = 0.0;
-  omp_atomic_t<f_t> total_lp_solve_time  = 0.0;
-  omp_atomic_t<int64_t> nodes_explored   = 0;
-  omp_atomic_t<int64_t> nodes_unexplored = 0;
-  omp_atomic_t<int64_t> total_lp_iters   = 0;
-  omp_atomic_t<i_t> nodes_since_last_log = 0;
-  omp_atomic_t<f_t> last_log             = 0.0;
-};
-
-template <typename i_t, typename f_t>
-class branch_and_bound_worker_t {
- public:
-  const i_t worker_id;
-  omp_atomic_t<search_strategy_t> search_strategy;
-  omp_atomic_t<bool> is_active;
-  omp_atomic_t<f_t> lower_bound;
-
-  lp_problem_t<i_t, f_t> leaf_problem;
-  lp_solution_t<i_t, f_t> leaf_solution;
-  std::vector<f_t> leaf_edge_norms;
-
-  basis_update_mpf_t<i_t, f_t> basis_factors;
-  std::vector<i_t> basic_list;
-  std::vector<i_t> nonbasic_list;
-
-  bounds_strengthening_t<i_t, f_t> node_presolver;
-  std::vector<bool> bounds_changed;
-
-  std::vector<f_t> start_lower;
-  std::vector<f_t> start_upper;
-  mip_node_t<i_t, f_t>* start_node;
-
-  pcgenerator_t rng;
-
-  bool recompute_basis  = true;
-  bool recompute_bounds = true;
-
-  branch_and_bound_worker_t(i_t worker_id,
-                            const lp_problem_t<i_t, f_t>& original_lp,
-                            const csr_matrix_t<i_t, f_t>& Arow,
-                            const std::vector<variable_type_t>& var_type,
-                            const simplex_solver_settings_t<i_t, f_t>& settings)
-    : worker_id(worker_id),
-      search_strategy(BEST_FIRST),
-      is_active(false),
-      lower_bound(-std::numeric_limits<f_t>::infinity()),
-      leaf_problem(original_lp),
-      leaf_solution(original_lp.num_rows, original_lp.num_cols),
-      basis_factors(original_lp.num_rows, settings.refactor_frequency),
-      basic_list(original_lp.num_rows),
-      nonbasic_list(),
-      node_presolver(leaf_problem, Arow, {}, var_type),
-      bounds_changed(original_lp.num_cols, false),
-      rng(settings.random_seed + pcgenerator_t::default_seed + worker_id,
-          pcgenerator_t::default_stream ^ worker_id)
-  {
-  }
-
-  // Set the `start_node` for best-first search.
-  void init_best_first(mip_node_t<i_t, f_t>* node, const lp_problem_t<i_t, f_t>& original_lp)
-  {
-    start_node      = node;
-    start_lower     = original_lp.lower;
-    start_upper     = original_lp.upper;
-    search_strategy = BEST_FIRST;
-    lower_bound     = node->lower_bound;
-    is_active       = true;
-  }
-
-  // Initialize the worker for diving, setting the `start_node`, `start_lower` and
-  // `start_upper`. Returns `true` if the starting node is feasible via
-  // bounds propagation.
-  bool init_diving(mip_node_t<i_t, f_t>* node,
-                   search_strategy_t type,
-                   const lp_problem_t<i_t, f_t>& original_lp,
-                   const simplex_solver_settings_t<i_t, f_t>& settings)
-  {
-    internal_node = node->detach_copy();
-    start_node    = &internal_node;
-
-    start_lower     = original_lp.lower;
-    start_upper     = original_lp.upper;
-    search_strategy = type;
-    lower_bound     = node->lower_bound;
-    is_active       = true;
-
-    std::fill(bounds_changed.begin(), bounds_changed.end(), false);
-    node->get_variable_bounds(start_lower, start_upper, bounds_changed);
-    return node_presolver.bounds_strengthening(settings, bounds_changed, start_lower, start_upper);
-  }
-
-  // Set the variables bounds for the LP relaxation of the current node.
-  bool set_lp_variable_bounds(mip_node_t<i_t, f_t>* node_ptr,
-                              const simplex_solver_settings_t<i_t, f_t>& settings)
-  {
-    // Reset the bound_changed markers
-    std::fill(bounds_changed.begin(), bounds_changed.end(), false);
-
-    // Set the correct bounds for the leaf problem
-    if (recompute_bounds) {
-      leaf_problem.lower = start_lower;
-      leaf_problem.upper = start_upper;
-      node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed);
-
-    } else {
-      node_ptr->update_branched_variable_bounds(
-        leaf_problem.lower, leaf_problem.upper, bounds_changed);
-    }
-
-    return node_presolver.bounds_strengthening(
-      settings, bounds_changed, leaf_problem.lower, leaf_problem.upper);
-  }
-
- private:
-  // For diving, we need to store the full node instead of
-  // of just a pointer, since it is not stored in the tree anymore.
-  // To keep the same interface across all worker types,
-  // this will be used as a temporary storage and
-  // will be pointed by `start_node`.
-  // For exploration, this will not be used.
-  mip_node_t<i_t, f_t> internal_node;
-};
-
-template <typename i_t, typename f_t>
-class branch_and_bound_worker_pool_t {
- public:
-  void init(i_t num_workers,
-            const lp_problem_t<i_t, f_t>& original_lp,
-            const csr_matrix_t<i_t, f_t>& Arow,
-            const std::vector<variable_type_t>& var_type,
-            const simplex_solver_settings_t<i_t, f_t>& settings)
-  {
-    workers_.resize(num_workers);
-    num_idle_workers_ = num_workers;
-    for (i_t i = 0; i < num_workers; ++i) {
-      workers_[i] = std::make_unique<branch_and_bound_worker_t<i_t, f_t>>(
-        i, original_lp, Arow, var_type, settings);
-      idle_workers_.push_front(i);
-    }
-
-    is_initialized = true;
-  }
-
-  // Here, we are assuming that the scheduler is the only
-  // thread that can retrieve/pop an idle worker.
-  branch_and_bound_worker_t<i_t, f_t>* get_idle_worker()
-  {
-    std::lock_guard<omp_mutex_t> lock(mutex_);
-    if (idle_workers_.empty()) {
-      return nullptr;
-    } else {
-      i_t idx = idle_workers_.front();
-      return workers_[idx].get();
-    }
-  }
-
-  // Here, we are assuming that the scheduler is the only
-  // thread that can retrieve/pop an idle worker.
-  void pop_idle_worker()
-  {
-    std::lock_guard<omp_mutex_t> lock(mutex_);
-    if (!idle_workers_.empty()) {
-      idle_workers_.pop_front();
-      num_idle_workers_--;
-    }
-  }
-
-  void return_worker_to_pool(branch_and_bound_worker_t<i_t, f_t>* worker)
-  {
-    worker->is_active = false;
-    std::lock_guard<omp_mutex_t> lock(mutex_);
-    idle_workers_.push_back(worker->worker_id);
-    num_idle_workers_++;
-  }
-
-  f_t get_lower_bound()
-  {
-    f_t lower_bound = std::numeric_limits<f_t>::infinity();
-
-    if (is_initialized) {
-      for (i_t i = 0; i < workers_.size(); ++i) {
-        if (workers_[i]->search_strategy == BEST_FIRST && workers_[i]->is_active) {
-          lower_bound = std::min(workers_[i]->lower_bound.load(), lower_bound);
-        }
-      }
-    }
-
-    return lower_bound;
-  }
-
-  i_t num_idle_workers() { return num_idle_workers_; }
-
- private:
-  // Worker pool
-  std::vector<std::unique_ptr<branch_and_bound_worker_t<i_t, f_t>>> workers_;
-  bool is_initialized = false;
-
-  omp_mutex_t mutex_;
-  std::deque<i_t> idle_workers_;
-  omp_atomic_t<i_t> num_idle_workers_;
-};
-
-template <typename f_t, typename i_t>
-std::vector<search_strategy_t> get_search_strategies(
-  diving_heuristics_settings_t<i_t, f_t> settings)
-{
-  std::vector<search_strategy_t> types;
-  types.reserve(num_search_strategies);
-  types.push_back(BEST_FIRST);
-  if (settings.pseudocost_diving != 0) { types.push_back(PSEUDOCOST_DIVING); }
-  if (settings.line_search_diving != 0) { types.push_back(LINE_SEARCH_DIVING); }
-  if (settings.guided_diving != 0) { types.push_back(GUIDED_DIVING); }
-  if (settings.coefficient_diving != 0) { types.push_back(COEFFICIENT_DIVING); }
-  return types;
-}
-
-template <typename i_t>
-std::array<i_t, num_search_strategies> get_max_workers(
-  i_t num_workers, const std::vector<search_strategy_t>& strategies)
-{
-  std::array<i_t, num_search_strategies> max_num_workers;
-  max_num_workers.fill(0);
-
-  i_t bfs_workers             = std::max(strategies.size() == 1 ? num_workers : num_workers / 4, 1);
-  max_num_workers[BEST_FIRST] = bfs_workers;
-
-  i_t diving_workers = (num_workers - bfs_workers);
-  i_t m              = strategies.size() - 1;
-
-  for (size_t i = 1, k = 0; i < strategies.size(); ++i) {
-    i_t start                      = (double)k * diving_workers / m;
-    i_t end                        = (double)(k + 1) * diving_workers / m;
-    max_num_workers[strategies[i]] = end - start;
-    ++k;
-  }
-
-  return max_num_workers;
-}
-
-}  // namespace cuopt::linear_programming::dual_simplex
diff --git a/cpp/src/branch_and_bound/deterministic_workers.hpp b/cpp/src/branch_and_bound/deterministic_workers.hpp
index 7a074051c6..50b8e905db 100644
--- a/cpp/src/branch_and_bound/deterministic_workers.hpp
+++ b/cpp/src/branch_and_bound/deterministic_workers.hpp
@@ -8,9 +8,9 @@
 #pragma once
 
 #include <branch_and_bound/bb_event.hpp>
-#include <branch_and_bound/branch_and_bound_worker.hpp>
 #include <branch_and_bound/diving_heuristics.hpp>
 #include <branch_and_bound/node_queue.hpp>
+#include <branch_and_bound/worker.hpp>
 
 #include <utilities/work_limit_context.hpp>
 
@@ -316,10 +316,13 @@ class deterministic_diving_worker_t
   void enqueue_dive_node(mip_node_t<i_t, f_t>* node, const lp_problem_t<i_t, f_t>& original_lp)
   {
     dive_queue_entry_t<i_t, f_t> entry;
-    entry.resolved_lower = original_lp.lower;
-    entry.resolved_upper = original_lp.upper;
     std::vector<bool> bounds_changed(original_lp.num_cols, false);
-    node->get_variable_bounds(entry.resolved_lower, entry.resolved_upper, bounds_changed);
+    bool feasible = node->get_variable_bounds(original_lp.lower,
+                                              original_lp.upper,
+                                              entry.resolved_lower,
+                                              entry.resolved_upper,
+                                              bounds_changed);
+    if (!feasible) { return; }
     entry.node = node->detach_copy();
     dive_queue.push_back(std::move(entry));
   }
diff --git a/cpp/src/branch_and_bound/mip_node.hpp b/cpp/src/branch_and_bound/mip_node.hpp
index 58e1a6f683..ca41e73653 100644
--- a/cpp/src/branch_and_bound/mip_node.hpp
+++ b/cpp/src/branch_and_bound/mip_node.hpp
@@ -98,22 +98,62 @@ class mip_node_t {
     children[1]      = nullptr;
   }
 
-  void get_variable_bounds(std::vector<f_t>& lower,
+  // Check if no node in the path violates the initial bounds after branching on a given variable.
+  // The bound violation can happen after the initial bounds is changed via reduced cost
+  // strengthening and one of the node in the path was created based on the old values.
+  bool check_variable_bounds(const std::vector<f_t>& start_lower,
+                             const std::vector<f_t>& start_upper)
+  {
+    if (branch_var_upper < start_lower[branch_var] || branch_var_lower > start_upper[branch_var]) {
+      return false;
+    }
+
+    mip_node_t* parent_ptr = parent;
+    while (parent_ptr != nullptr && parent_ptr->node_id != 0) {
+      if (parent_ptr->branch_var_upper < start_lower[parent_ptr->branch_var] ||
+          parent_ptr->branch_var_lower > start_upper[parent_ptr->branch_var]) {
+        return false;
+      }
+      parent_ptr = parent_ptr->parent;
+    }
+
+    return true;
+  }
+
+  // Get the variable bounds starting at the current node and then traversing it back until the
+  // the root node. The bounds are initially set based on the `start_lower` and `start_upper`.
+  // Return true if all bounds are valid (i.e., no node in the path violates the initial bounds
+  // after branching on a given variable).
+  bool get_variable_bounds(const std::vector<f_t>& start_lower,
+                           const std::vector<f_t>& start_upper,
+                           std::vector<f_t>& lower,
                            std::vector<f_t>& upper,
                            std::vector<bool>& bounds_changed) const
   {
-    update_branched_variable_bounds(lower, upper, bounds_changed);
+    lower = start_lower;
+    upper = start_upper;
+
+    bool feasible =
+      update_branched_variable_bounds(start_lower, start_upper, lower, upper, bounds_changed);
+    if (!feasible) { return false; }
 
-    mip_node_t<i_t, f_t>* parent_ptr = parent;
+    mip_node_t* parent_ptr = parent;
     while (parent_ptr != nullptr && parent_ptr->node_id != 0) {
-      parent_ptr->update_branched_variable_bounds(lower, upper, bounds_changed);
+      feasible = parent_ptr->update_branched_variable_bounds(
+        start_lower, start_upper, lower, upper, bounds_changed);
+      if (!feasible) { return false; }
       parent_ptr = parent_ptr->parent;
     }
+
+    return true;
   }
 
   // Here we assume that we are traversing from the deepest node to the
-  // root of the tree
-  void update_branched_variable_bounds(std::vector<f_t>& lower,
+  // root of the tree. Return true if no bounds were violated in this node
+  // considering the starting bounds
+  bool update_branched_variable_bounds(const std::vector<f_t>& start_lower,
+                                       const std::vector<f_t>& start_upper,
+                                       std::vector<f_t>& lower,
                                        std::vector<f_t>& upper,
                                        std::vector<bool>& bounds_changed) const
   {
@@ -122,15 +162,22 @@ class mip_node_t {
     assert(upper.size() > branch_var);
     assert(bounds_changed.size() > branch_var);
 
+    // If the start bounds has changed (via reduced cost strengthening), check if the
+    // bounds in the node is still valid.
+    if (branch_var_upper < start_lower[branch_var] || branch_var_lower > start_upper[branch_var]) {
+      return false;
+    }
+
     // If the bounds have already been updated on another node,
     // skip this node as it contains looser bounds, since we
     // are traversing up the tree toward the root
-    if (bounds_changed[branch_var]) { return; }
+    if (bounds_changed[branch_var]) { return true; }
 
     // Apply the bounds at the current node
     lower[branch_var]          = branch_var_lower;
     upper[branch_var]          = branch_var_upper;
     bounds_changed[branch_var] = true;
+    return true;
   }
 
   mip_node_t* get_down_child() const { return children[0].get(); }
diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp
index 6b6c6917b6..08970c0fcf 100644
--- a/cpp/src/branch_and_bound/pseudo_costs.hpp
+++ b/cpp/src/branch_and_bound/pseudo_costs.hpp
@@ -7,8 +7,8 @@
 
 #pragma once
 
-#include <branch_and_bound/branch_and_bound_worker.hpp>
 #include <branch_and_bound/mip_node.hpp>
+#include <branch_and_bound/worker.hpp>
 
 #include <dual_simplex/basis_updates.hpp>
 #include <dual_simplex/logger.hpp>
diff --git a/cpp/src/branch_and_bound/reduced_cost_fixing.hpp b/cpp/src/branch_and_bound/reduced_cost_fixing.hpp
new file mode 100644
index 0000000000..53fd1737c4
--- /dev/null
+++ b/cpp/src/branch_and_bound/reduced_cost_fixing.hpp
@@ -0,0 +1,79 @@
+/* clang-format off */
+/*
+ * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-License-Identifier: Apache-2.0
+ */
+/* clang-format on */
+
+#pragma once
+
+#include <dual_simplex/presolve.hpp>
+#include <vector>
+
+namespace cuopt::linear_programming::dual_simplex {
+
+// Applies reduced cost fixing over the lower and bounds. Stores the bounds changes
+// for applying bound strengthening later. Returns {num_fixed, num_improved}.
+template <typename i_t, typename f_t>
+std::pair<i_t, i_t> reduced_cost_fixing(const std::vector<f_t>& reduced_costs,
+                                        const std::vector<variable_type_t>& var_types,
+                                        const simplex_solver_settings_t<i_t, f_t>& settings,
+                                        f_t obj,
+                                        f_t upper_bound,
+                                        std::vector<f_t>& lower_bounds,
+                                        std::vector<f_t>& upper_bounds,
+                                        std::vector<bool>& bounds_changed)
+{
+  const f_t threshold   = 100.0 * settings.integer_tol;
+  const f_t weaken      = settings.integer_tol;
+  const f_t fixed_tol   = settings.fixed_tol;
+  i_t num_improved      = 0;
+  i_t num_fixed         = 0;
+  i_t num_cols_to_check = reduced_costs.size();  // Reduced costs will be smaller than the original
+                                                 // problem because we have added slacks for cuts
+
+  bounds_changed.assign(lower_bounds.size(), false);
+
+  for (i_t j = 0; j < num_cols_to_check; j++) {
+    if (std::isfinite(reduced_costs[j]) && std::abs(reduced_costs[j]) > threshold) {
+      const f_t lower_j            = lower_bounds[j];
+      const f_t upper_j            = upper_bounds[j];
+      const f_t abs_gap            = upper_bound - obj;
+      f_t reduced_cost_upper_bound = upper_j;
+      f_t reduced_cost_lower_bound = lower_j;
+      if (lower_j > -inf && reduced_costs[j] > 0) {
+        const f_t new_upper_bound = lower_j + abs_gap / reduced_costs[j];
+        reduced_cost_upper_bound  = var_types[j] == variable_type_t::INTEGER
+                                      ? std::floor(new_upper_bound + weaken)
+                                      : new_upper_bound;
+        if (reduced_cost_upper_bound < upper_j && var_types[j] == variable_type_t::INTEGER) {
+          ++num_improved;
+          upper_bounds[j]   = reduced_cost_upper_bound;
+          bounds_changed[j] = true;
+        }
+      }
+      if (upper_j < inf && reduced_costs[j] < 0) {
+        const f_t new_lower_bound = upper_j + abs_gap / reduced_costs[j];
+        reduced_cost_lower_bound  = var_types[j] == variable_type_t::INTEGER
+                                      ? std::ceil(new_lower_bound - weaken)
+                                      : new_lower_bound;
+        if (reduced_cost_lower_bound > lower_j && var_types[j] == variable_type_t::INTEGER) {
+          ++num_improved;
+          lower_bounds[j]   = reduced_cost_lower_bound;
+          bounds_changed[j] = true;
+        }
+      }
+      if (var_types[j] == variable_type_t::INTEGER &&
+          reduced_cost_upper_bound <= reduced_cost_lower_bound + fixed_tol) {
+        ++num_fixed;
+      }
+    }
+  }
+
+  if (num_fixed > 0 || num_improved > 0) {
+    settings.log.printf(
+      "Reduced costs: Found %d improved bounds and %d fixed variables\n", num_improved, num_fixed);
+  }
+  return {num_fixed, num_improved};
+}
+}  // namespace cuopt::linear_programming::dual_simplex
diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp
new file mode 100644
index 0000000000..53ef5ee541
--- /dev/null
+++ b/cpp/src/branch_and_bound/worker.hpp
@@ -0,0 +1,191 @@
+/* clang-format off */
+/*
+ * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-License-Identifier: Apache-2.0
+ */
+/* clang-format on */
+
+#pragma once
+
+#include <branch_and_bound/mip_node.hpp>
+
+#include <dual_simplex/basis_updates.hpp>
+#include <dual_simplex/bounds_strengthening.hpp>
+#include <dual_simplex/phase2.hpp>
+
+#include <utilities/pcgenerator.hpp>
+
+#include <array>
+#include <deque>
+#include <mutex>
+#include <vector>
+
+namespace cuopt::linear_programming::dual_simplex {
+
+constexpr int num_search_strategies = 5;
+
+// Indicate the search and variable selection algorithms used by each thread
+// in B&B (See [1]).
+//
+// [1] T. Achterberg, “Constraint Integer Programming,” PhD, Technischen Universität Berlin,
+// Berlin, 2007. doi: 10.14279/depositonce-1634.
+enum search_strategy_t : int {
+  BEST_FIRST         = 0,  // Best-First + Plunging.
+  PSEUDOCOST_DIVING  = 1,  // Pseudocost diving (9.2.5)
+  LINE_SEARCH_DIVING = 2,  // Line search diving (9.2.4)
+  GUIDED_DIVING      = 3,  // Guided diving (9.2.3).
+  COEFFICIENT_DIVING = 4   // Coefficient diving (9.2.1)
+};
+
+template <typename i_t, typename f_t>
+struct branch_and_bound_stats_t {
+  f_t start_time                         = 0.0;
+  omp_atomic_t<f_t> total_lp_solve_time  = 0.0;
+  omp_atomic_t<int64_t> nodes_explored   = 0;
+  omp_atomic_t<int64_t> nodes_unexplored = 0;
+  omp_atomic_t<int64_t> total_lp_iters   = 0;
+  omp_atomic_t<i_t> nodes_since_last_log = 0;
+  omp_atomic_t<f_t> last_log             = 0.0;
+};
+
+template <typename i_t, typename f_t>
+class branch_and_bound_worker_t {
+ public:
+  const i_t worker_id;
+  omp_atomic_t<search_strategy_t> search_strategy;
+  omp_atomic_t<bool> is_active;
+  omp_atomic_t<f_t> lower_bound;
+
+  lp_problem_t<i_t, f_t> leaf_problem;
+  lp_solution_t<i_t, f_t> leaf_solution;
+  std::vector<f_t> leaf_edge_norms;
+
+  basis_update_mpf_t<i_t, f_t> basis_factors;
+  std::vector<i_t> basic_list;
+  std::vector<i_t> nonbasic_list;
+
+  bounds_strengthening_t<i_t, f_t> node_presolver;
+  std::vector<bool> bounds_changed;
+
+  // We need to maintain a copy in each worker for 2 reasons:
+  //
+  // - The root LP may be modified by multiple threads and
+  // require a mutex for accessing its variable bounds.
+  // Since we are maintain a copy here, we can access
+  // without having to acquire the mutex. Only if the
+  // bounds is modified, then we acquire the lock and update the copy.
+  //
+  // - When diving, we are working on a separated subtree. Hence, we cannot
+  // retrieve the bounds from the main tree. Instead, we copy the bounds until
+  // the starting node before it is detached from the main tree and use it
+  // as the starting bounds.
+  std::vector<f_t> start_lower;
+  std::vector<f_t> start_upper;
+  mip_node_t<i_t, f_t>* start_node;
+
+  pcgenerator_t rng;
+
+  bool recompute_basis                    = true;
+  bool recompute_bounds                   = true;
+  omp_atomic_t<bool> start_bounds_updated = false;
+
+  branch_and_bound_worker_t(i_t worker_id,
+                            const lp_problem_t<i_t, f_t>& original_lp,
+                            const csr_matrix_t<i_t, f_t>& Arow,
+                            const std::vector<variable_type_t>& var_type,
+                            const simplex_solver_settings_t<i_t, f_t>& settings)
+    : worker_id(worker_id),
+      search_strategy(BEST_FIRST),
+      is_active(false),
+      lower_bound(-std::numeric_limits<f_t>::infinity()),
+      leaf_problem(original_lp),
+      leaf_solution(original_lp.num_rows, original_lp.num_cols),
+      basis_factors(original_lp.num_rows, settings.refactor_frequency),
+      basic_list(original_lp.num_rows),
+      nonbasic_list(),
+      node_presolver(leaf_problem, Arow, {}, var_type),
+      bounds_changed(original_lp.num_cols, false),
+      start_node(nullptr),
+      rng(settings.random_seed + pcgenerator_t::default_seed + worker_id,
+          pcgenerator_t::default_stream ^ worker_id)
+  {
+  }
+
+  // Initialize the worker for plunging, setting the `start_node`, `start_lower` and
+  // `start_upper`. Returns `true` if no bounds were violated in any of the previous nodes
+  bool init_best_first(mip_node_t<i_t, f_t>* node, lp_problem_t<i_t, f_t>& original_lp)
+  {
+    bool feasible = node->check_variable_bounds(original_lp.lower, original_lp.upper);
+    if (!feasible) { return false; }
+
+    start_node      = node;
+    start_lower     = original_lp.lower;
+    start_upper     = original_lp.upper;
+    search_strategy = BEST_FIRST;
+    lower_bound     = node->lower_bound;
+    is_active       = true;
+    return true;
+  }
+
+  // Initialize the worker for diving, setting the `start_node`, `start_lower` and
+  // `start_upper`. Returns `true` if the starting node is feasible via
+  // bounds propagation and no bounds were violated in any of the previous nodes
+  bool init_diving(mip_node_t<i_t, f_t>* node_ptr,
+                   search_strategy_t type,
+                   const lp_problem_t<i_t, f_t>& original_lp,
+                   const simplex_solver_settings_t<i_t, f_t>& settings)
+  {
+    internal_node   = node_ptr->detach_copy();
+    start_node      = &internal_node;
+    start_lower     = original_lp.lower;
+    start_upper     = original_lp.upper;
+    search_strategy = type;
+    lower_bound     = node_ptr->lower_bound;
+    std::fill(bounds_changed.begin(), bounds_changed.end(), false);
+
+    bool feasible = node_ptr->get_variable_bounds(
+      original_lp.lower, original_lp.upper, start_lower, start_upper, bounds_changed);
+    if (feasible) {
+      feasible =
+        node_presolver.bounds_strengthening(settings, bounds_changed, start_lower, start_upper);
+    }
+    is_active = feasible;
+    return feasible;
+  }
+
+  // Set the variables bounds for the LP relaxation in the current node.
+  bool set_lp_variable_bounds(mip_node_t<i_t, f_t>* node_ptr,
+                              const simplex_solver_settings_t<i_t, f_t>& settings)
+  {
+    bool feasible = false;
+
+    // Reset the bound_changed markers
+    std::fill(bounds_changed.begin(), bounds_changed.end(), false);
+
+    // Set the correct bounds for the leaf problem
+    if (recompute_bounds) {
+      feasible = node_ptr->get_variable_bounds(
+        start_lower, start_upper, leaf_problem.lower, leaf_problem.upper, bounds_changed);
+    } else {
+      feasible = node_ptr->update_branched_variable_bounds(
+        start_lower, start_upper, leaf_problem.lower, leaf_problem.upper, bounds_changed);
+    }
+
+    if (feasible) {
+      feasible = node_presolver.bounds_strengthening(
+        settings, bounds_changed, leaf_problem.lower, leaf_problem.upper);
+    }
+    return feasible;
+  }
+
+ private:
+  // For diving, we need to store the full node instead
+  // of just a pointer, since it is not stored in the tree anymore.
+  // To keep the same interface across all worker types,
+  // this will be used as a temporary storage and
+  // will be pointed by `start_node`.
+  // For exploration, this will not be used.
+  mip_node_t<i_t, f_t> internal_node;
+};
+
+}  // namespace cuopt::linear_programming::dual_simplex
diff --git a/cpp/src/branch_and_bound/worker_pool.hpp b/cpp/src/branch_and_bound/worker_pool.hpp
new file mode 100644
index 0000000000..7ea1a0c184
--- /dev/null
+++ b/cpp/src/branch_and_bound/worker_pool.hpp
@@ -0,0 +1,139 @@
+/* clang-format off */
+/*
+ * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-License-Identifier: Apache-2.0
+ */
+/* clang-format on */
+
+#pragma once
+
+#include <branch_and_bound/worker.hpp>
+
+namespace cuopt::linear_programming::dual_simplex {
+
+template <typename i_t, typename f_t>
+class branch_and_bound_worker_pool_t {
+ public:
+  void init(i_t num_workers,
+            const lp_problem_t<i_t, f_t>& original_lp,
+            const csr_matrix_t<i_t, f_t>& Arow,
+            const std::vector<variable_type_t>& var_type,
+            const simplex_solver_settings_t<i_t, f_t>& settings)
+  {
+    workers_.resize(num_workers);
+    num_idle_workers_ = num_workers;
+    for (i_t i = 0; i < num_workers; ++i) {
+      workers_[i] = std::make_unique<branch_and_bound_worker_t<i_t, f_t>>(
+        i, original_lp, Arow, var_type, settings);
+      idle_workers_.push_front(i);
+    }
+
+    is_initialized = true;
+  }
+
+  // Here, we are assuming that the scheduler is the only
+  // thread that can retrieve/pop an idle worker.
+  branch_and_bound_worker_t<i_t, f_t>* get_idle_worker()
+  {
+    std::lock_guard<omp_mutex_t> lock(mutex_);
+    if (idle_workers_.empty()) {
+      return nullptr;
+    } else {
+      i_t idx = idle_workers_.front();
+      return workers_[idx].get();
+    }
+  }
+
+  // Here, we are assuming that the scheduler is the only
+  // thread that can retrieve/pop an idle worker.
+  void pop_idle_worker()
+  {
+    std::lock_guard<omp_mutex_t> lock(mutex_);
+    if (!idle_workers_.empty()) {
+      idle_workers_.pop_front();
+      num_idle_workers_--;
+    }
+  }
+
+  void return_worker_to_pool(branch_and_bound_worker_t<i_t, f_t>* worker)
+  {
+    worker->is_active = false;
+    std::lock_guard<omp_mutex_t> lock(mutex_);
+    idle_workers_.push_back(worker->worker_id);
+    num_idle_workers_++;
+  }
+
+  f_t get_lower_bound()
+  {
+    f_t lower_bound = std::numeric_limits<f_t>::infinity();
+
+    if (is_initialized) {
+      for (i_t i = 0; i < workers_.size(); ++i) {
+        if (workers_[i]->search_strategy == BEST_FIRST && workers_[i]->is_active) {
+          lower_bound = std::min(workers_[i]->lower_bound.load(), lower_bound);
+        }
+      }
+    }
+
+    return lower_bound;
+  }
+
+  i_t num_idle_workers() { return num_idle_workers_; }
+
+  void broadcast_root_bounds_change()
+  {
+    if (is_initialized) {
+      for (i_t i = 0; i < workers_.size(); ++i) {
+        if (workers_[i]->is_active) { workers_[i]->start_bounds_updated = true; }
+      }
+    }
+  }
+
+ private:
+  // Worker pool
+  std::vector<std::unique_ptr<branch_and_bound_worker_t<i_t, f_t>>> workers_;
+  bool is_initialized = false;
+
+  omp_mutex_t mutex_;
+  std::deque<i_t> idle_workers_;
+  omp_atomic_t<i_t> num_idle_workers_;
+};
+
+template <typename i_t, typename f_t>
+std::vector<search_strategy_t> get_search_strategies(
+  diving_heuristics_settings_t<i_t, f_t> settings)
+{
+  std::vector<search_strategy_t> types;
+  types.reserve(num_search_strategies);
+  types.push_back(BEST_FIRST);
+  if (settings.pseudocost_diving != 0) { types.push_back(PSEUDOCOST_DIVING); }
+  if (settings.line_search_diving != 0) { types.push_back(LINE_SEARCH_DIVING); }
+  if (settings.guided_diving != 0) { types.push_back(GUIDED_DIVING); }
+  if (settings.coefficient_diving != 0) { types.push_back(COEFFICIENT_DIVING); }
+  return types;
+}
+
+template <typename i_t>
+std::array<i_t, num_search_strategies> get_max_workers(
+  i_t num_workers, const std::vector<search_strategy_t>& strategies)
+{
+  std::array<i_t, num_search_strategies> max_num_workers;
+  max_num_workers.fill(0);
+
+  i_t bfs_workers             = std::max(strategies.size() == 1 ? num_workers : num_workers / 4, 1);
+  max_num_workers[BEST_FIRST] = bfs_workers;
+
+  i_t diving_workers = (num_workers - bfs_workers);
+  i_t m              = strategies.size() - 1;
+
+  for (size_t i = 1, k = 0; i < strategies.size(); ++i) {
+    i_t start                      = (double)k * diving_workers / m;
+    i_t end                        = (double)(k + 1) * diving_workers / m;
+    max_num_workers[strategies[i]] = end - start;
+    ++k;
+  }
+
+  return max_num_workers;
+}
+
+}  // namespace cuopt::linear_programming::dual_simplex
diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp
index eadd93040c..0372771a3d 100644
--- a/cpp/src/dual_simplex/simplex_solver_settings.hpp
+++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp
@@ -183,10 +183,16 @@ struct simplex_solver_settings_t {
   i_t clique_cuts;                 // -1 automatic, 0 to disable, >0 to enable clique cuts
   i_t strong_chvatal_gomory_cuts;  // -1 automatic, 0 to disable, >0 to enable strong Chvatal Gomory
                                    // cuts
-  i_t reduced_cost_strengthening;  // -1 automatic, 0 to disable, >0 to enable reduced cost
-                                   // strengthening
-  f_t cut_change_threshold;        // threshold for cut change
-  f_t cut_min_orthogonality;       // minimum orthogonality for cuts
+  // >= 0 adds additional reduced cost strengthening passes at each level
+  // -1 - automatic
+  // 0 - disable
+  // 1 - apply to root after each cut pass
+  // 2 - apply to root after all cuts
+  // 3 - apply to root after each incumbent update
+  i_t reduced_cost_strengthening;
+
+  f_t cut_change_threshold;                // threshold for cut change
+  f_t cut_min_orthogonality;               // minimum orthogonality for cuts
   i_t mip_batch_pdlp_strong_branching{0};  // 0 if not using batch PDLP for strong branching, 1 if
                                            // using batch PDLP for strong branching
 
diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu
index d7d7601014..de4b2557fb 100644
--- a/cpp/src/mip_heuristics/diversity/lns/rins.cu
+++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu
@@ -264,12 +264,13 @@ void rins_t<i_t, f_t>::run_rins()
     std::min(current_mip_gap, (f_t)settings.target_mip_gap);
   branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance;
   branch_and_bound_settings.num_threads = 1;
-  branch_and_bound_settings.reliability_branching = 0;
-  branch_and_bound_settings.max_cut_passes        = 0;
-  branch_and_bound_settings.clique_cuts           = 0;
-  branch_and_bound_settings.sub_mip               = 1;
-  branch_and_bound_settings.log.log               = false;
-  branch_and_bound_settings.log.log_prefix        = "[RINS] ";
+  branch_and_bound_settings.reliability_branching      = 0;
+  branch_and_bound_settings.reduced_cost_strengthening = 0;
+  branch_and_bound_settings.max_cut_passes             = 0;
+  branch_and_bound_settings.clique_cuts                = 0;
+  branch_and_bound_settings.sub_mip                    = 1;
+  branch_and_bound_settings.log.log                    = false;
+  branch_and_bound_settings.log.log_prefix             = "[RINS] ";
   branch_and_bound_settings.solution_callback = [&rins_solution_queue](std::vector<f_t>& solution,
                                                                        f_t objective) {
     rins_solution_queue.push_back(solution);
diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh
index a867141d0a..43aae14204 100644
--- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh
+++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh
@@ -106,11 +106,12 @@ class sub_mip_recombiner_t : public recombiner_t<i_t, f_t> {
       branch_and_bound_settings.relative_mip_gap_tol = context.settings.tolerances.relative_mip_gap;
       branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance;
       branch_and_bound_settings.num_threads = 1;
-      branch_and_bound_settings.reliability_branching = 0;
-      branch_and_bound_settings.max_cut_passes        = 0;
-      branch_and_bound_settings.clique_cuts           = 0;
-      branch_and_bound_settings.sub_mip               = 1;
-      branch_and_bound_settings.solution_callback     = [this](std::vector<f_t>& solution,
+      branch_and_bound_settings.reliability_branching      = 0;
+      branch_and_bound_settings.reduced_cost_strengthening = 0;
+      branch_and_bound_settings.max_cut_passes             = 0;
+      branch_and_bound_settings.clique_cuts                = 0;
+      branch_and_bound_settings.sub_mip                    = 1;
+      branch_and_bound_settings.solution_callback          = [this](std::vector<f_t>& solution,
                                                            f_t objective) {
         this->solution_callback(solution, objective);
       };
diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu
index f25c093afb..b428940fb3 100644
--- a/cpp/src/mip_heuristics/solver.cu
+++ b/cpp/src/mip_heuristics/solver.cu
@@ -224,15 +224,12 @@ solution_t<i_t, f_t> mip_solver_t<i_t, f_t>::run_solver()
     branch_and_bound_settings.strong_chvatal_gomory_cuts =
       context.settings.strong_chvatal_gomory_cuts;
     branch_and_bound_settings.reduced_cost_strengthening =
-      context.settings.reduced_cost_strengthening;
+      context.settings.reduced_cost_strengthening >= 0 ? context.settings.reduced_cost_strengthening
+                                                       : 3;
     branch_and_bound_settings.cut_change_threshold  = context.settings.cut_change_threshold;
     branch_and_bound_settings.cut_min_orthogonality = context.settings.cut_min_orthogonality;
     branch_and_bound_settings.mip_batch_pdlp_strong_branching =
       context.settings.mip_batch_pdlp_strong_branching;
-    branch_and_bound_settings.reduced_cost_strengthening =
-      context.settings.reduced_cost_strengthening == -1
-        ? 2
-        : context.settings.reduced_cost_strengthening;
 
     if (context.settings.num_cpu_threads < 0) {
       branch_and_bound_settings.num_threads = std::max(1, omp_get_max_threads() - 1);