diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp
index e786f35a59..c5a31cb4ad 100644
--- a/cpp/src/branch_and_bound/branch_and_bound.cpp
+++ b/cpp/src/branch_and_bound/branch_and_bound.cpp
@@ -222,6 +222,7 @@ branch_and_bound_t<i_t, f_t>::branch_and_bound_t(
   const user_problem_t<i_t, f_t>& user_problem,
   const simplex_solver_settings_t<i_t, f_t>& solver_settings,
   f_t start_time,
+  std::atomic<int>* restart_concurrent_halt,
   const probing_implied_bound_t<i_t, f_t>& probing_implied_bound,
   std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table,
   mip_symmetry_t<i_t, f_t>* symmetry)
@@ -235,6 +236,7 @@ branch_and_bound_t<i_t, f_t>::branch_and_bound_t(
     incumbent_(1),
     root_relax_soln_(1, 1),
     root_crossover_soln_(1, 1),
+    restart_concurrent_halt_(restart_concurrent_halt),
     pc_(1, solver_settings),
     solver_status_(mip_status_t::UNSET)
 {
@@ -271,6 +273,7 @@ branch_and_bound_t<i_t, f_t>::branch_and_bound_t(
   upper_bound_                 = inf;
   root_objective_              = std::numeric_limits<f_t>::quiet_NaN();
   root_lp_current_lower_bound_ = -inf;
+  pc_.concurrent_halt_         = &node_concurrent_halt_;
 }
 
 template <typename i_t, typename f_t>
@@ -714,7 +717,7 @@ void branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t, f_t>&
   bool is_maximization = original_lp_.obj_scale < 0.0;
 
   settings_.log.printf("Explored %d nodes in %.2fs.\n",
-                       exploration_stats_.nodes_explored,
+                       exploration_stats_.total_nodes_explored,
                        toc(exploration_stats_.start_time));
   if (exploration_stats_.orbital_fixing_nodes.load() > 0 ||
       exploration_stats_.orbital_conflict_nodes.load() > 0) {
@@ -775,7 +778,7 @@ void branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t, f_t>&
     solution.objective = incumbent_.objective;
   }
   solution.lower_bound        = lower_bound;
-  solution.nodes_explored     = exploration_stats_.nodes_explored;
+  solution.nodes_explored     = exploration_stats_.total_nodes_explored;
   solution.simplex_iterations = exploration_stats_.total_lp_iters;
 }
 
@@ -1517,10 +1520,60 @@ dual::status_t branch_and_bound_t<i_t, f_t>::solve_node_lp(
   return lp_status;
 }
 
+template <typename i_t, typename f_t>
+bool branch_and_bound_t<i_t, f_t>::should_restart(f_t current_abs_gap)
+{
+  if (settings_.sub_mip || restart_count_ >= settings_.max_restarts) return false;
+
+  i_t num_nodes   = exploration_stats_.nodes_explored;
+  i_t total_nodes = exploration_stats_.total_nodes_explored;
+  if (num_nodes < settings_.restart_min_nodes) return false;
+
+  i_t nodes_since_last_check = num_nodes - exploration_stats_.restart_nodes_at_last_check;
+  if (nodes_since_last_check < settings_.restart_check_freq) return false;
+
+  f_t current_progress = search_tree_.progress;
+  f_t progress_since_last_check =
+    std::max(current_progress - exploration_stats_.restart_progress_at_last_check, 1E-6);
+  i_t tree_size_estimate =
+    exploration_stats_.restart_nodes_at_last_check +
+    nodes_since_last_check * (1.0 - current_progress) / progress_since_last_check;
+
+  f_t gap_reduction = exploration_stats_.restart_gap_at_last_check / current_abs_gap;
+
+  settings_.log.debug(
+    "[Restart] Current: explored=%d, progress=%.4f, gap=%.4f. Since last: explored=%d, "
+    "progress=%.4f, gap=%.4f. Tree size estimate=%d",
+    num_nodes,
+    current_progress,
+    current_abs_gap,
+    nodes_since_last_check,
+    progress_since_last_check,
+    gap_reduction,
+    tree_size_estimate);
+
+  if (gap_reduction < 1.05 &&
+      tree_size_estimate >= settings_.restart_tree_size_factor * total_nodes) {
+    ++exploration_stats_.restart_large_tree_count;
+    i_t min_count =
+      settings_.restart_min_estimates + total_nodes * settings_.restart_threshold_grow_per_node;
+    return exploration_stats_.restart_large_tree_count >=
+           min_count * std::pow(settings_.restart_threshold_grow_per_restart, restart_count_);
+  }
+
+  exploration_stats_.restart_large_tree_count       = 0;
+  exploration_stats_.restart_gap_at_last_check      = current_abs_gap;
+  exploration_stats_.restart_progress_at_last_check = current_progress;
+  exploration_stats_.restart_nodes_at_last_check    = num_nodes;
+  return false;
+}
+
 template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::plunge_with(bfs_worker_t<i_t, f_t>* worker,
                                                mip_node_t<i_t, f_t>* start_node)
 {
+  raft::common::nvtx::range scope_guess("BB::plunge");
+
   assert(worker != nullptr && worker->is_active.load());
   assert(start_node != nullptr);
 
@@ -1540,16 +1593,20 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(bfs_worker_t<i_t, f_t>* worker,
 
   while (stack.size() > 0 && (solver_status_ == mip_status_t::UNSET && is_running_) &&
          rel_gap > settings_.relative_mip_gap_tol && abs_gap > settings_.absolute_mip_gap_tol) {
-    if (worker->worker_id == 0) { repair_heuristic_solutions(); }
+    if (worker->worker_id == 0) {
+      if (should_restart(abs_gap)) {
+        mip_node_t<i_t, f_t>* node = stack.front();
+        report(' ', upper_bound_, lower_bound, node->depth, node->integer_infeasible);
+        solver_status_            = mip_status_t::RESTART;
+        node_concurrent_halt_     = 1;
+        *restart_concurrent_halt_ = 1;
+        break;
+      }
 
-    if (worker->total_active_diving_workers < worker->total_max_diving_workers &&
-        worker->node_queue.diving_queue_size() > 0) {
-      launch_diving_worker(worker);
+      repair_heuristic_solutions();
     }
 
-    if (bfs_worker_pool_.num_idle() > 0 && worker->node_queue.best_first_queue_size() > 0) {
-      launch_bfs_worker(worker);
-    }
+    if (*restart_concurrent_halt_ == 1) { break; }
 
     assert(stack.size() <= 2);
     mip_node_t<i_t, f_t>* node_ptr = stack.front();
@@ -1597,7 +1654,7 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(bfs_worker_t<i_t, f_t>* worker,
       break;
     }
 
-    if (exploration_stats_.nodes_explored + exploration_stats_.nodes_being_solved >
+    if (exploration_stats_.total_nodes_explored + exploration_stats_.nodes_being_solved >
         settings_.node_limit) {
       solver_status_ = mip_status_t::NODE_LIMIT;
       stack.push_front(node_ptr);
@@ -1608,6 +1665,7 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(bfs_worker_t<i_t, f_t>* worker,
     dual::status_t lp_status = solve_node_lp(node_ptr, worker, exploration_stats_, settings_.log);
     ++exploration_stats_.nodes_since_last_log;
     ++exploration_stats_.nodes_explored;
+    ++exploration_stats_.total_nodes_explored;
     --exploration_stats_.nodes_unexplored;
     --exploration_stats_.nodes_being_solved;
 
@@ -1669,6 +1727,15 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(bfs_worker_t<i_t, f_t>* worker,
     upper_bound = upper_bound_;
     rel_gap     = user_relative_gap(original_lp_, upper_bound, lower_bound);
     abs_gap     = compute_user_abs_gap(original_lp_, upper_bound, lower_bound);
+
+    if (worker->total_active_diving_workers < worker->total_max_diving_workers &&
+        worker->node_queue.diving_queue_size() > 0) {
+      launch_diving_worker(worker);
+    }
+
+    if (bfs_worker_pool_.num_idle() > 0 && worker->node_queue.best_first_queue_size() > 0) {
+      launch_bfs_worker(worker);
+    }
   }
 
   // If the solver exits early without consuming the local stack, or converged according to
@@ -1688,6 +1755,8 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(bfs_worker_t<i_t, f_t>* worker,
 template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::launch_bfs_worker(bfs_worker_t<i_t, f_t>* worker)
 {
+  raft::common::nvtx::range scope_guess("BB::launch_bfs_worker");
+
   bfs_worker_t<i_t, f_t>* idle_worker = bfs_worker_pool_.pop_idle_worker();
   if (!idle_worker) return;
 
@@ -1721,6 +1790,8 @@ void branch_and_bound_t<i_t, f_t>::launch_bfs_worker(bfs_worker_t<i_t, f_t>* wor
 template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::work_stealing(bfs_worker_t<i_t, f_t>* worker)
 {
+  raft::common::nvtx::range scope_guess("BB::work_stealing");
+
   i_t nodes_to_steal = settings_.bnb_nodes_per_steal >= 0 ? settings_.bnb_nodes_per_steal
                                                           : MIP_DEFAULT_NODES_PER_STEAL;
   i_t max_attempts   = settings_.bnb_max_steal_attempts >= 0 ? settings_.bnb_max_steal_attempts
@@ -1735,6 +1806,8 @@ void branch_and_bound_t<i_t, f_t>::work_stealing(bfs_worker_t<i_t, f_t>* worker)
 template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::best_first_search_with(bfs_worker_t<i_t, f_t>* worker)
 {
+  raft::common::nvtx::range scope_guess("BB::bfs_worker");
+
   f_t lower_bound = get_lower_bound();
   f_t abs_gap     = compute_user_abs_gap(original_lp_, upper_bound_.load(), lower_bound);
   f_t rel_gap     = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
@@ -1768,6 +1841,8 @@ void branch_and_bound_t<i_t, f_t>::best_first_search_with(bfs_worker_t<i_t, f_t>
       break;
     }
 
+    if (*restart_concurrent_halt_ == 1) { break; }
+
     // Pre-emptively set the lower bound of the worker
     worker->lower_bound              = node_queue.get_lower_bound();
     mip_node_t<i_t, f_t>* start_node = node_queue.pop();
@@ -1795,6 +1870,8 @@ void branch_and_bound_t<i_t, f_t>::best_first_search_with(bfs_worker_t<i_t, f_t>
       break;
     }
 
+    if (solver_status_ == mip_status_t::RESTART) { break; }
+
     // Steal a node with some probability or when it is empty. The victim is determined at random.
     if (node_queue.best_first_queue_size() == 0 || worker->rng.next_double() < steal_chance) {
       work_stealing(worker);
@@ -1812,7 +1889,7 @@ void branch_and_bound_t<i_t, f_t>::best_first_search_with(bfs_worker_t<i_t, f_t>
 template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::dive_with(diving_worker_t<i_t, f_t>* worker)
 {
-  raft::common::nvtx::range scope("BB::diving_thread");
+  raft::common::nvtx::range scope("BB::diving_worker");
   if (worker->orbital_fixing) { worker->orbital_fixing->disable(); }
   logger_t log;
   log.log = false;
@@ -1861,6 +1938,7 @@ void branch_and_bound_t<i_t, f_t>::dive_with(diving_worker_t<i_t, f_t>* worker)
       break;
     }
     if (dive_stats.nodes_explored >= diving_node_limit) { break; }
+    if (*restart_concurrent_halt_ == 1) { break; }
 
     dual::status_t lp_status = solve_node_lp(node_ptr, worker, dive_stats, log);
     ++dive_stats.nodes_explored;
@@ -1906,6 +1984,8 @@ void branch_and_bound_t<i_t, f_t>::dive_with(diving_worker_t<i_t, f_t>* worker)
 template <typename i_t, typename f_t>
 bool branch_and_bound_t<i_t, f_t>::launch_diving_worker(bfs_worker_t<i_t, f_t>* bfs_worker)
 {
+  raft::common::nvtx::range scope_guess("BB::launch_diving_worker");
+
   // Get an idle worker.
   diving_worker_t<i_t, f_t>* diving_worker = diving_worker_pool_.pop_idle_worker();
   if (diving_worker == nullptr) { return false; }
@@ -1975,6 +2055,8 @@ lp_status_t branch_and_bound_t<i_t, f_t>::solve_root_relaxation(
   std::vector<i_t>& nonbasic_list,
   std::vector<f_t>& edge_norms)
 {
+  raft::common::nvtx::range scope_guess("BB::solve_root_relaxation");
+
   f_t start_time          = tic();
   f_t user_objective      = 0;
   i_t iter                = 0;
@@ -2146,6 +2228,8 @@ auto branch_and_bound_t<i_t, f_t>::do_cut_pass(
   i_t& cut_pool_size,
   [[maybe_unused]] const std::vector<f_t>& saved_solution) -> cut_pass_result_t
 {
+  raft::common::nvtx::range scope_guess("BB::cut_pass");
+
 #ifdef PRINT_FRACTIONAL_INFO
   settings_.log.printf("Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass);
   for (i_t j : fractional) {
@@ -2413,13 +2497,16 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   raft::common::nvtx::range scope("BB::solve");
 
   logger_t log;
-  log.log                             = false;
-  log.log_prefix                      = settings_.log.log_prefix;
-  solver_status_                      = mip_status_t::UNSET;
-  is_running_                         = false;
-  root_lp_current_lower_bound_        = -inf;
-  exploration_stats_.nodes_unexplored = 0;
-  exploration_stats_.nodes_explored   = 0;
+  log.log                                 = false;
+  log.log_prefix                          = settings_.log.log_prefix;
+  solver_status_                          = mip_status_t::UNSET;
+  is_running_                             = false;
+  restart_count_                          = 0;
+  min_node_queue_size_                    = 20;
+  root_lp_current_lower_bound_            = -inf;
+  exploration_stats_.nodes_unexplored     = 0;
+  exploration_stats_.total_nodes_explored = 0;
+  exploration_stats_.nodes_explored       = 0;
   original_lp_.A.to_compressed_row(Arow_);
 
   settings_.log.printf("Reduced cost strengthening enabled: %d\n",
@@ -2807,128 +2894,151 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
     return solver_status_;
   }
 
-  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;
+  if (settings_.diving_settings.coefficient_diving != 0) {
+    calculate_variable_locks(original_lp_, var_up_locks_, var_down_locks_);
+  }
 
-      bounds_strengthening_t<i_t, f_t> node_presolve(original_lp_, Arow_, row_sense, var_types_);
+  const i_t num_workers        = settings_.num_threads;
+  const i_t num_bfs_workers    = std::max(settings_.num_threads / 2, 1);
+  const i_t num_diving_workers = num_workers - num_bfs_workers;
+  bfs_worker_pool_.init(num_bfs_workers, original_lp_, Arow_, var_types_, symmetry_, settings_);
 
-      mutex_original_lp_.lock();
-      original_lp_.lower = lower_bounds;
-      original_lp_.upper = upper_bounds;
-      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
-                                         // strengthening thinks we are infeasible.
-      }
-      // Go through and check the fractional variables and remove any that are now fixed to their
-      // bounds
-      std::vector<i_t> to_remove(fractional.size(), 0);
-      i_t num_to_remove = 0;
-      for (i_t k = 0; k < fractional.size(); k++) {
-        const i_t j = fractional[k];
-        if (std::abs(original_lp_.upper[j] - original_lp_.lower[j]) < settings_.fixed_tol) {
-          to_remove[k] = 1;
-          num_to_remove++;
+  if (num_diving_workers > 0) {
+    diving_worker_pool_.init(
+      num_diving_workers, original_lp_, Arow_, var_types_, symmetry_, settings_, num_bfs_workers);
+  }
+
+  do {
+    if (settings_.reduced_cost_strengthening >= 2 && upper_bound_.load() < last_upper_bound) {
+      last_upper_bound = upper_bound_.load();
+
+      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;
+
+        bounds_strengthening_t<i_t, f_t> node_presolve(original_lp_, Arow_, row_sense, var_types_);
+
+        mutex_original_lp_.lock();
+        original_lp_.lower = lower_bounds;
+        original_lp_.upper = upper_bounds;
+        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
+                                           // strengthening thinks we are infeasible.
         }
-      }
-      if (num_to_remove > 0) {
-        std::vector<i_t> new_fractional;
-        new_fractional.reserve(fractional.size() - num_to_remove);
+        // Go through and check the fractional variables and remove any that are now fixed to their
+        // bounds
+        std::vector<i_t> to_remove(fractional.size(), 0);
+        i_t num_to_remove = 0;
         for (i_t k = 0; k < fractional.size(); k++) {
-          if (!to_remove[k]) { new_fractional.push_back(fractional[k]); }
+          const i_t j = fractional[k];
+          if (std::abs(original_lp_.upper[j] - original_lp_.lower[j]) < settings_.fixed_tol) {
+            to_remove[k] = 1;
+            num_to_remove++;
+          }
+        }
+        if (num_to_remove > 0) {
+          std::vector<i_t> new_fractional;
+          new_fractional.reserve(fractional.size() - num_to_remove);
+          for (i_t k = 0; k < fractional.size(); k++) {
+            if (!to_remove[k]) { new_fractional.push_back(fractional[k]); }
+          }
+          fractional     = new_fractional;
+          num_fractional = fractional.size();
         }
-        fractional     = new_fractional;
-        num_fractional = fractional.size();
       }
     }
-  }
 
-  // Choose variable to branch on
-  i_t branch_var = pc_.variable_selection(fractional, root_relax_soln_.x);
-
-  search_tree_.root      = std::move(mip_node_t<i_t, f_t>(root_objective_, root_vstatus_));
-  search_tree_.num_nodes = 0;
-  search_tree_.graphviz_node(settings_.log, &search_tree_.root, "lower bound", root_objective_);
-  search_tree_.branch(&search_tree_.root,
-                      branch_var,
-                      root_relax_soln_.x[branch_var],
-                      num_fractional,
-                      root_vstatus_,
-                      original_lp_,
-                      log);
-
-  if (symmetry_ != nullptr) {
-    i_t removed =
-      symmetry_->generators.template prune_by_bounds<f_t>(original_lp_.lower, original_lp_.upper);
-    if (removed > 0) {
-      symmetry_->num_generators = static_cast<int>(symmetry_->generators.num_generators());
-      settings_.log.printf(
-        "Pruned %d generators invalidated by root-level bound tightening, %d remain\n",
-        removed,
-        symmetry_->num_generators);
+    if (symmetry_ != nullptr) {
+      i_t removed =
+        symmetry_->generators.template prune_by_bounds<f_t>(original_lp_.lower, original_lp_.upper);
+      if (removed > 0) {
+        symmetry_->num_generators = static_cast<int>(symmetry_->generators.num_generators());
+        settings_.log.printf(
+          "Pruned %d generators invalidated by root-level bound tightening, %d remain\n",
+          removed,
+          symmetry_->num_generators);
+      }
     }
-  }
 
-  settings_.log.printf("Exploring the B&B tree using %d threads\n\n", settings_.num_threads);
-  node_concurrent_halt_ = 0;
-
-  exploration_stats_.nodes_explored       = 0;
-  exploration_stats_.nodes_unexplored     = 2;
-  exploration_stats_.nodes_since_last_log = 0;
-  exploration_stats_.last_log             = tic();
-  min_node_queue_size_                    = 20;
+    if (toc(exploration_stats_.start_time) > settings_.time_limit) {
+      solver_status_ = mip_status_t::TIME_LIMIT;
+      break;
+    }
 
-  if (settings_.diving_settings.coefficient_diving != 0) {
-    calculate_variable_locks(original_lp_, var_up_locks_, var_down_locks_);
-  }
+    if (solver_status_ == mip_status_t::RESTART) {
+      settings_.log.printf("\nRestarting B&B after %.2fs and %d nodes\n",
+                           toc(exploration_stats_.start_time),
+                           exploration_stats_.nodes_explored.load());
+      search_tree_.clean();
+      bfs_worker_pool_.reset();
+      diving_worker_pool_.reset();
+      solver_status_ = mip_status_t::UNSET;
+      ++restart_count_;
+    }
 
-  if (settings_.deterministic) {
-    settings_.log.printf(
-      " | Explored | Unexplored |    Objective    |     Bound     | IntInf | Depth | Iter/Node "
-      "|   Gap    |  Work |  Time  |\n");
-  } else {
-    settings_.log.printf(
-      " | Explored | Unexplored |    Objective    |     Bound     | IntInf | Depth | Iter/Node "
-      "|   Gap    |  Time  |\n");
-  }
+    // Choose variable to branch on
+    i_t branch_var = pc_.variable_selection(fractional, root_relax_soln_.x);
+
+    search_tree_.root      = std::move(mip_node_t<i_t, f_t>(root_objective_, root_vstatus_));
+    search_tree_.num_nodes = 0;
+    search_tree_.graphviz_node(settings_.log, &search_tree_.root, "lower bound", root_objective_);
+    search_tree_.branch(&search_tree_.root,
+                        branch_var,
+                        root_relax_soln_.x[branch_var],
+                        num_fractional,
+                        root_vstatus_,
+                        original_lp_,
+                        log);
+
+    if (!settings_.sub_mip) *restart_concurrent_halt_ = 0;
+    node_concurrent_halt_ = 0;
+
+    exploration_stats_.nodes_explored           = 0;
+    exploration_stats_.nodes_unexplored         = 2;
+    exploration_stats_.nodes_since_last_log     = 0;
+    exploration_stats_.last_log                 = 0;
+    exploration_stats_.restart_large_tree_count = 0;
+    exploration_stats_.restart_gap_at_last_check =
+      compute_user_abs_gap(original_lp_, upper_bound_.load(), root_relax_objective);
+    exploration_stats_.restart_nodes_at_last_check    = 0;
+    exploration_stats_.restart_progress_at_last_check = 0;
 
-#pragma omp taskgroup
-  {
     if (settings_.deterministic) {
-      run_deterministic_coordinator(Arow_);
+      settings_.log.printf(
+        " | Explored | Unexplored |    Objective    |     Bound     | IntInf | Depth | Iter/Node "
+        "|   Gap    |  Work |  Time  |\n");
     } else {
-      const i_t num_workers        = settings_.num_threads;
-      const i_t num_bfs_workers    = std::max(settings_.num_threads / 2, 1);
-      const i_t num_diving_workers = num_workers - num_bfs_workers;
-      bfs_worker_pool_.init(num_bfs_workers, original_lp_, Arow_, var_types_, symmetry_, settings_);
-
-      if (num_diving_workers > 0) {
-        diving_worker_pool_.init(num_diving_workers,
-                                 original_lp_,
-                                 Arow_,
-                                 var_types_,
-                                 symmetry_,
-                                 settings_,
-                                 num_bfs_workers);
-      }
+      settings_.log.printf(
+        " | Explored | Unexplored |    Objective    |     Bound     | IntInf | Depth | Iter/Node "
+        "|   Gap    |  Time  |\n");
+    }
 
-      bfs_worker_t<i_t, f_t>* initial_worker = bfs_worker_pool_.pop_idle_worker();
-      node_queue_t<i_t, f_t>& node_queue     = initial_worker->node_queue;
-      node_queue.push_lockfree(search_tree_.root.get_down_child());
-      node_queue.push_lockfree(search_tree_.root.get_up_child());
-      initial_worker->lower_bound = initial_worker->node_queue.get_lower_bound();
-      initial_worker->set_active();
-      best_first_search_with(initial_worker);
+    if (restart_count_ == 0) {
+      settings_.log.printf("Exploring the B&B tree using %d threads\n\n", settings_.num_threads);
     }
-  }  // Implicit barrier for all tasks created within the group (RINS, B&B workers)
+
+#pragma omp taskgroup
+    {
+      if (settings_.deterministic) {
+        run_deterministic_coordinator(Arow_);
+      } else {
+        bfs_worker_t<i_t, f_t>* initial_worker = bfs_worker_pool_.pop_idle_worker();
+        node_queue_t<i_t, f_t>& node_queue     = initial_worker->node_queue;
+        node_queue.push_lockfree(search_tree_.root.get_down_child());
+        node_queue.push_lockfree(search_tree_.root.get_up_child());
+        initial_worker->lower_bound = initial_worker->node_queue.get_lower_bound();
+        initial_worker->set_active();
+        best_first_search_with(initial_worker);
+      }
+    }  // Implicit barrier for all tasks created within the group (RINS, B&B workers)
+  } while (solver_status_ == mip_status_t::RESTART);
 
   is_running_ = false;
 
@@ -3534,6 +3644,7 @@ node_status_t branch_and_bound_t<i_t, f_t>::solve_node_deterministic(
   exploration_stats_.total_lp_solve_time += toc(lp_start_time);
   exploration_stats_.total_lp_iters += node_iter;
   ++exploration_stats_.nodes_explored;
+  ++exploration_stats_.total_nodes_explored;
   --exploration_stats_.nodes_unexplored;
 
   deterministic_bfs_policy_t<i_t, f_t> policy{*this, worker};
diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp
index 0b32b2ece5..e891f676d9 100644
--- a/cpp/src/branch_and_bound/branch_and_bound.hpp
+++ b/cpp/src/branch_and_bound/branch_and_bound.hpp
@@ -12,6 +12,7 @@
 #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/search_tree.hpp>
 #include <branch_and_bound/worker.hpp>
 #include <branch_and_bound/worker_pool.hpp>
 
@@ -59,6 +60,7 @@ enum class mip_status_t {
   NUMERICAL  = 5,  // The solver encountered a numerical error
   UNSET      = 6,  // The status is not set
   WORK_LIMIT = 7,  // The solver reached a deterministic work limit
+  RESTART    = 8,  // The solver triggered a restart
 };
 
 template <typename i_t, typename f_t>
@@ -82,6 +84,7 @@ class branch_and_bound_t {
   branch_and_bound_t(const user_problem_t<i_t, f_t>& user_problem,
                      const simplex_solver_settings_t<i_t, f_t>& solver_settings,
                      f_t start_time,
+                     std::atomic<int>* restart_concurrent_halt,
                      const probing_implied_bound_t<i_t, f_t>& probing_implied_bound,
                      std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table = nullptr,
                      mip_symmetry_t<i_t, f_t>* symmetry                             = nullptr);
@@ -234,6 +237,7 @@ class branch_and_bound_t {
   bool enable_concurrent_lp_root_solve_{false};
   std::atomic<int> root_concurrent_halt_{0};
   std::atomic<int> node_concurrent_halt_{0};
+  std::atomic<int>* restart_concurrent_halt_{nullptr};
   bool is_root_solution_set{false};
 
   // Pseudocosts
@@ -263,6 +267,8 @@ class branch_and_bound_t {
   omp_atomic_t<f_t> lower_bound_numerical_;
   std::function<void(f_t)> user_bound_callback_;
 
+  i_t restart_count_;
+
   void report_heuristic(f_t obj);
   void report(char symbol,
               f_t obj,
@@ -314,6 +320,8 @@ class branch_and_bound_t {
   // Repairs low-quality solutions from the heuristics, if it is applicable.
   void repair_heuristic_solutions();
 
+  bool should_restart(f_t current_abs_gap);
+
   // Launch a new diving worker from a given best-first worker.
   bool launch_diving_worker(bfs_worker_t<i_t, f_t>* bfs_worker);
 
diff --git a/cpp/src/branch_and_bound/mip_node.hpp b/cpp/src/branch_and_bound/mip_node.hpp
index 2c0968ffce..ad53aaabb9 100644
--- a/cpp/src/branch_and_bound/mip_node.hpp
+++ b/cpp/src/branch_and_bound/mip_node.hpp
@@ -41,43 +41,8 @@ inline bool inactive_status(node_status_t status)
 template <typename i_t, typename f_t>
 class mip_node_t {
  public:
-  ~mip_node_t()
-  {
-    // Iterative teardown to avoid stack overflow on deep trees.
-    // Detach all descendants breadth-first, then destroy them as leaves.
-    // vector::push_back can throw bad_alloc; the catch-all keeps the destructor
-    // exception-free. Under OOM, any not-yet-detached descendants are destroyed
-    // via the recursive unique_ptr chain in `children` as this frame unwinds.
-    try {
-      std::vector<std::unique_ptr<mip_node_t>> nodes;
-      for (auto& c : children) {
-        if (c) { nodes.push_back(std::move(c)); }
-      }
-      // nodes.size() grows so that this loop only terminates when only leaves remain
-      for (size_t i = 0; i < nodes.size(); ++i) {
-        for (auto& c : nodes[i]->children) {
-          if (c) { nodes.push_back(std::move(c)); }
-        }
-      }
-
-      // scope-exit ensure destruction of all detached leaves
-    } catch (const std::exception& e) {
-      // fprintf to stderr is allocation-free and cannot throw; using the
-      // project logger here would risk a secondary bad_alloc that would
-      // escape the destructor and re-introduce std::terminate.
-      std::fprintf(stderr,
-                   "mip_node_t destructor: iterative teardown failed (%s); falling back to "
-                   "recursive unique_ptr destruction.\n",
-                   e.what());
-    } catch (...) {
-      std::fprintf(stderr,
-                   "mip_node_t destructor: iterative teardown failed (unknown exception); "
-                   "falling back to recursive unique_ptr destruction.\n");
-    }
-  }
-
-  mip_node_t(mip_node_t&&) noexcept            = default;
-  mip_node_t& operator=(mip_node_t&&) noexcept = default;
+  mip_node_t(mip_node_t&&)            = default;
+  mip_node_t& operator=(mip_node_t&&) = default;
 
   mip_node_t()
     : status(node_status_t::PENDING),
@@ -364,99 +329,4 @@ void remove_fathomed_nodes(std::vector<mip_node_t<i_t, f_t>*>& stack)
   }
 }
 
-template <typename i_t, typename f_t>
-class search_tree_t {
- public:
-  search_tree_t() : num_nodes(0) {}
-
-  search_tree_t(mip_node_t<i_t, f_t>&& node) : root(std::move(node)), num_nodes(0) {}
-
-  void update(mip_node_t<i_t, f_t>* node_ptr, node_status_t status)
-  {
-    std::lock_guard<omp_mutex_t> lock(mutex);
-    std::vector<mip_node_t<i_t, f_t>*> stack;
-    node_ptr->set_status(status, stack);
-    remove_fathomed_nodes(stack);
-  }
-
-  void branch(mip_node_t<i_t, f_t>* parent_node,
-              const i_t branch_var,
-              const f_t fractional_val,
-              const i_t integer_infeasible,
-              const std::vector<variable_status_t>& parent_vstatus,
-              const lp_problem_t<i_t, f_t>& original_lp,
-              logger_t& log)
-  {
-    i_t id = num_nodes.fetch_add(2);
-
-    auto down_child = std::make_unique<mip_node_t<i_t, f_t>>(original_lp,
-                                                             parent_node,
-                                                             ++id,
-                                                             branch_var,
-                                                             branch_direction_t::DOWN,
-                                                             fractional_val,
-                                                             integer_infeasible,
-                                                             parent_vstatus);
-    graphviz_edge(log,
-                  parent_node,
-                  down_child.get(),
-                  branch_var,
-                  branch_direction_t::DOWN,
-                  std::floor(fractional_val));
-
-    auto up_child = std::make_unique<mip_node_t<i_t, f_t>>(original_lp,
-                                                           parent_node,
-                                                           ++id,
-                                                           branch_var,
-                                                           branch_direction_t::UP,
-                                                           fractional_val,
-                                                           integer_infeasible,
-                                                           parent_vstatus);
-
-    graphviz_edge(log,
-                  parent_node,
-                  up_child.get(),
-                  branch_var,
-                  branch_direction_t::UP,
-                  std::ceil(fractional_val));
-
-    assert(parent_vstatus.size() == original_lp.num_cols);
-    parent_node->add_children(std::move(down_child),
-                              std::move(up_child));  // child pointers moved into the tree
-  }
-
-  void graphviz_node(logger_t& log,
-                     const mip_node_t<i_t, f_t>* node_ptr,
-                     const std::string label,
-                     const f_t val)
-  {
-    if (write_graphviz) {
-      log.printf("Node%d [label=\"%s %.16e\"]\n", node_ptr->node_id, label.c_str(), val);
-    }
-  }
-
-  void graphviz_edge(logger_t& log,
-                     const mip_node_t<i_t, f_t>* origin_ptr,
-                     const mip_node_t<i_t, f_t>* dest_ptr,
-                     const i_t branch_var,
-                     branch_direction_t branch_dir,
-                     const f_t bound)
-  {
-    if (write_graphviz) {
-      log.printf("Node%d -> Node%d [label=\"x%d %s %e\"]\n",
-                 origin_ptr->node_id,
-                 dest_ptr->node_id,
-                 branch_var,
-                 branch_dir == branch_direction_t::DOWN ? "<=" : ">=",
-                 bound);
-    }
-  }
-
-  mip_node_t<i_t, f_t> root;
-  omp_mutex_t mutex;
-  omp_atomic_t<i_t> num_nodes;
-
-  static constexpr bool write_graphviz = false;
-};
-
 }  // namespace cuopt::linear_programming::dual_simplex
diff --git a/cpp/src/branch_and_bound/node_queue.hpp b/cpp/src/branch_and_bound/node_queue.hpp
index b07b371b81..f63aa4587f 100644
--- a/cpp/src/branch_and_bound/node_queue.hpp
+++ b/cpp/src/branch_and_bound/node_queue.hpp
@@ -62,7 +62,7 @@ class heap_t {
 
   void clear()
   {
-    buffer.clear();
+    buffer       = {};
     num_entries_ = 0;
   }
 
@@ -157,6 +157,13 @@ class node_queue_t {
     return best_first_heap_.empty() ? std::numeric_limits<f_t>::infinity() : lower_bound_.load();
   }
 
+  void clear()
+  {
+    std::lock_guard lock(mutex_);
+    best_first_heap_.clear();
+    diving_heap_.clear();
+  }
+
  private:
   struct heap_entry_t {
     mip_node_t<i_t, f_t>* node = nullptr;
diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp
index fc3ee3bb3f..0ecf719bba 100644
--- a/cpp/src/branch_and_bound/pseudo_costs.cpp
+++ b/cpp/src/branch_and_bound/pseudo_costs.cpp
@@ -462,6 +462,7 @@ std::pair<f_t, dual::status_t> trial_branching(const lp_problem_t<i_t, f_t>& ori
                                                const basis_update_mpf_t<i_t, f_t>& basis_factors,
                                                const std::vector<i_t>& basic_list,
                                                const std::vector<i_t>& nonbasic_list,
+                                               std::atomic<int>* concurrent_halt,
                                                i_t branch_var,
                                                f_t branch_var_lower,
                                                f_t branch_var_upper,
@@ -482,6 +483,7 @@ std::pair<f_t, dual::status_t> trial_branching(const lp_problem_t<i_t, f_t>& ori
   child_settings.scale_columns   = false;
   child_settings.cut_off =
     objective_upper_bound(child_problem, upper_bound, child_settings.dual_tol);
+  child_settings.concurrent_halt = concurrent_halt;
 
   lp_solution_t<i_t, f_t> solution(original_lp.num_rows, original_lp.num_cols);
   iter                                             = 0;
@@ -1707,7 +1709,7 @@ i_t pseudo_costs_t<i_t, f_t>::reliable_variable_selection(
   std::vector<f_t> pdlp_obj_down(num_candidates, std::numeric_limits<f_t>::quiet_NaN());
   std::vector<f_t> pdlp_obj_up(num_candidates, std::numeric_limits<f_t>::quiet_NaN());
 
-  std::atomic<int> concurrent_halt{0};
+  std::atomic<int> pdlp_concurrent_halt{0};
 
   if (use_pdlp) {
 #pragma omp task default(shared) priority(CUOPT_HIGH_TASK_PRIORITY)
@@ -1715,7 +1717,7 @@ i_t pseudo_costs_t<i_t, f_t>::reliable_variable_selection(
                                           rb_mode,
                                           num_candidates,
                                           start_time,
-                                          concurrent_halt,
+                                          pdlp_concurrent_halt,
                                           original_lp,
                                           new_slacks,
                                           leaf_solution.x,
@@ -1731,7 +1733,7 @@ i_t pseudo_costs_t<i_t, f_t>::reliable_variable_selection(
   if (toc(start_time) > settings.time_limit) {
     settings.log.debug("Time limit reached\n");
     if (use_pdlp) {
-      concurrent_halt.store(1);
+      pdlp_concurrent_halt.store(1);
 #pragma omp taskwait  // Wait for the batch PDLP task to finish
     }
     return branch_var;
@@ -1748,6 +1750,8 @@ i_t pseudo_costs_t<i_t, f_t>::reliable_variable_selection(
 #pragma omp taskloop if (num_tasks > 1) priority(CUOPT_HIGH_TASK_PRIORITY) \
   num_tasks(num_tasks) default(shared)
     for (i_t i = 0; i < num_candidates; ++i) {
+      if (*concurrent_halt_ == 1) continue;  // OpenMP does not allow to break out of the loop
+
       auto [score, j] = unreliable_list[i];
 
       if (toc(start_time) > settings.time_limit) { continue; }
@@ -1768,6 +1772,7 @@ i_t pseudo_costs_t<i_t, f_t>::reliable_variable_selection(
                                                      worker->basis_factors,
                                                      worker->basic_list,
                                                      worker->nonbasic_list,
+                                                     concurrent_halt_,
                                                      j,
                                                      worker->leaf_problem.lower[j],
                                                      std::floor(leaf_solution.x[j]),
@@ -1813,6 +1818,7 @@ i_t pseudo_costs_t<i_t, f_t>::reliable_variable_selection(
                                                      worker->basis_factors,
                                                      worker->basic_list,
                                                      worker->nonbasic_list,
+                                                     concurrent_halt_,
                                                      j,
                                                      std::ceil(leaf_solution.x[j]),
                                                      worker->leaf_problem.upper[j],
@@ -1849,7 +1855,7 @@ i_t pseudo_costs_t<i_t, f_t>::reliable_variable_selection(
       score_mutex.unlock();
     }
 
-    concurrent_halt.store(1);
+    pdlp_concurrent_halt.store(1);
   }
 
   f_t dual_simplex_elapsed = toc(dual_simplex_start_time);
diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp
index 6fd7dc5bb6..4c72f16d3d 100644
--- a/cpp/src/branch_and_bound/pseudo_costs.hpp
+++ b/cpp/src/branch_and_bound/pseudo_costs.hpp
@@ -214,6 +214,7 @@ class pseudo_costs_t {
 
   std::shared_ptr<csc_matrix_t<i_t, f_t>> AT;  // Transpose of the constraint matrix A
   std::shared_ptr<batch_pdlp_warm_cache_t<i_t, f_t>> pdlp_warm_cache;
+  std::atomic<int>* concurrent_halt_;
 
   reliability_branching_settings_t<i_t, f_t> reliability_branching_settings;
   simplex_solver_settings_t<i_t, f_t> settings;
diff --git a/cpp/src/branch_and_bound/search_tree.hpp b/cpp/src/branch_and_bound/search_tree.hpp
new file mode 100644
index 0000000000..0f5fc2258e
--- /dev/null
+++ b/cpp/src/branch_and_bound/search_tree.hpp
@@ -0,0 +1,165 @@
+/* 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/mip_node.hpp>
+
+#include <utilities/omp_helpers.hpp>
+
+#include <cmath>
+#include <cstdio>
+#include <memory>
+#include <mutex>
+#include <vector>
+
+namespace cuopt::linear_programming::dual_simplex {
+
+template <typename i_t, typename f_t>
+class search_tree_t {
+ public:
+  search_tree_t() = default;
+
+  search_tree_t(mip_node_t<i_t, f_t>&& node) : search_tree_t() { root = std::move(node); }
+
+  ~search_tree_t() { clean(); }
+
+  void update(mip_node_t<i_t, f_t>* node_ptr, node_status_t status)
+  {
+    std::lock_guard lock(mutex);
+
+    --num_open_nodes;
+    if (status == node_status_t::HAS_CHILDREN) {
+      ++num_inner_nodes;
+    } else {
+      ++num_final_nodes;
+      progress += std::ldexp(f_t(1), -node_ptr->depth);
+    }
+
+    std::vector<mip_node_t<i_t, f_t>*> stack;
+    node_ptr->set_status(status, stack);
+    remove_fathomed_nodes(stack);
+  }
+
+  void branch(mip_node_t<i_t, f_t>* parent_node,
+              const i_t branch_var,
+              const f_t fractional_val,
+              const i_t integer_infeasible,
+              const std::vector<variable_status_t>& parent_vstatus,
+              const lp_problem_t<i_t, f_t>& original_lp,
+              logger_t& log)
+  {
+    i_t id = num_nodes.fetch_add(2);
+
+    auto down_child = std::make_unique<mip_node_t<i_t, f_t>>(original_lp,
+                                                             parent_node,
+                                                             ++id,
+                                                             branch_var,
+                                                             branch_direction_t::DOWN,
+                                                             fractional_val,
+                                                             integer_infeasible,
+                                                             parent_vstatus);
+    graphviz_edge(log,
+                  parent_node,
+                  down_child.get(),
+                  branch_var,
+                  branch_direction_t::DOWN,
+                  std::floor(fractional_val));
+
+    auto up_child = std::make_unique<mip_node_t<i_t, f_t>>(original_lp,
+                                                           parent_node,
+                                                           ++id,
+                                                           branch_var,
+                                                           branch_direction_t::UP,
+                                                           fractional_val,
+                                                           integer_infeasible,
+                                                           parent_vstatus);
+
+    graphviz_edge(log,
+                  parent_node,
+                  up_child.get(),
+                  branch_var,
+                  branch_direction_t::UP,
+                  std::ceil(fractional_val));
+
+    assert(parent_vstatus.size() == original_lp.num_cols);
+    parent_node->add_children(std::move(down_child),
+                              std::move(up_child));  // child pointers moved into the tree
+    num_open_nodes += 2;
+  }
+
+  static void graphviz_node(logger_t& log,
+                            const mip_node_t<i_t, f_t>* node_ptr,
+                            const std::string label,
+                            const f_t val)
+  {
+    if (write_graphviz) {
+      log.printf("Node%d [label=\"%s %.16e\"]\n", node_ptr->node_id, label.c_str(), val);
+    }
+  }
+
+  static void graphviz_edge(logger_t& log,
+                            const mip_node_t<i_t, f_t>* origin_ptr,
+                            const mip_node_t<i_t, f_t>* dest_ptr,
+                            const i_t branch_var,
+                            branch_direction_t branch_dir,
+                            const f_t bound)
+  {
+    if (write_graphviz) {
+      log.printf("Node%d -> Node%d [label=\"x%d %s %e\"]\n",
+                 origin_ptr->node_id,
+                 dest_ptr->node_id,
+                 branch_var,
+                 branch_dir == branch_direction_t::DOWN ? "<=" : ">=",
+                 bound);
+    }
+  }
+
+  // Clean the tree using a depth first scheme
+  void clean()
+  {
+    std::vector<std::unique_ptr<mip_node_t<i_t, f_t>>> stack;
+
+    if (root.children[0]) stack.push_back(std::move(root.children[0]));
+    if (root.children[1]) stack.push_back(std::move(root.children[1]));
+
+    while (!stack.empty()) {
+      auto node = std::move(stack.back());
+      stack.pop_back();
+      if (node->children[0]) stack.push_back(std::move(node->children[0]));
+      if (node->children[1]) stack.push_back(std::move(node->children[1]));
+      // Implicitly call destructor for `node`
+    }
+
+    num_nodes       = 0;
+    num_open_nodes  = 0;
+    num_final_nodes = 0;
+    num_inner_nodes = 0;
+    progress        = 0;
+  }
+
+  mip_node_t<i_t, f_t> root;
+  omp_mutex_t mutex;
+  omp_atomic_t<i_t> num_nodes = 0;
+
+  // Number of nodes that still needs to be explored
+  omp_atomic_t<i_t> num_open_nodes = 0;
+
+  // Number of integer feasible, infeasible or fathomed nodes
+  omp_atomic_t<i_t> num_final_nodes = 0;
+
+  // Number of inner nodes
+  omp_atomic_t<i_t> num_inner_nodes = 0;
+
+  // Track the solver progress based on how much the tree was explored
+  // using the tree weight metric
+  omp_atomic_t<f_t> progress = 0;
+
+  static constexpr bool write_graphviz = false;
+};
+
+}  // namespace cuopt::linear_programming::dual_simplex
diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp
index bcce5ece2a..231e1e0902 100644
--- a/cpp/src/branch_and_bound/symmetry.hpp
+++ b/cpp/src/branch_and_bound/symmetry.hpp
@@ -319,7 +319,7 @@ class orbital_fixing_t {
       node = node->parent;
     }
 
-    surviving_generators_.resize(max_generators_);
+    surviving_generators_.resize(symmetry->num_generators);
     std::iota(surviving_generators_.begin(), surviving_generators_.end(), 0);
 
     // Seed cumulative fixings from the parent's stored orbital fixings.
diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp
index 62bd0cde6b..739d4cc6b2 100644
--- a/cpp/src/branch_and_bound/worker.hpp
+++ b/cpp/src/branch_and_bound/worker.hpp
@@ -24,17 +24,25 @@ namespace cuopt::linear_programming::dual_simplex {
 
 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;
+  f_t start_time = 0.0;
+
+  omp_atomic_t<f_t> total_lp_solve_time = 0.0;
+  omp_atomic_t<int64_t> total_lp_iters  = 0;
+
+  omp_atomic_t<int64_t> nodes_explored       = 0;
+  omp_atomic_t<int64_t> total_nodes_explored = 0;
+  omp_atomic_t<int64_t> nodes_unexplored     = 0;
   // Tracks the number of nodes being solved by the workers at a given time
   omp_atomic_t<i_t> nodes_being_solved = 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;
 
+  i_t restart_nodes_at_last_check    = 0;
+  f_t restart_progress_at_last_check = 0;
+  f_t restart_gap_at_last_check      = 0;
+  i_t restart_large_tree_count       = 0;
+
   omp_atomic_t<int64_t> orbital_fixing_nodes              = 0;
   omp_atomic_t<int64_t> orbital_fixings_applied           = 0;
   omp_atomic_t<int64_t> orbital_conflict_nodes            = 0;
@@ -203,6 +211,16 @@ class bfs_worker_t : public branch_and_bound_worker_t<i_t, f_t> {
     }
   }
 
+  void reset_state()
+  {
+    node_queue.clear();
+    total_max_diving_workers    = 0;
+    total_active_diving_workers = 0;
+    active_diving_workers.fill(0);
+    this->is_active   = false;
+    this->lower_bound = -std::numeric_limits<f_t>::infinity();
+  }
+
   // The worker-local node heap.
   node_queue_t<i_t, f_t> node_queue;
 
@@ -249,6 +267,13 @@ class diving_worker_t : public branch_and_bound_worker_t<i_t, f_t> {
 
   f_t get_lower_bound() { return this->lower_bound; }
 
+  void reset_state()
+  {
+    this->is_active   = false;
+    this->lower_bound = -std::numeric_limits<f_t>::infinity();
+    bfs_worker        = nullptr;
+  }
+
   mip_node_t<i_t, f_t> start_node;
 
   // The best-first worker that is associated with this diving worker. Used for controlling the
diff --git a/cpp/src/branch_and_bound/worker_pool.hpp b/cpp/src/branch_and_bound/worker_pool.hpp
index d75794bbf0..858d3949be 100644
--- a/cpp/src/branch_and_bound/worker_pool.hpp
+++ b/cpp/src/branch_and_bound/worker_pool.hpp
@@ -91,6 +91,17 @@ class worker_pool_t {
   i_t num_idle() const { return num_idle_workers_; }
   i_t size() const { return workers_.size(); }
 
+  void reset()
+  {
+    std::lock_guard lock(mutex_);
+    num_idle_workers_ = workers_.size();
+    idle_workers_.clear_resize(workers_.size());
+    for (i_t i = 0; i < workers_.size(); ++i) {
+      workers_[i]->reset_state();
+      idle_workers_.push_back(i);
+    }
+  }
+
  private:
   std::vector<std::unique_ptr<WorkerType>> workers_;
   bool is_initialized_ = false;
diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp
index 286ac7364f..ef94f6db72 100644
--- a/cpp/src/dual_simplex/simplex_solver_settings.hpp
+++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp
@@ -207,6 +207,14 @@ struct simplex_solver_settings_t {
   i_t bnb_nodes_per_steal;
   i_t bnb_max_steal_attempts;
 
+  i_t restart_min_nodes                  = 1000;
+  i_t restart_min_estimates              = 10;
+  f_t restart_threshold_grow_per_node    = 0.001;
+  f_t restart_threshold_grow_per_restart = 1.5;
+  i_t restart_tree_size_factor           = 50;
+  i_t restart_check_freq                 = 100;
+  i_t max_restarts                       = 50;
+
   // Settings for the reliability branching.
   // - -1: automatic
   // - 0: disable (use pseudocost branching instead)
diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp
index d81265358a..124567115a 100644
--- a/cpp/src/dual_simplex/solve.cpp
+++ b/cpp/src/dual_simplex/solve.cpp
@@ -725,8 +725,10 @@ i_t solve(const user_problem_t<i_t, f_t>& problem,
 {
   i_t status;
   if (is_mip(problem) && !settings.relaxation) {
+    std::atomic<int> restart_concurrent_halt;
     probing_implied_bound_t<i_t, f_t> empty_probing(problem.num_cols);
-    branch_and_bound_t branch_and_bound(problem, settings, tic(), empty_probing);
+    branch_and_bound_t branch_and_bound(
+      problem, settings, tic(), &restart_concurrent_halt, empty_probing);
     mip_solution_t<i_t, f_t> mip_solution(problem.num_cols);
     mip_status_t mip_status = branch_and_bound.solve(mip_solution);
     if (mip_status == mip_status_t::OPTIMAL) {
@@ -766,7 +768,9 @@ i_t solve_mip_with_guess(const user_problem_t<i_t, f_t>& problem,
   i_t status;
   if (is_mip(problem)) {
     probing_implied_bound_t<i_t, f_t> empty_probing(problem.num_cols);
-    branch_and_bound_t branch_and_bound(problem, settings, tic(), empty_probing);
+    std::atomic<int> restart_concurrent_halt;
+    branch_and_bound_t branch_and_bound(
+      problem, settings, tic(), &restart_concurrent_halt, empty_probing);
     branch_and_bound.set_initial_guess(guess);
     mip_status_t mip_status = branch_and_bound.solve(solution);
     if (mip_status == mip_status_t::OPTIMAL) {
diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu
index e1318edf4d..d70ee86520 100644
--- a/cpp/src/mip_heuristics/diversity/lns/rins.cu
+++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu
@@ -214,6 +214,7 @@ void rins_t<i_t, f_t>::run_rins()
   fj_solution.copy_new_assignment(cuopt::host_copy(fixed_assignment, rins_handle.get_stream()));
   std::vector<f_t> default_weights(fixed_problem.n_constraints, 1.);
 
+  std::atomic<bool> fj_halt_flag = false;
   std::unique_ptr<fj_cpu_climber_t<i_t, f_t>> fj_cpu =
     fj.create_cpu_climber(fj_solution,
                           default_weights,
@@ -222,7 +223,8 @@ void rins_t<i_t, f_t>::run_rins()
                           context.preempt_heuristic_solver_,
                           fj_settings_t{},
                           true);
-  fj_cpu->log_prefix = "[RINS] ";
+  fj_cpu->log_prefix      = "[RINS] ";
+  fj_cpu->preemption_flag = &fj_halt_flag;
 
   CUOPT_LOG_DEBUG("Launching CPUFJ (RINS) task");
 #pragma omp task shared(fj_cpu) firstprivate(time_limit) \
@@ -263,8 +265,11 @@ void rins_t<i_t, f_t>::run_rins()
     rins_solution_queue.push_back(solution);
   };
   dual_simplex::probing_implied_bound_t<i_t, f_t> empty_probing(branch_and_bound_problem.num_cols);
-  dual_simplex::branch_and_bound_t<i_t, f_t> branch_and_bound(
-    branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic(), empty_probing);
+  dual_simplex::branch_and_bound_t<i_t, f_t> branch_and_bound(branch_and_bound_problem,
+                                                              branch_and_bound_settings,
+                                                              dual_simplex::tic(),
+                                                              &context.restart_concurrent_halt,
+                                                              empty_probing);
   branch_and_bound.set_initial_guess(cuopt::host_copy(fixed_assignment, rins_handle.get_stream()));
   branch_and_bound_status = branch_and_bound.solve(branch_and_bound_solution);
 
@@ -298,6 +303,7 @@ void rins_t<i_t, f_t>::run_rins()
                           static_cast<f_t>(context.settings.heuristic_params.rins_max_time_limit));
   }
 
+  fj_halt_flag = true;
 #pragma omp taskwait  // Wait for the CPU FJ (RINS) to finish
   CUOPT_LOG_DEBUG("CPUFJ (RINS) task was stopped");
 
diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh
index 1d0b9245d7..a9667cf0a3 100644
--- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh
+++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh
@@ -122,8 +122,12 @@ class sub_mip_recombiner_t : public recombiner_t<i_t, f_t> {
       branch_and_bound_settings.log.log = false;
       dual_simplex::probing_implied_bound_t<i_t, f_t> empty_probing(
         branch_and_bound_problem.num_cols);
-      dual_simplex::branch_and_bound_t<i_t, f_t> branch_and_bound(
-        branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic(), empty_probing);
+      std::atomic<int> concurrent_halt = 0;
+      dual_simplex::branch_and_bound_t<i_t, f_t> branch_and_bound(branch_and_bound_problem,
+                                                                  branch_and_bound_settings,
+                                                                  dual_simplex::tic(),
+                                                                  &concurrent_halt,
+                                                                  empty_probing);
       branch_and_bound_status = branch_and_bound.solve(branch_and_bound_solution);
       if (solution_vector.size() > 0) {
         cuopt_assert(fixed_assignment.size() == branch_and_bound_solution.x.size(),
diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu
index a537426457..0382d700bf 100644
--- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu
+++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu
@@ -1614,7 +1614,7 @@ void cpufj_solve(fj_cpu_climber_t<i_t, f_t>* fj_cpu, f_t in_time_limit, double w
   fj_cpu->prev_best_objective   = fj_cpu->h_best_objective;
   fj_cpu->iterations_since_best = 0;
 
-  while (!fj_cpu->halted && !fj_cpu->preemption_flag.load()) {
+  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<f_t>::infinity() &&
@@ -1823,9 +1823,9 @@ template <typename i_t, typename f_t>
 void stop_fj_cpu_task(fj_cpu_task_t<i_t, f_t>& task)
 {
   if (task.fj_cpu) {
-    auto& fj_cpu           = *task.fj_cpu;
-    fj_cpu.preemption_flag = true;
-    fj_cpu.halted          = true;
+    auto& fj_cpu = *task.fj_cpu;
+    fj_cpu.preemption_flag->store(true);
+    fj_cpu.halted = true;
   }
 }
 
diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh
index cdf3a2f58a..47383ac759 100644
--- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh
+++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh
@@ -24,7 +24,7 @@ namespace cuopt::linear_programming::detail {
 // Maintaining a single source of truth for all members would be nice
 template <typename i_t, typename f_t>
 struct fj_cpu_climber_t {
-  fj_cpu_climber_t(std::atomic<bool>& preemption_flag) : preemption_flag(preemption_flag)
+  fj_cpu_climber_t(std::atomic<bool>& preemption_flag) : preemption_flag(&preemption_flag)
   {
 #define ADD_INSTRUMENTED(var) \
   std::make_pair(#var, std::ref(static_cast<memory_instrumentation_base_t&>(var)))
@@ -187,7 +187,7 @@ struct fj_cpu_climber_t {
   // Memory instrumentation aggregator
   instrumentation_aggregator_t memory_aggregator;
   // TODO atomic ref? c++20
-  std::atomic<bool>& preemption_flag;
+  std::atomic<bool>* preemption_flag;
 };
 
 template <typename i_t, typename f_t>
diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu
index 858e865bc3..23b72675cc 100644
--- a/cpp/src/mip_heuristics/solver.cu
+++ b/cpp/src/mip_heuristics/solver.cu
@@ -430,6 +430,7 @@ solution_t<i_t, f_t> mip_solver_t<i_t, f_t>::run_solver()
       branch_and_bound_problem,
       branch_and_bound_settings,
       timer_.get_tic_start(),
+      &context.restart_concurrent_halt,
       probing_implied_bound,
       context.problem_ptr->clique_table,
       context.symmetry.get());
diff --git a/cpp/src/mip_heuristics/solver_context.cuh b/cpp/src/mip_heuristics/solver_context.cuh
index 739f7d130a..6a5b84c636 100644
--- a/cpp/src/mip_heuristics/solver_context.cuh
+++ b/cpp/src/mip_heuristics/solver_context.cuh
@@ -77,6 +77,8 @@ struct mip_solver_context_t {
 
   // Symmetry information for orbital fixing during B&B. Null if no exploitable symmetry.
   std::unique_ptr<dual_simplex::mip_symmetry_t<i_t, f_t>> symmetry;
+
+  std::atomic<int> restart_concurrent_halt{0};
 };
 
 }  // namespace cuopt::linear_programming::detail