diff --git a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp
index b56d4f409a..2bbc9871e3 100644
--- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp
+++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp
@@ -26,6 +26,7 @@
 namespace cuopt {
 namespace CUOPT_EXPORT mathematical_optimization {
 
+template <typename i_t, typename f_t>
 class barrier_cache_t;
 
 // Forward declare solver_settings_t for friend class
@@ -384,7 +385,7 @@ class pdlp_solver_settings_t {
   /** When true, the first GPU barrier/QCQP solve retains cache state for later reuse. */
   bool sequence_solve{false};
   /** Non-owning cache pointer set by ``call_solve`` for barrier cache reuse. */
-  barrier_cache_t* barrier_cache{nullptr};
+  barrier_cache_t<i_t, f_t>* barrier_cache{nullptr};
 
  private:
   /** Initial primal solution */
diff --git a/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp b/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp
index 7993783b91..4c5490db2f 100644
--- a/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp
+++ b/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp
@@ -30,6 +30,7 @@ void apply_barrier_rhs(iteration_data_t<int, double>& data, double const* barrie
 namespace cuopt {
 namespace CUOPT_EXPORT mathematical_optimization {
 
+template <typename i_t, typename f_t>
 struct barrier_transform_t;
 
 /**
@@ -39,6 +40,7 @@ struct barrier_transform_t;
  * The update APIs crush new user data into that workspace and marks the cache dirty so the
  * next solve reuses it (skip convert/presolve/scaling).
  */
+template <typename i_t, typename f_t>
 class barrier_cache_t {
  public:
   static std::unique_ptr<barrier_cache_t> create(unsigned stream_flags);
@@ -56,16 +58,16 @@ class barrier_cache_t {
   /**
    * @brief Take ownership of barrier iteration workspace. @p data may be null (clears).
    */
-  void store_iteration_data(barrier::iteration_data_t<int, double>* data);
+  void store_iteration_data(barrier::iteration_data_t<i_t, f_t>* data);
 
   /**
    * @brief Release ownership of cached iteration workspace; caller must delete or wrap it.
    */
-  barrier::iteration_data_t<int, double>* release_iteration_data();
+  barrier::iteration_data_t<i_t, f_t>* release_iteration_data();
 
-  void store_transform(std::unique_ptr<barrier_transform_t> transform);
-  [[nodiscard]] barrier_transform_t* transform();
-  [[nodiscard]] barrier_transform_t const* transform() const;
+  void store_transform(std::unique_ptr<barrier_transform_t<i_t, f_t>> transform);
+  [[nodiscard]] barrier_transform_t<i_t, f_t>* transform();
+  [[nodiscard]] barrier_transform_t<i_t, f_t> const* transform() const;
   /** True when an update API has staged new data that the next solve should reuse. */
   [[nodiscard]] bool dirty() const;
   void mark_clean();
@@ -77,13 +79,13 @@ class barrier_cache_t {
    * Crush the input linear objective into cached iteration_data_t.c / d_c_ and mark dirty.
    * Requires a stored transform and iteration_data from an Optimal solve.
    */
-  void update_linear_objective(double const* c, int n);
+  void update_linear_objective(f_t const* c, i_t n);
 
   /**
    * Crush the input constraint RHS into cached iteration_data_t.b / d_b_ and mark dirty.
    * Requires a stored transform and iteration_data from an Optimal solve.
    */
-  void update_rhs(double const* b, int m);
+  void update_rhs(f_t const* b, i_t m);
 
  private:
   barrier_cache_t(std::unique_ptr<rmm::cuda_stream> stream, std::unique_ptr<raft::handle_t> handle);
diff --git a/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp b/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp
index 57042c1e05..705a1e064c 100644
--- a/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp
+++ b/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp
@@ -54,9 +54,9 @@ mathematical_optimization::mip_solution_interface_t<int, double>* call_solve_mip
 std::unique_ptr<solver_ret_t> call_solve(
   cuopt::mathematical_optimization::io::data_model_view_t<int, double>*,
   mathematical_optimization::solver_settings_t<int, double>*,
-  unsigned int flags                                   = cudaStreamNonBlocking,
-  bool is_batch_mode                                   = false,
-  mathematical_optimization::barrier_cache_t* cache_in = nullptr);
+  unsigned int flags                                                = cudaStreamNonBlocking,
+  bool is_batch_mode                                                = false,
+  mathematical_optimization::barrier_cache_t<int, double>* cache_in = nullptr);
 
 std::pair<std::vector<std::unique_ptr<solver_ret_t>>, double> solve_batch_remote(
   std::vector<cuopt::mathematical_optimization::io::data_model_view_t<int, double>*>,
diff --git a/cpp/include/cuopt/mathematical_optimization/utilities/cython_types.hpp b/cpp/include/cuopt/mathematical_optimization/utilities/cython_types.hpp
index 7555e8b19b..73cb351c71 100644
--- a/cpp/include/cuopt/mathematical_optimization/utilities/cython_types.hpp
+++ b/cpp/include/cuopt/mathematical_optimization/utilities/cython_types.hpp
@@ -22,6 +22,7 @@ namespace cuopt {
 namespace CUOPT_EXPORT mathematical_optimization {
 // Forward declared, not included: these structs are also compiled into cuopt_client, which
 // is CPU-only and cannot link the GPU-side barrier_cache_t destructor.
+template <typename i_t, typename f_t>
 class barrier_cache_t;
 }  // namespace CUOPT_EXPORT mathematical_optimization
 
@@ -101,7 +102,7 @@ struct linear_programming_ret_t {
   /** GPU barrier cache (stream + handle + iteration workspace), non-owning. call_solve hands
    * ownership to the caller, which wraps it in a Python capsule and deletes it there.
    */
-  mathematical_optimization::barrier_cache_t* barrier_cache{nullptr};
+  mathematical_optimization::barrier_cache_t<int, double>* barrier_cache{nullptr};
 
   bool is_gpu() const { return std::holds_alternative<lp_gpu_ptr>(solutions_); }
 };
diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu
index 28b5b56c70..5a2eca9ee7 100644
--- a/cpp/src/barrier/barrier.cu
+++ b/cpp/src/barrier/barrier.cu
@@ -932,7 +932,7 @@ class iteration_data_t {
         {
           raft::common::nvtx::range form_scope("Barrier: LP Data: form augmented");
           // Build the sparsity pattern of the augmented system
-          form_augmented(true);
+          form_augmented(augmented_form_t::build);
         }
         if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; }
         symbolic_status = chol->analyze(device_augmented);
@@ -948,9 +948,10 @@ class iteration_data_t {
   }
 
   // Attach this solve's settings and rewind iterate-dependent state so barrier can
-  // start with the new c / b. A and Q are unchanged; the previous solve
-  // left D and the KKT values at its last iterate. Reuse is QP-only (no cones),
-  // so form_*(false) updates values in the existing CSR; no symbolic rebuild.
+  // start with the new c / b. A and Q are unchanged. The previous solve left D
+  // and the KKT values at its last iterate. reset_for_reuse (or form_adat(false)) rewrites
+  // values in the existing CSR. The cone block is put back to the initial diagonal a cold
+  // start factorizes; Nesterov-Todd scaling is recomputed from the new point.
   bool reset_iterate_state(const simplex_solver_settings_t<i_t, f_t>& settings)
   {
     if (chol == nullptr || symbolic_status != 0) { return false; }
@@ -987,7 +988,7 @@ class iteration_data_t {
     }
 
     if (use_augmented) {
-      form_augmented(false);
+      form_augmented(augmented_form_t::reset_for_reuse);
     } else {
       form_adat(false);
     }
@@ -1063,7 +1064,12 @@ class iteration_data_t {
     return degree;
   }
 
-  void form_augmented(bool first_call = false)
+  // build: first call, device CSR and metadata.
+  // reset_for_reuse: values only; cone block back to the cold-start matrix.
+  // update: values only; cone block from the current Nesterov-Todd scaling.
+  enum class augmented_form_t { build, reset_for_reuse, update };
+
+  void form_augmented(augmented_form_t mode = augmented_form_t::update)
   {
     i_t n    = A.n;
     i_t m    = A.m;
@@ -1074,7 +1080,7 @@ class iteration_data_t {
     const i_t p            = augmented_expansion_count();
     i_t factorization_size = augmented_system_size(n, m);
 
-    if (first_call) {
+    if (mode == augmented_form_t::build) {
       raft::common::nvtx::range scope("Barrier: augmented: device CSR build");
 
       const size_t n_sparse_cone_entries =
@@ -1165,7 +1171,40 @@ class iteration_data_t {
                          });
       RAFT_CHECK_CUDA(handle_ptr->get_stream().get());
 
-      if (has_soc) {
+      if (has_soc && mode == augmented_form_t::reset_for_reuse) {
+        // Cold initial_point factorizes this diagonal, then the first Newton step
+        // rebuilds the Nesterov-Todd Hessian. Zero w and eta so a dense block
+        // scatter and the matrix-free product both see that same initial matrix.
+        auto stream = handle_ptr->get_stream();
+        thrust::fill(rmm::exec_policy(stream), cones().w.begin(), cones().w.end(), f_t(0));
+        thrust::fill(rmm::exec_policy(stream), cones().eta.begin(), cones().eta.end(), f_t(0));
+        if (cones().has_sparse_cones()) {
+          restore_initial_sparse_cone_block(cones(),
+                                            device_augmented.x,
+                                            cone_kkt_data_.sparse_Hs_diag,
+                                            cone_kkt_data_.sparse_hessian_diag,
+                                            cone_kkt_data_.sparse_hessian_Q,
+                                            cone_kkt_data_.sparse_exp_v_col,
+                                            cone_kkt_data_.sparse_exp_u_col,
+                                            cone_kkt_data_.sparse_exp_v_row,
+                                            cone_kkt_data_.sparse_exp_u_row,
+                                            cone_kkt_data_.sparse_expansion_D,
+                                            stream,
+                                            dual_perturb);
+          RAFT_CHECK_CUDA(stream.get());
+        }
+        if (cones().n_dense_cones() > 0) {
+          scatter_dense_hessian_into_augmented(cones(),
+                                               device_augmented.x,
+                                               cone_kkt_data_.cone_csr_indices,
+                                               cone_kkt_data_.cone_Q_values,
+                                               cone_kkt_data_.dense_block_offsets,
+                                               cone_kkt_data_.dense_cone_ids,
+                                               stream,
+                                               dual_perturb);
+          RAFT_CHECK_CUDA(stream.get());
+        }
+      } else if (has_soc) {
         if (cones().has_sparse_cones()) {
           scatter_sparse_hessian_into_augmented(cones(),
                                                 device_augmented.x,
@@ -4884,7 +4923,7 @@ template <typename i_t, typename f_t>
 lp_status_t barrier_solver_t<i_t, f_t>::solve_with_cache(
   f_t start_time,
   lp_solution_t<i_t, f_t>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache)
+  cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache)
 {
   settings.log.printf("Barrier solver started at %.2f seconds\n", toc(start_time));
   try {
@@ -4941,7 +4980,7 @@ template <typename i_t, typename f_t>
 lp_status_t barrier_solver_t<i_t, f_t>::solve(
   f_t start_time,
   lp_solution_t<i_t, f_t>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache)
+  cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache)
 {
   settings.log.printf("Barrier solver started at %.2f seconds\n", toc(start_time));
   try {
@@ -5010,7 +5049,7 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(
 
 // Optimal: persist iteration_data_t on the cache. Otherwise drop it.
 template <typename i_t, typename f_t>
-lp_status_t store_or_clear_cache(cuopt::mathematical_optimization::barrier_cache_t* cache,
+lp_status_t store_or_clear_cache(cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache,
                                  std::unique_ptr<iteration_data_t<i_t, f_t>>& owned_data,
                                  lp_status_t status)
 {
diff --git a/cpp/src/barrier/barrier.hpp b/cpp/src/barrier/barrier.hpp
index 9f7dbb2e2f..6de3d2b1be 100644
--- a/cpp/src/barrier/barrier.hpp
+++ b/cpp/src/barrier/barrier.hpp
@@ -24,6 +24,7 @@
 #include <utility>
 
 namespace cuopt::mathematical_optimization {
+template <typename i_t, typename f_t>
 class barrier_cache_t;
 }
 
@@ -65,15 +66,17 @@ class barrier_solver_t {
                    const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
                    device_csc_matrix_ptr_t<i_t, f_t> device_A = nullptr,
                    device_csc_matrix_ptr_t<i_t, f_t> device_Q = nullptr);
-  simplex::lp_status_t solve(f_t start_time,
-                             simplex::lp_solution_t<i_t, f_t>& solution,
-                             cuopt::mathematical_optimization::barrier_cache_t* cache = nullptr);
+  simplex::lp_status_t solve(
+    f_t start_time,
+    simplex::lp_solution_t<i_t, f_t>& solution,
+    cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache = nullptr);
   // Cache reuse: cached iteration_data_t already has the updated linear objective.
   // Reset iterate state, compute a new initial point, run barrier. Same status/solution contract as
   // solve().
-  simplex::lp_status_t solve_with_cache(f_t start_time,
-                                        simplex::lp_solution_t<i_t, f_t>& solution,
-                                        cuopt::mathematical_optimization::barrier_cache_t* cache);
+  simplex::lp_status_t solve_with_cache(
+    f_t start_time,
+    simplex::lp_solution_t<i_t, f_t>& solution,
+    cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache);
 
  private:
   simplex::lp_status_t barrier_advanced_solve(f_t start_time,
diff --git a/cpp/src/barrier/barrier_cache.cu b/cpp/src/barrier/barrier_cache.cu
index 1d80b8d389..52ed0f7289 100644
--- a/cpp/src/barrier/barrier_cache.cu
+++ b/cpp/src/barrier/barrier_cache.cu
@@ -19,12 +19,9 @@
 
 namespace cuopt::mathematical_optimization {
 
-using barrier_iteration_data_t = barrier::iteration_data_t<int, double>;
-using barrier_iteration_data_ptr =
-  std::unique_ptr<barrier_iteration_data_t, void (*)(barrier_iteration_data_t*)>;
-
-static void require_cache(barrier_transform_t const* transform,
-                          barrier_iteration_data_t const* data,
+template <typename i_t, typename f_t>
+static void require_cache(barrier_transform_t<i_t, f_t> const* transform,
+                          barrier::iteration_data_t<i_t, f_t> const* data,
                           char const* api)
 {
   cuopt_expects(transform != nullptr,
@@ -38,7 +35,8 @@ static void require_cache(barrier_transform_t const* transform,
 }
 
 // Re-adds the first solve's barrier-minus-crush shift so the update lands in the presolved model.
-static void add_shift(std::vector<double>& crushed, std::vector<double> const& shift)
+template <typename f_t>
+static void add_shift(std::vector<f_t>& crushed, std::vector<f_t> const& shift)
 {
   cuopt_expects(shift.size() == crushed.size(),
                 error_type_t::ValidationError,
@@ -48,7 +46,11 @@ static void add_shift(std::vector<double>& crushed, std::vector<double> const& s
   }
 }
 
-struct barrier_cache_t::impl {
+template <typename i_t, typename f_t>
+struct barrier_cache_t<i_t, f_t>::impl {
+  using iteration_data_t   = barrier::iteration_data_t<i_t, f_t>;
+  using iteration_data_ptr = std::unique_ptr<iteration_data_t, void (*)(iteration_data_t*)>;
+
   impl(std::unique_ptr<rmm::cuda_stream> stream_in, std::unique_ptr<raft::handle_t> handle_in)
     : stream(std::move(stream_in)),
       handle(std::move(handle_in)),
@@ -59,91 +61,126 @@ struct barrier_cache_t::impl {
   std::unique_ptr<rmm::cuda_stream> stream;
   std::unique_ptr<raft::handle_t> handle;
   // Destroy iteration_data before transform: it may const-ref A/Q stored on the transform.
-  std::unique_ptr<barrier_transform_t> transform;
-  barrier_iteration_data_ptr iteration_data;
+  std::unique_ptr<barrier_transform_t<i_t, f_t>> transform;
+  iteration_data_ptr iteration_data;
   bool linear_objective_dirty{false};
   bool rhs_dirty{false};
   bool rhs_infeasible{false};
 };
 
-barrier_cache_t::barrier_cache_t(std::unique_ptr<rmm::cuda_stream> stream,
-                                 std::unique_ptr<raft::handle_t> handle)
+template <typename i_t, typename f_t>
+barrier_cache_t<i_t, f_t>::barrier_cache_t(std::unique_ptr<rmm::cuda_stream> stream,
+                                           std::unique_ptr<raft::handle_t> handle)
   : impl_(std::make_unique<impl>(std::move(stream), std::move(handle)))
 {
 }
 
-barrier_cache_t::~barrier_cache_t() = default;
+template <typename i_t, typename f_t>
+barrier_cache_t<i_t, f_t>::~barrier_cache_t() = default;
 
-barrier_cache_t::barrier_cache_t(barrier_cache_t&&) noexcept            = default;
-barrier_cache_t& barrier_cache_t::operator=(barrier_cache_t&&) noexcept = default;
+template <typename i_t, typename f_t>
+barrier_cache_t<i_t, f_t>::barrier_cache_t(barrier_cache_t&&) noexcept = default;
 
-std::unique_ptr<barrier_cache_t> barrier_cache_t::create(unsigned stream_flags)
+template <typename i_t, typename f_t>
+barrier_cache_t<i_t, f_t>& barrier_cache_t<i_t, f_t>::operator=(barrier_cache_t&&) noexcept =
+  default;
+
+template <typename i_t, typename f_t>
+std::unique_ptr<barrier_cache_t<i_t, f_t>> barrier_cache_t<i_t, f_t>::create(unsigned stream_flags)
 {
   auto stream =
     std::make_unique<rmm::cuda_stream>(static_cast<rmm::cuda_stream::flags>(stream_flags));
   auto handle = std::make_unique<raft::handle_t>(*stream);
-  return std::unique_ptr<barrier_cache_t>(
+  return std::unique_ptr<barrier_cache_t<i_t, f_t>>(
     new barrier_cache_t(std::move(stream), std::move(handle)));
 }
 
-raft::handle_t* barrier_cache_t::handle_ptr() { return impl_->handle.get(); }
+template <typename i_t, typename f_t>
+raft::handle_t* barrier_cache_t<i_t, f_t>::handle_ptr()
+{
+  return impl_->handle.get();
+}
 
-raft::handle_t const* barrier_cache_t::handle_ptr() const { return impl_->handle.get(); }
+template <typename i_t, typename f_t>
+raft::handle_t const* barrier_cache_t<i_t, f_t>::handle_ptr() const
+{
+  return impl_->handle.get();
+}
 
-void barrier_cache_t::clear()
+template <typename i_t, typename f_t>
+void barrier_cache_t<i_t, f_t>::clear()
 {
   impl_->iteration_data.reset();
   impl_->transform.reset();
   mark_clean();
 }
 
-void barrier_cache_t::store_iteration_data(barrier_iteration_data_t* data)
+template <typename i_t, typename f_t>
+void barrier_cache_t<i_t, f_t>::store_iteration_data(barrier::iteration_data_t<i_t, f_t>* data)
 {
   impl_->iteration_data.reset(data);
 }
 
-barrier_iteration_data_t* barrier_cache_t::release_iteration_data()
+template <typename i_t, typename f_t>
+barrier::iteration_data_t<i_t, f_t>* barrier_cache_t<i_t, f_t>::release_iteration_data()
 {
   return impl_->iteration_data.release();
 }
 
-void barrier_cache_t::store_transform(std::unique_ptr<barrier_transform_t> transform)
+template <typename i_t, typename f_t>
+void barrier_cache_t<i_t, f_t>::store_transform(
+  std::unique_ptr<barrier_transform_t<i_t, f_t>> transform)
 {
   impl_->transform = std::move(transform);
 }
 
-barrier_transform_t* barrier_cache_t::transform() { return impl_->transform.get(); }
+template <typename i_t, typename f_t>
+barrier_transform_t<i_t, f_t>* barrier_cache_t<i_t, f_t>::transform()
+{
+  return impl_->transform.get();
+}
 
-barrier_transform_t const* barrier_cache_t::transform() const { return impl_->transform.get(); }
+template <typename i_t, typename f_t>
+barrier_transform_t<i_t, f_t> const* barrier_cache_t<i_t, f_t>::transform() const
+{
+  return impl_->transform.get();
+}
 
-bool barrier_cache_t::dirty() const
+template <typename i_t, typename f_t>
+bool barrier_cache_t<i_t, f_t>::dirty() const
 {
   return (impl_->linear_objective_dirty || impl_->rhs_dirty) && impl_->transform != nullptr &&
          impl_->iteration_data.get() != nullptr;
 }
 
-void barrier_cache_t::mark_clean()
+template <typename i_t, typename f_t>
+void barrier_cache_t<i_t, f_t>::mark_clean()
 {
   impl_->linear_objective_dirty = false;
   impl_->rhs_dirty              = false;
   impl_->rhs_infeasible         = false;
 }
 
-bool barrier_cache_t::rhs_infeasible() const { return impl_->rhs_infeasible; }
+template <typename i_t, typename f_t>
+bool barrier_cache_t<i_t, f_t>::rhs_infeasible() const
+{
+  return impl_->rhs_infeasible;
+}
 
-void barrier_cache_t::update_linear_objective(double const* c, int n)
+template <typename i_t, typename f_t>
+void barrier_cache_t<i_t, f_t>::update_linear_objective(f_t const* c, i_t n)
 {
   require_cache(impl_->transform.get(), impl_->iteration_data.get(), "update_linear_objective");
   // Cached Q and c are in minimization space.
-  std::vector<double> user_objective;
+  std::vector<f_t> user_objective;
   if (impl_->transform->maximize && c != nullptr && n > 0) {
     user_objective.assign(c, c + n);
-    for (double& value : user_objective) {
+    for (f_t& value : user_objective) {
       value = -value;
     }
     c = user_objective.data();
   }
-  std::vector<double> crushed;
+  std::vector<f_t> crushed;
   try {
     crushed = crush_user_linear_objective(*impl_->transform, c, n);
   } catch (std::invalid_argument const& e) {
@@ -154,13 +191,13 @@ void barrier_cache_t::update_linear_objective(double const* c, int n)
   auto const& linear_obj_shift = impl_->transform->linear_obj_shift;
   auto const& column_scales    = impl_->transform->column_scales;
   auto const& translated_lower = impl_->transform->presolve_info.removed_lower_bounds;
-  simplex::lp_problem_t<int, double>& barrier_lp = *impl_->transform->barrier_lp;
+  simplex::lp_problem_t<i_t, f_t>& barrier_lp = *impl_->transform->barrier_lp;
   if (!translated_lower.empty() && linear_obj_shift.size() == crushed.size() &&
       column_scales.size() == crushed.size() && barrier_lp.objective.size() == crushed.size()) {
-    double obj_constant_delta = 0.0;
+    f_t obj_constant_delta    = 0.0;
     std::size_t const n_lower = std::min(translated_lower.size(), crushed.size());
     for (std::size_t j = 0; j < n_lower; ++j) {
-      double const crushed_before = barrier_lp.objective[j] - linear_obj_shift[j];
+      f_t const crushed_before = barrier_lp.objective[j] - linear_obj_shift[j];
       obj_constant_delta += (crushed[j] - crushed_before) * column_scales[j] * translated_lower[j];
     }
     barrier_lp.obj_constant += obj_constant_delta;
@@ -168,21 +205,22 @@ void barrier_cache_t::update_linear_objective(double const* c, int n)
   add_shift(crushed, linear_obj_shift);
   // The next solve builds its solver from barrier_lp, so keep its objective and the cached
   // iteration workspace on the same c.
-  std::vector<double>& barrier_objective = barrier_lp.objective;
+  std::vector<f_t>& barrier_objective = barrier_lp.objective;
   cuopt_expects(barrier_objective.size() == crushed.size(),
                 error_type_t::ValidationError,
                 "update_linear_objective: crushed objective size does not match the cached "
                 "barrier LP.");
   barrier_objective = crushed;
   barrier::apply_barrier_linear_objective(
-    *impl_->iteration_data, crushed.data(), static_cast<int>(crushed.size()));
+    *impl_->iteration_data, crushed.data(), static_cast<i_t>(crushed.size()));
   impl_->linear_objective_dirty = true;
 }
 
-void barrier_cache_t::update_rhs(double const* b, int m)
+template <typename i_t, typename f_t>
+void barrier_cache_t<i_t, f_t>::update_rhs(f_t const* b, i_t m)
 {
   require_cache(impl_->transform.get(), impl_->iteration_data.get(), "update_rhs");
-  std::vector<double> crushed;
+  std::vector<f_t> crushed;
   std::string error;
   crush_rhs_status_t const status = crush_user_rhs(*impl_->transform, b, m, crushed, error);
   if (status == crush_rhs_status_t::infeasible) {
@@ -198,14 +236,18 @@ void barrier_cache_t::update_rhs(double const* b, int m)
   add_shift(crushed, impl_->transform->rhs_shift);
   // barrier_lp->rhs also seeds the next solve's Mehrotra start, so keep it and the cached
   // workspace on the same b.
-  std::vector<double>& barrier_rhs = impl_->transform->barrier_lp->rhs;
+  std::vector<f_t>& barrier_rhs = impl_->transform->barrier_lp->rhs;
   cuopt_expects(barrier_rhs.size() == crushed.size(),
                 error_type_t::ValidationError,
                 "update_rhs: crushed RHS size does not match the cached barrier LP.");
   barrier_rhs = crushed;
   barrier::apply_barrier_rhs(
-    *impl_->iteration_data, crushed.data(), static_cast<int>(crushed.size()));
+    *impl_->iteration_data, crushed.data(), static_cast<i_t>(crushed.size()));
   impl_->rhs_dirty = true;
 }
 
+#ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE
+template class barrier_cache_t<int, double>;
+#endif
+
 }  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/barrier/barrier_transform.hpp b/cpp/src/barrier/barrier_transform.hpp
index 8212a6ca3d..c190b74ace 100644
--- a/cpp/src/barrier/barrier_transform.hpp
+++ b/cpp/src/barrier/barrier_transform.hpp
@@ -8,73 +8,256 @@
 #pragma once
 
 #include <dual_simplex/presolve.hpp>
+#include <dual_simplex/user_problem.hpp>
 #include <linear_algebra/sparse_matrix.hpp>
 
+#include <algorithm>
+#include <cmath>
+#include <limits>
 #include <memory>
 #include <stdexcept>
 #include <string>
+#include <utility>
 #include <vector>
 
 namespace cuopt::mathematical_optimization {
 
 /**
- * User space to presolve space transform retained on barrier_cache_t after Optimal:
+ * User space to presolved space transform retained on barrier_cache_t after Optimal:
  * convert / presolve / scaling, plus the scaled LP.
  * Enough to crush new linear objective or RHS data from the original problem into
  * presolved space and to uncrush a solution without rerunning those algorithms.
  */
+template <typename i_t, typename f_t>
 struct barrier_transform_t {
-  int user_num_cols{0};
-  int user_num_rows{0};
-  int original_num_cols{0};
-  int original_num_rows{0};
-  double obj_scale{1.0};
-  double obj_constant{0.0};
+  i_t user_num_cols{0};
+  i_t user_num_rows{0};
+  i_t original_num_cols{0};
+  i_t original_num_rows{0};
+  f_t obj_scale{1.0};
+  f_t obj_constant{0.0};
   bool maximize{false};  // first-solve sense; cached Q/c were negated if true
 
   // Enough of the user problem for reuse uncrush without rebuilding A.
   std::vector<char> row_sense;
-  int cone_var_start{0};
-  std::vector<int> second_order_cone_dims;
-  int expanded_original_num_cols{0};
-  std::vector<int> original_col_to_expanded_col;
-
-  cuopt::mathematical_optimization::simplex::presolve_info_t<int, double> presolve_info;
-  std::vector<double> column_scales;
-  std::vector<double> row_scales;
+  i_t cone_var_start{0};
+  // cone_var_start after convert, which inserts inequality slacks ahead of the cone block.
+  i_t converted_cone_var_start{0};
+  std::vector<i_t> second_order_cone_dims;
+  // Quadratic constraint count the cached expansion was built from.
+  i_t num_quadratic_constraints{0};
+  // Dimensions before the QCMATRIX->SOC expansion, which permutes columns and appends rows.
+  // Updates arrive in these coordinates, not the expanded ones. Zero when no expansion ran.
+  i_t pre_expansion_num_cols{0};
+  i_t pre_expansion_num_rows{0};
+  std::vector<i_t> original_col_to_expanded_col;
+  // RHS of the rows the expansion appended. Fixed by the quadratic constraints, so an RHS
+  // update keeps them and only overwrites the problem's own rows.
+  std::vector<f_t> cone_row_rhs;
+  std::vector<simplex::cone_head_bound_t<i_t, f_t>> cone_head_bounds;
+
+  cuopt::mathematical_optimization::simplex::presolve_info_t<i_t, f_t> presolve_info;
+  std::vector<f_t> column_scales;
+  std::vector<f_t> row_scales;
   // Barrier linear objective minus crush(user c) from the first solve (Q*ell shift, etc.).
-  std::vector<double> linear_obj_shift;
+  std::vector<f_t> linear_obj_shift;
   // Barrier RHS minus crush(user b) from the first solve (fixed/lower-bound shifts).
-  std::vector<double> rhs_shift;
+  std::vector<f_t> rhs_shift;
   // False when range rows or folding put the user RHS somewhere other than barrier_lp->rhs.
   bool rhs_update_supported{false};
-  std::unique_ptr<cuopt::mathematical_optimization::simplex::lp_problem_t<int, double>> barrier_lp;
+  std::unique_ptr<cuopt::mathematical_optimization::simplex::lp_problem_t<i_t, f_t>> barrier_lp;
   // CSC Q with slack columns, as consumed by iteration_data_t. Not the same object as
   // barrier_lp->Q.
-  std::unique_ptr<csc_matrix_t<int, double>> barrier_Q;
+  std::unique_ptr<csc_matrix_t<i_t, f_t>> barrier_Q;
 };
 
-inline bool can_reuse_barrier_cache(barrier_transform_t const* xf,
-                                    int bound_free_variables,
-                                    int num_cols,
-                                    int num_rows,
+// Shared reuse gate from the update-API work. A cone problem skips the size
+// check: the caller compares either expanded user counts or pre-expansion
+// problem counts, which are not the same number.
+template <typename i_t, typename f_t>
+inline bool can_reuse_barrier_cache(barrier_transform_t<i_t, f_t> const* xf,
+                                    i_t bound_free_variables,
+                                    i_t num_cols,
+                                    i_t num_rows,
                                     bool has_quadratic_objective,
                                     bool user_has_soc)
 {
-  return xf != nullptr && xf->barrier_lp != nullptr && has_quadratic_objective && !user_has_soc &&
-         xf->second_order_cone_dims.empty() && xf->barrier_lp->second_order_cone_dims.empty() &&
-         bound_free_variables == 0 && xf->presolve_info.bounded_free_variables.empty() &&
-         static_cast<int>(xf->row_sense.size()) == xf->user_num_rows &&
-         num_cols == xf->user_num_cols && num_rows == xf->user_num_rows;
+  if (xf == nullptr || xf->barrier_lp == nullptr) { return false; }
+  if (bound_free_variables != 0 || !xf->presolve_info.bounded_free_variables.empty()) {
+    return false;
+  }
+  if (user_has_soc) { return true; }
+  // LP and QP reuse when the sizes match. A cached quadratic objective does not match an LP.
+  return xf->second_order_cone_dims.empty() && xf->barrier_lp->second_order_cone_dims.empty() &&
+         static_cast<i_t>(xf->row_sense.size()) == xf->user_num_rows &&
+         num_cols == xf->user_num_cols && num_rows == xf->user_num_rows &&
+         (has_quadratic_objective || xf->barrier_lp->Q.n == 0);
+}
+
+// Column count an update is sized in: the pre-expansion count when the SOC expansion grew the
+// problem, otherwise the cached user count.
+template <typename i_t, typename f_t>
+inline i_t problem_num_cols(barrier_transform_t<i_t, f_t> const& xf)
+{
+  return xf.pre_expansion_num_cols > 0 ? xf.pre_expansion_num_cols : xf.user_num_cols;
+}
+
+// Row count an update is sized in.
+template <typename i_t, typename f_t>
+inline i_t problem_num_rows(barrier_transform_t<i_t, f_t> const& xf)
+{
+  return xf.original_col_to_expanded_col.empty() ? xf.user_num_rows : xf.pre_expansion_num_rows;
+}
+
+// The expansion permutes problem columns into a [linear | cone] layout. An empty map means no
+// expansion ran, so the layouts coincide.
+template <typename i_t, typename f_t>
+inline i_t problem_col_to_expanded_col(barrier_transform_t<i_t, f_t> const& xf, i_t problem_col)
+{
+  return xf.original_col_to_expanded_col.empty() ? problem_col
+                                                 : xf.original_col_to_expanded_col[problem_col];
+}
+
+// convert inserts the inequality slacks ahead of the cone block, pushing every cone column
+// right. Mirrors user_col_to_problem_col in presolve.cpp.
+template <typename i_t, typename f_t>
+inline i_t expanded_col_to_converted_col(barrier_transform_t<i_t, f_t> const& xf, i_t expanded_col)
+{
+  if (xf.second_order_cone_dims.empty() || xf.converted_cone_var_start <= xf.cone_var_start ||
+      expanded_col < xf.cone_var_start) {
+    return expanded_col;
+  }
+  return xf.converted_cone_var_start + (expanded_col - xf.cone_var_start);
+}
+
+// Move the problem's own coefficients into the expanded layout the cached problem is sized for.
+template <typename i_t, typename f_t, typename value_t>
+std::vector<value_t> scatter_problem_objective(barrier_transform_t<i_t, f_t> const& xf,
+                                               std::vector<value_t> const& problem_objective)
+{
+  std::vector<value_t> expanded(static_cast<std::size_t>(xf.user_num_cols), value_t(0));
+  for (i_t j = 0; j < static_cast<i_t>(problem_objective.size()); ++j) {
+    expanded[problem_col_to_expanded_col(xf, j)] = problem_objective[j];
+  }
+  return expanded;
+}
+
+// Inverse of scatter_problem_objective, giving crush_user_linear_objective the problem-sized input
+// it expects. The expansion adds no objective coefficients, so nothing is lost.
+template <typename i_t, typename f_t>
+std::vector<f_t> gather_problem_objective(barrier_transform_t<i_t, f_t> const& xf,
+                                          std::vector<f_t> const& expanded_objective)
+{
+  std::vector<f_t> problem_objective(static_cast<std::size_t>(problem_num_cols(xf)));
+  for (i_t j = 0; j < static_cast<i_t>(problem_objective.size()); ++j) {
+    problem_objective[j] = expanded_objective[problem_col_to_expanded_col(xf, j)];
+  }
+  return problem_objective;
+}
+
+// Reuse never re-runs the expansion, so the cone block must be laid out exactly as the cached
+// one left it.
+template <typename i_t, typename f_t>
+bool cone_layout_matches(barrier_transform_t<i_t, f_t> const& xf,
+                         simplex::user_problem_t<i_t, f_t> const& user_problem)
+{
+  return user_problem.cone_var_start == xf.cone_var_start &&
+         user_problem.second_order_cone_dims.size() == xf.second_order_cone_dims.size() &&
+         std::equal(user_problem.second_order_cone_dims.begin(),
+                    user_problem.second_order_cone_dims.end(),
+                    xf.second_order_cone_dims.begin(),
+                    [](i_t dim, i_t cached) { return dim == cached; });
+}
+
+// Equality substitution is the last presolve step and is not represented by remaining_variables
+// or free_variable_pairs. It drops zero-cost free columns (and their pivot rows). Indices in the
+// elimination records are into the vector produced by the earlier presolve maps.
+template <typename i_t, typename f_t>
+inline void drop_substituted_free_variables(barrier_transform_t<i_t, f_t> const& xf,
+                                            std::vector<f_t>& values,
+                                            char const* what)
+{
+  auto const& eliminations = xf.presolve_info.free_variable_eliminations;
+  if (eliminations.empty()) { return; }
+  auto const& keep = xf.presolve_info.free_elimination_remaining_variables;
+  for (auto const& elimination : eliminations) {
+    if (elimination.variable < 0 ||
+        static_cast<std::size_t>(elimination.variable) >= values.size()) {
+      throw std::invalid_argument(std::string(what) +
+                                  ": eliminated free variable index is out of range.");
+    }
+    // Substitution is valid only while this coefficient stays zero. A nonzero cost would have
+    // kept the column, so the cached factorization cannot be reused.
+    if (values[elimination.variable] != 0.0) {
+      throw std::invalid_argument(
+        std::string(what) +
+        ": a free variable removed by equality substitution has a nonzero objective "
+        "coefficient; run a full Solve.");
+    }
+  }
+  std::vector<f_t> reduced(keep.size());
+  for (std::size_t k = 0; k < keep.size(); ++k) {
+    i_t const column = keep[k];
+    if (column < 0 || static_cast<std::size_t>(column) >= values.size()) {
+      throw std::invalid_argument(std::string(what) +
+                                  ": free-elimination column index is out of range.");
+    }
+    reduced[k] = values[column];
+  }
+  values = std::move(reduced);
 }
 
-inline std::vector<double> crush_user_linear_objective(barrier_transform_t const& xf,
-                                                       double const* c,
-                                                       int n)
+// Replay rhs[row] -= factor * rhs[pivot] in elimination order, then drop the pivot rows.
+// Factors come from the matrix and do not depend on the right-hand side.
+template <typename i_t, typename f_t>
+inline void substitute_free_variable_rows(barrier_transform_t<i_t, f_t> const& xf,
+                                          std::vector<f_t>& values,
+                                          char const* what)
 {
-  if (c == nullptr || n != xf.user_num_cols) {
+  auto const& eliminations = xf.presolve_info.free_variable_eliminations;
+  if (eliminations.empty()) { return; }
+  for (auto const& elimination : eliminations) {
+    if (elimination.pivot_row < 0 ||
+        static_cast<std::size_t>(elimination.pivot_row) >= values.size()) {
+      throw std::invalid_argument(std::string(what) +
+                                  ": free-elimination pivot row is out of range.");
+    }
+    if (elimination.affected_rows.size() != elimination.factors.size()) {
+      throw std::invalid_argument(std::string(what) +
+                                  ": free-elimination factors do not match affected rows.");
+    }
+    f_t const pivot_rhs = values[elimination.pivot_row];
+    for (std::size_t k = 0; k < elimination.affected_rows.size(); ++k) {
+      i_t const row = elimination.affected_rows[k];
+      if (row < 0 || static_cast<std::size_t>(row) >= values.size()) {
+        throw std::invalid_argument(std::string(what) +
+                                    ": free-elimination affected row is out of range.");
+      }
+      values[row] -= elimination.factors[k] * pivot_rhs;
+    }
+  }
+  auto const& keep = xf.presolve_info.free_elimination_remaining_constraints;
+  std::vector<f_t> reduced(keep.size());
+  for (std::size_t k = 0; k < keep.size(); ++k) {
+    i_t const row = keep[k];
+    if (row < 0 || static_cast<std::size_t>(row) >= values.size()) {
+      throw std::invalid_argument(std::string(what) +
+                                  ": free-elimination row index is out of range.");
+    }
+    reduced[k] = values[row];
+  }
+  values = std::move(reduced);
+}
+
+template <typename i_t, typename f_t>
+inline std::vector<f_t> crush_user_linear_objective(barrier_transform_t<i_t, f_t> const& xf,
+                                                    f_t const* c,
+                                                    i_t n)
+{
+  if (c == nullptr || n != problem_num_cols(xf)) {
     throw std::invalid_argument(
-      "update_linear_objective: linear objective length must match the cached user column count.");
+      "update_linear_objective: linear objective length must match the cached problem column "
+      "count.");
   }
   if (xf.original_num_cols < xf.user_num_cols) {
     throw std::invalid_argument(
@@ -83,20 +266,32 @@ inline std::vector<double> crush_user_linear_objective(barrier_transform_t const
   if (xf.barrier_lp == nullptr) {
     throw std::invalid_argument("update_linear_objective: cached barrier LP is missing.");
   }
+  if (!xf.original_col_to_expanded_col.empty() &&
+      static_cast<i_t>(xf.original_col_to_expanded_col.size()) != n) {
+    throw std::invalid_argument(
+      "update_linear_objective: cached column map does not cover the problem columns.");
+  }
 
-  std::vector<double> orig(static_cast<std::size_t>(xf.original_num_cols), 0.0);
-  for (int j = 0; j < n; ++j) {
-    orig[static_cast<std::size_t>(j)] = c[j];
+  // The expansion leaves the variables it adds out of the objective, so only the positions of
+  // the problem's own coefficients move.
+  std::vector<f_t> orig(static_cast<std::size_t>(xf.original_num_cols), 0.0);
+  for (i_t j = 0; j < n; ++j) {
+    i_t const converted_col = expanded_col_to_converted_col(xf, problem_col_to_expanded_col(xf, j));
+    if (converted_col < 0 || converted_col >= xf.original_num_cols) {
+      throw std::invalid_argument(
+        "update_linear_objective: cached column map points outside the converted problem.");
+    }
+    orig[converted_col] = c[j];
   }
-  for (int j : xf.presolve_info.negated_variables) {
-    orig[static_cast<std::size_t>(j)] *= -1.0;
+  for (i_t j : xf.presolve_info.negated_variables) {
+    orig[j] *= -1.0;
   }
 
-  std::vector<double> presolved;
+  std::vector<f_t> presolved;
   if (!xf.presolve_info.remaining_variables.empty()) {
     presolved.resize(xf.presolve_info.remaining_variables.size());
     for (std::size_t k = 0; k < xf.presolve_info.remaining_variables.size(); ++k) {
-      presolved[k] = orig[static_cast<std::size_t>(xf.presolve_info.remaining_variables[k])];
+      presolved[k] = orig[xf.presolve_info.remaining_variables[k]];
     }
   } else {
     presolved = std::move(orig);
@@ -110,13 +305,15 @@ inline std::vector<double> crush_user_linear_objective(barrier_transform_t const
     std::size_t extra = pairs.size() / 2;
     presolved.resize(presolved.size() + extra);
     for (std::size_t k = 0; k < extra; ++k) {
-      int u                                  = pairs[2 * k];
-      int v                                  = pairs[2 * k + 1];
-      presolved[static_cast<std::size_t>(v)] = -presolved[static_cast<std::size_t>(u)];
+      i_t u        = pairs[2 * k];
+      i_t v        = pairs[2 * k + 1];
+      presolved[v] = -presolved[u];
     }
   }
 
-  if (static_cast<int>(presolved.size()) != xf.barrier_lp->num_cols ||
+  drop_substituted_free_variables<i_t>(xf, presolved, "update_linear_objective");
+
+  if (static_cast<i_t>(presolved.size()) != xf.barrier_lp->num_cols ||
       xf.column_scales.size() != presolved.size()) {
     throw std::invalid_argument(
       "update_linear_objective: crushed objective size does not match barrier columns / "
@@ -134,15 +331,18 @@ inline std::vector<double> crush_user_linear_objective(barrier_transform_t const
 enum class crush_rhs_status_t { success = 0, invalid = -1, infeasible = -2 };
 
 template <typename i_t, typename f_t>
-inline crush_rhs_status_t crush_user_rhs(
-  barrier_transform_t const& xf, f_t const* b, i_t m, std::vector<f_t>& crushed, std::string& error)
+inline crush_rhs_status_t crush_user_rhs(barrier_transform_t<i_t, f_t> const& xf,
+                                         f_t const* b,
+                                         i_t m,
+                                         std::vector<f_t>& crushed,
+                                         std::string& error)
 {
   auto invalid = [&](char const* message) {
     error = message;
     return crush_rhs_status_t::invalid;
   };
-  if (b == nullptr || m != xf.user_num_rows) {
-    return invalid("update_rhs: RHS length must match the cached user row count.");
+  if (m != problem_num_rows(xf) || (b == nullptr && m != 0)) {
+    return invalid("update_rhs: RHS length must match the cached problem row count.");
   }
   if (!xf.rhs_update_supported) {
     return invalid("update_rhs: cached convert used range rows or folding; run a full Solve.");
@@ -155,19 +355,46 @@ inline crush_rhs_status_t crush_user_rhs(
   }
   if (xf.barrier_lp == nullptr) { return invalid("update_rhs: cached barrier LP is missing."); }
 
+  // The quadratic constraints fix the RHS of the appended rows, so an update overwrites the
+  // problem's own rows and keeps the cached tail. The tail is empty without an expansion.
+  std::vector<f_t> expanded(m);
+  for (i_t i = 0; i < m; ++i) {
+    expanded[i] = b[i];
+  }
+  for (f_t rhs : xf.cone_row_rhs) {
+    expanded.push_back(rhs);
+  }
+  if (static_cast<i_t>(expanded.size()) != xf.user_num_rows) {
+    return invalid("update_rhs: cached cone-row RHS does not span the expanded rows.");
+  }
+
+  // The expansion proved these heads nonnegative from the old RHS. A full solve rejects the
+  // problem once that no longer holds, so re-prove it here rather than trust the cached verdict.
+  for (simplex::cone_head_bound_t<i_t, f_t> const& bound : xf.cone_head_bounds) {
+    f_t implied = -std::numeric_limits<f_t>::infinity();
+    for (auto const& [row, coefficient] : bound.rows) {
+      implied = std::max(implied, expanded[row] / coefficient);
+    }
+    if (!(implied >= 0.0)) {
+      error = "update_rhs: new RHS no longer implies second-order cone head variable " +
+              std::to_string(bound.head_col) + " is nonnegative.";
+      return crush_rhs_status_t::invalid;
+    }
+  }
+
   // convert turns 'G' rows into 'L' rows by negating the row and its RHS.
   std::vector<f_t> original(static_cast<std::size_t>(xf.original_num_rows));
-  for (i_t i = 0; i < m; ++i) {
-    original[static_cast<std::size_t>(i)] =
-      xf.row_sense[static_cast<std::size_t>(i)] == 'G' ? -b[i] : b[i];
+  const i_t user_num_rows = xf.user_num_rows;
+  for (i_t i = 0; i < user_num_rows; ++i) {
+    original[i] = xf.row_sense[i] == 'G' ? -expanded[i] : expanded[i];
   }
 
   // Presolve drops only empty equalities, and only when the RHS is exactly 0.
   for (i_t i : xf.presolve_info.removed_constraints) {
-    if (i < 0 || i >= m) {
+    if (i < 0 || i >= xf.user_num_rows) {
       return invalid("update_rhs: removed constraint index is out of range.");
     }
-    if (original[i] != f_t(0)) { return crush_rhs_status_t::infeasible; }
+    if (original[i] != 0.0) { return crush_rhs_status_t::infeasible; }
   }
 
   // Empty remaining_constraints means either no empty-row pass ran, or every row was dropped
@@ -176,12 +403,14 @@ inline crush_rhs_status_t crush_user_rhs(
   if (!xf.presolve_info.remaining_constraints.empty()) {
     presolved.resize(xf.presolve_info.remaining_constraints.size());
     for (std::size_t k = 0; k < xf.presolve_info.remaining_constraints.size(); ++k) {
-      presolved[k] = original[static_cast<std::size_t>(xf.presolve_info.remaining_constraints[k])];
+      presolved[k] = original[xf.presolve_info.remaining_constraints[k]];
     }
   } else if (xf.presolve_info.removed_constraints.empty()) {
     presolved = std::move(original);
   }
 
+  substitute_free_variable_rows<i_t>(xf, presolved, "update_rhs");
+
   if (static_cast<i_t>(presolved.size()) != xf.barrier_lp->num_rows ||
       xf.row_scales.size() != presolved.size()) {
     return invalid("update_rhs: crushed RHS size does not match barrier rows / row_scales.");
@@ -189,7 +418,7 @@ inline crush_rhs_status_t crush_user_rhs(
   for (std::size_t i = 0; i < presolved.size(); ++i) {
     presolved[i] /= xf.row_scales[i];
   }
-  crushed = std::move(presolved);
+  crushed.assign(presolved.begin(), presolved.end());
   return crush_rhs_status_t::success;
 }
 
diff --git a/cpp/src/barrier/second_order_cone_kernels.cuh b/cpp/src/barrier/second_order_cone_kernels.cuh
index 06d5b2cbe1..c5a26853f1 100644
--- a/cpp/src/barrier/second_order_cone_kernels.cuh
+++ b/cpp/src/barrier/second_order_cone_kernels.cuh
@@ -19,6 +19,7 @@
 #include <raft/core/device_span.hpp>
 
 #include <thrust/binary_search.h>
+#include <thrust/fill.h>
 #include <thrust/gather.h>
 #include <thrust/iterator/counting_iterator.h>
 #include <thrust/iterator/permutation_iterator.h>
@@ -34,7 +35,6 @@
 #include <cub/device/device_reduce.cuh>
 #include <cub/warp/warp_reduce.cuh>
 
-#include <concepts>
 #include <cstddef>
 #include <numeric>
 #include <span>
@@ -65,7 +65,7 @@ inline constexpr int soc_block_size = 256;
 /**
  * Tail aggregates for the cone step-length reduction (CUB value type).
  */
-template <std::floating_point f_t>
+template <typename f_t>
 struct step_tail_sums_t {
   f_t du_tail_sq{};
   f_t u_tail_du_tail{};
@@ -84,7 +84,7 @@ struct step_tail_sums_t {
  * slots and `temp_cone` sequentially inside a higher-level operation, but no
  * persistent NT scaling or iterate state is stored here.
  */
-template <std::integral i_t, std::floating_point f_t, int n_slots = 3>
+template <typename i_t, typename f_t, int n_slots = 3>
 struct cone_scratch_t {
   i_t n_cones;            // number of SOC blocks
   size_t n_cone_entries;  // total packed cone dimension
@@ -156,7 +156,7 @@ struct to_size_t_t {
  * solver's global x/z vectors. The caller must keep the underlying storage
  * alive for the lifetime of this object.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 struct cone_data_t {
   // Topology. This is immutable after construction.
   i_t n_cones;            // number of SOC blocks
@@ -293,7 +293,7 @@ struct cone_data_t {
   i_t expansion_var_count() const { return 2 * n_sparse_cones; }
 };
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   nt_finalize_scaling_scalars_kernel(raft::device_span<const f_t> x,
                                      raft::device_span<const f_t> z,
@@ -317,7 +317,7 @@ __global__ void __launch_bounds__(soc_block_size)
   eta[cone]     = sqrt(z_scale[cone] / x_scale[cone]);
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   nt_finalize_w_scale_kernel(raft::device_span<const f_t> w,
                              raft::device_span<const f_t> tail_sq,
@@ -341,7 +341,7 @@ __global__ void __launch_bounds__(soc_block_size)
  *   w_0 = z_0 / z_scale + x_0 / x_scale
  *   w_tail = z_tail / z_scale - x_tail / x_scale.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   nt_write_w_kernel(raft::device_span<const f_t> x,
                     raft::device_span<const f_t> z,
@@ -364,7 +364,7 @@ __global__ void __launch_bounds__(soc_block_size)
   w[idx] = z[idx] / z_scale[cone] - x[idx] / x_scale[cone];
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   nt_normalize_w_kernel(raft::device_span<f_t> w,
                         raft::device_span<const f_t> w_scale,
@@ -377,7 +377,7 @@ __global__ void __launch_bounds__(soc_block_size)
   w[idx] /= w_scale[cone];
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   nt_finalize_head_kernel(raft::device_span<f_t> w,
                           raft::device_span<const f_t> normalized_tail_sq,
@@ -390,7 +390,7 @@ __global__ void __launch_bounds__(soc_block_size)
   w[cone_offsets[cone]] = sqrt(1 + normalized_tail_sq[cone]);
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   nt_write_lambda_kernel(raft::device_span<const f_t> x,
                          raft::device_span<const f_t> z,
@@ -420,7 +420,7 @@ __global__ void __launch_bounds__(soc_block_size)
 
   const f_t x_head  = x[cone_off];
   const f_t z_head  = z[cone_off];
-  const f_t denom   = z_head / z_scale_cone + x_head / x_scale_cone + static_cast<f_t>(2) * gamma;
+  const f_t denom   = z_head / z_scale_cone + x_head / x_scale_cone + 2.0 * gamma;
   const f_t coeff_z = (gamma + x_head / x_scale_cone) / z_scale_cone;
   const f_t coeff_x = (gamma + z_head / z_scale_cone) / x_scale_cone;
 
@@ -447,7 +447,7 @@ __global__ void __launch_bounds__(soc_block_size)
  *   0: ||x_tail||^2 -> x_scale
  *   1: ||z_tail||^2 -> z_scale
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void launch_nt_scaling(cone_data_t<i_t, f_t>& cones, cuda::stream_ref stream)
 {
   auto x_scale = cones.scratch.template get_slot<0>();
@@ -533,7 +533,7 @@ void launch_nt_scaling(cone_data_t<i_t, f_t>& cones, cuda::stream_ref stream)
 // One block per sparse cone. Recompute the rank-2 factors (corner d and the
 // vectors v, u, both scaled by eta^2) from the current NT direction w so that
 // the implicit block reproduces the dense H = eta^2 (2 w w^T - J).
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void update_scaling_sparse_kernel(raft::device_span<const f_t> w,
                                              raft::device_span<const f_t> eta,
                                              raft::device_span<f_t> d,
@@ -556,18 +556,18 @@ __global__ void update_scaling_sparse_kernel(raft::device_span<const f_t> w,
   const i_t block_start = sparse_entry_offsets[sparse_idx];
 
   if (threadIdx.x == 0) {
-    const f_t alpha    = f_t(2) * w[cone_off];
-    const f_t wsq      = f_t(2) * w[cone_off] * w[cone_off] - f_t(1);
-    const f_t wsq_safe = f_t(0.5) * (wsq + sqrt(wsq * wsq + f_t(1)));
-    const f_t wsqinv   = f_t(1) / wsq_safe;
+    const f_t alpha    = 2.0 * w[cone_off];
+    const f_t wsq      = 2.0 * w[cone_off] * w[cone_off] - 1.0;
+    const f_t wsq_safe = f_t(0.5) * (wsq + sqrt(wsq * wsq + 1.0));
+    const f_t wsqinv   = 1.0 / wsq_safe;
     const f_t di       = f_t(0.5) * wsqinv;
     d[sparse_idx]      = di;
     const f_t radicand = wsq_safe - di;
     const f_t u0       = sqrt(max(radicand, f_t(0)));
-    const f_t u1       = (u0 > f_t(0)) ? alpha / u0 : f_t(0);
-    const f_t v0       = f_t(0);
-    const f_t denom    = f_t(2) * wsq_safe - wsqinv;
-    const f_t v1_arg   = (abs(denom) > f_t(1e-12)) ? f_t(2) * (f_t(2) + wsqinv) / denom : f_t(2);
+    const f_t u1       = (u0 > 0.0) ? alpha / u0 : 0.0;
+    const f_t v0       = 0.0;
+    const f_t denom    = 2.0 * wsq_safe - wsqinv;
+    const f_t v1_arg   = (abs(denom) > f_t(1e-12)) ? 2.0 * (2.0 + wsqinv) / denom : 2.0;
     const f_t v1       = sqrt(max(v1_arg, f_t(0)));
     const f_t eta_sq   = eta[cone_idx] * eta[cone_idx];
     s_mem[0]           = eta_sq * u0;
@@ -598,7 +598,7 @@ __global__ void update_scaling_sparse_kernel(raft::device_span<const f_t> w,
  * scaling, so the implicit sparse block matches the dense Hessian for this
  * iteration. Call after `launch_nt_scaling` has updated w and eta.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void launch_update_scaling_sparse(cone_data_t<i_t, f_t>& cones, cuda::stream_ref stream)
 {
   if (!cones.has_sparse_cones()) { return; }
@@ -618,7 +618,7 @@ void launch_update_scaling_sparse(cone_data_t<i_t, f_t>& cones, cuda::stream_ref
   RAFT_CUDA_TRY(cudaPeekAtLastError());
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   apply_w_inv_write_kernel(raft::device_span<const f_t> v,
                            raft::device_span<f_t> out,
@@ -638,18 +638,18 @@ __global__ void __launch_bounds__(soc_block_size)
   const f_t w0      = w[cone_off];
   const f_t zeta    = tail_dot[cone];
   const f_t v0      = v[cone_off];
-  const f_t inv_eta = f_t(1) / eta[cone];
+  const f_t inv_eta = 1.0 / eta[cone];
 
   if (local_idx == 0) {
     out[idx] = inv_eta * (w0 * v0 - zeta);
     return;
   }
 
-  const f_t coeff = -v0 + zeta / (f_t(1) + w0);
+  const f_t coeff = -v0 + zeta / (1.0 + w0);
   out[idx]        = inv_eta * (v[idx] + coeff * w[idx]);
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   apply_w_write_kernel(raft::device_span<const f_t> v,
                        raft::device_span<f_t> out,
@@ -676,11 +676,11 @@ __global__ void __launch_bounds__(soc_block_size)
     return;
   }
 
-  const f_t coeff = v0 + zeta / (f_t(1) + w0);
+  const f_t coeff = v0 + zeta / (1.0 + w0);
   out[idx]        = cone_eta * (v[idx] + coeff * w[idx]);
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   apply_hessian_kernel(raft::device_span<const f_t> v,
                        raft::device_span<f_t> out,
@@ -712,7 +712,7 @@ __global__ void __launch_bounds__(soc_block_size)
   out[idx] = bias.empty() ? h_value : bias_scale * bias[idx] + h_value;
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   gather_cone_heads_kernel(raft::device_span<const f_t> values,
                            raft::device_span<f_t> heads,
@@ -737,7 +737,7 @@ __global__ void __launch_bounds__(soc_block_size)
  *   d_0    = <scaled_dx, scaled_dz> - sigma_mu
  *   d_tail = scaled_dx_0 * scaled_dz_tail + scaled_dz_0 * scaled_dx_tail.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   combined_cone_shift_write_kernel(raft::device_span<f_t> shift,
                                    raft::device_span<const f_t> scaled_dx,
@@ -773,7 +773,7 @@ __global__ void __launch_bounds__(soc_block_size)
  * A second flat kernel writes `-p`, which lets the final W^{-1} call produce
  * q = -W^{-1} p without adding an output-scale argument to W^{-1}.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   jordan_divide_by_lambda_scalar_kernel(raft::device_span<const f_t> shift,
                                         raft::device_span<const f_t> nt_point,
@@ -797,7 +797,7 @@ __global__ void __launch_bounds__(soc_block_size)
   inv_lambda0[cone] = 1 / lambda0;
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   jordan_divide_by_lambda_write_kernel(raft::device_span<const f_t> shift,
                                        raft::device_span<const f_t> nt_point,
@@ -830,7 +830,7 @@ __global__ void __launch_bounds__(soc_block_size)
  *   (W^{-1}v)_0 = inv_eta * (w_0 v_0 - zeta)
  *   (W^{-1}v)_tail = inv_eta * (v_tail + (-v_0 + zeta / (1 + w_0)) w_tail)
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void apply_w_inv(raft::device_span<const f_t> v,
                  raft::device_span<f_t> out,
                  cone_data_t<i_t, f_t>& cones,
@@ -866,7 +866,7 @@ void apply_w_inv(raft::device_span<const f_t> v,
  *   (W * v)_tail =
  *     eta * (v_tail + (v_0 + zeta / (1 + w_0)) w_tail)
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void apply_w(raft::device_span<const f_t> v,
              raft::device_span<f_t> out,
              cone_data_t<i_t, f_t>& cones,
@@ -899,7 +899,7 @@ void apply_w(raft::device_span<const f_t> v,
  *   (Hv)_0 = eta^{2} (2 w_0 rho - v_0)
  *   (Hv)_tail = eta^{2} (2 w_tail rho + v_tail)
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void apply_hessian(raft::device_span<const f_t> v,
                    raft::device_span<f_t> out,
                    cone_data_t<i_t, f_t>& cones,
@@ -944,7 +944,7 @@ void apply_hessian(raft::device_span<const f_t> v,
  * This function applies the cone block H = S^2 and writes:
  *   dz = cone_target - H dx.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void recover_cone_dz_from_target(raft::device_span<const f_t> dx,
                                  cone_data_t<i_t, f_t>& cones,
                                  raft::device_span<const f_t> cone_target,
@@ -958,7 +958,7 @@ void recover_cone_dz_from_target(raft::device_span<const f_t> dx,
  * Accumulate the dense SOC cone-block matvec into an existing output vector:
  *   out += H x, where H = S^2, applied to dense cones only.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void launch_dense_hessian_matvec(raft::device_span<const f_t> x,
                                  cone_data_t<i_t, f_t>& cones,
                                  raft::device_span<f_t> out,
@@ -970,7 +970,7 @@ void launch_dense_hessian_matvec(raft::device_span<const f_t> x,
 
 // Bucket index owning `entry` via upper-bound search on the exclusive prefix ends in
 // offsets[1..n_buckets].
-template <std::integral i_t, typename offset_t>
+template <typename i_t, typename offset_t>
 __device__ i_t bucket_index(raft::device_span<const offset_t> offsets,
                             offset_t entry,
                             i_t n_buckets)
@@ -988,7 +988,7 @@ __device__ i_t bucket_index(raft::device_span<const offset_t> offsets,
   return lo;
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size) scatter_sparse_hessian_into_augmented_kernel(
   raft::device_span<f_t> augmented_x,
   raft::device_span<f_t> Hs_diag,
@@ -1043,7 +1043,7 @@ __global__ void __launch_bounds__(soc_block_size) scatter_sparse_hessian_into_au
  *
  * `Hs_diag` is left populated for downstream matrix-free matvec / iterative refinement use.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void scatter_sparse_hessian_into_augmented(cone_data_t<i_t, f_t>& cones,
                                            rmm::device_uvector<f_t>& augmented_x,
                                            rmm::device_uvector<f_t>& Hs_diag,
@@ -1083,6 +1083,89 @@ void scatter_sparse_hessian_into_augmented(cone_data_t<i_t, f_t>& cones,
   RAFT_CUDA_TRY(cudaPeekAtLastError());
 }
 
+// Write the cone block a cold start factorizes, before any Nesterov-Todd scaling
+// exists: Hessian diagonal -Q - dual_perturb, and sparse expansion couplings and
+// diagonals at zero. The rank-2 matvec state is cleared to match that matrix.
+template <typename i_t, typename f_t>
+__global__ void __launch_bounds__(soc_block_size)
+  restore_initial_sparse_cone_block_kernel(raft::device_span<f_t> augmented_x,
+                                           raft::device_span<f_t> Hs_diag,
+                                           raft::device_span<f_t> sparse_v,
+                                           raft::device_span<f_t> sparse_u,
+                                           raft::device_span<f_t> d,
+                                           raft::device_span<const i_t> sparse_entry_offsets,
+                                           i_t n_sparse_cones,
+                                           raft::device_span<const i_t> hessian_diag_csr_indices,
+                                           raft::device_span<const f_t> q_values,
+                                           raft::device_span<const i_t> exp_v_col,
+                                           raft::device_span<const i_t> exp_u_col,
+                                           raft::device_span<const i_t> exp_v_row,
+                                           raft::device_span<const i_t> exp_u_row,
+                                           raft::device_span<const i_t> sparse_expansion_D,
+                                           f_t dual_perturb)
+{
+  const size_t idx = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
+  if (idx >= Hs_diag.size()) { return; }
+
+  const i_t idx_i      = static_cast<i_t>(idx);
+  const i_t sparse_idx = bucket_index(sparse_entry_offsets, idx_i, n_sparse_cones);
+  const bool is_head   = idx_i == sparse_entry_offsets[sparse_idx];
+
+  Hs_diag[idx]  = 0.0;
+  sparse_v[idx] = 0.0;
+  sparse_u[idx] = 0.0;
+
+  augmented_x[hessian_diag_csr_indices[idx]] = -q_values[idx] - dual_perturb;
+  augmented_x[exp_v_col[idx]]                = 0.0;
+  augmented_x[exp_u_col[idx]]                = 0.0;
+  augmented_x[exp_v_row[idx]]                = 0.0;
+  augmented_x[exp_u_row[idx]]                = 0.0;
+
+  if (is_head) {
+    d[sparse_idx]                                       = 0.0;
+    augmented_x[sparse_expansion_D[2 * sparse_idx]]     = 0.0;
+    augmented_x[sparse_expansion_D[2 * sparse_idx + 1]] = 0.0;
+  }
+}
+
+template <typename i_t, typename f_t>
+void restore_initial_sparse_cone_block(cone_data_t<i_t, f_t>& cones,
+                                       rmm::device_uvector<f_t>& augmented_x,
+                                       rmm::device_uvector<f_t>& Hs_diag,
+                                       const rmm::device_uvector<i_t>& hessian_diag_csr_indices,
+                                       const rmm::device_uvector<f_t>& q_values,
+                                       const rmm::device_uvector<i_t>& exp_v_col,
+                                       const rmm::device_uvector<i_t>& exp_u_col,
+                                       const rmm::device_uvector<i_t>& exp_v_row,
+                                       const rmm::device_uvector<i_t>& exp_u_row,
+                                       const rmm::device_uvector<i_t>& sparse_expansion_D,
+                                       cuda::stream_ref stream,
+                                       f_t dual_perturb)
+{
+  if (!cones.has_sparse_cones()) { return; }
+
+  const i_t n_sparse      = cones.n_sparse_cones;
+  const size_t E          = cones.n_sparse_cone_entries;
+  const size_t entry_grid = raft::ceildiv<size_t>(E, soc_block_size);
+  restore_initial_sparse_cone_block_kernel<i_t, f_t>
+    <<<entry_grid, soc_block_size, 0, stream.get()>>>(cuopt::make_span(augmented_x),
+                                                      cuopt::make_span(Hs_diag),
+                                                      cuopt::make_span(cones.sparse_v),
+                                                      cuopt::make_span(cones.sparse_u),
+                                                      cuopt::make_span(cones.d),
+                                                      cuopt::make_span(cones.sparse_entry_offsets),
+                                                      n_sparse,
+                                                      cuopt::make_span(hessian_diag_csr_indices),
+                                                      cuopt::make_span(q_values),
+                                                      cuopt::make_span(exp_v_col),
+                                                      cuopt::make_span(exp_u_col),
+                                                      cuopt::make_span(exp_v_row),
+                                                      cuopt::make_span(exp_u_row),
+                                                      cuopt::make_span(sparse_expansion_D),
+                                                      dual_perturb);
+  RAFT_CUDA_TRY(cudaPeekAtLastError());
+}
+
 /**
  * Accumulate the sparse-SOC expanded KKT block into a matrix-free product.
  *
@@ -1095,7 +1178,7 @@ void scatter_sparse_hessian_into_augmented(cone_data_t<i_t, f_t>& cones,
  *
  * `v` and `u` are the rank-2 vectors in cones.sparse_v / cones.sparse_u.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void sparse_augmented_matvec_kernel(raft::device_span<const f_t> x,
                                                raft::device_span<f_t> r1,
                                                raft::device_span<f_t> y_exp,
@@ -1126,8 +1209,8 @@ __global__ void sparse_augmented_matvec_kernel(raft::device_span<const f_t> x,
   const f_t x_exp_u = x[exp_u];
   const f_t eta_sq  = eta[cone] * eta[cone];
 
-  f_t partial_dot_v = f_t(0);
-  f_t partial_dot_u = f_t(0);
+  f_t partial_dot_v = 0.0;
+  f_t partial_dot_u = 0.0;
   for (i_t j = threadIdx.x; j < q; j += blockDim.x) {
     const f_t xj = x[base + j];
     const f_t vj = sparse_v[flat + j];
@@ -1162,7 +1245,7 @@ __global__ void sparse_augmented_matvec_kernel(raft::device_span<const f_t> x,
  * `sparse_augmented_matvec_kernel`). Accumulates the cone-row contribution into
  * `r1` and the expansion-row contribution into `y_exp`, one block per sparse cone.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void launch_sparse_augmented_matvec(raft::device_span<const f_t> x,
                                     raft::device_span<f_t> r1,
                                     raft::device_span<f_t> y_exp,
@@ -1201,7 +1284,7 @@ void launch_sparse_augmented_matvec(raft::device_span<const f_t> x,
 
 // Scatter one nonzero of a dense cone's q x q NT Hessian block into the
 // augmented value buffer. Only cones with dim <= soc_threshold take this path.
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void __launch_bounds__(soc_block_size)
   scatter_dense_hessian_into_augmented_kernel(raft::device_span<f_t> augmented_x,
                                               raft::device_span<const i_t> csr_indices,
@@ -1230,7 +1313,7 @@ __global__ void __launch_bounds__(soc_block_size)
   const f_t w0     = w[off];
   const f_t u_r    = (r == 0) ? w0 : w[off + r];
   const f_t u_c    = (c == 0) ? w0 : w[off + c];
-  const f_t val    = f_t{2} * u_r * eta_sq * u_c;
+  const f_t val    = 2.0 * u_r * eta_sq * u_c;
 
   f_t entry = -val - q_values[e];
   if (r == c) {
@@ -1244,7 +1327,7 @@ __global__ void __launch_bounds__(soc_block_size)
  * Write the full q x q NT Hessian blocks of the dense cones (dim <= soc_threshold)
  * into the augmented system value buffer at the precomputed CSR positions.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void scatter_dense_hessian_into_augmented(const cone_data_t<i_t, f_t>& cones,
                                           rmm::device_uvector<f_t>& augmented_x,
                                           const rmm::device_uvector<i_t>& csr_indices,
@@ -1283,7 +1366,7 @@ void scatter_dense_hessian_into_augmented(const cone_data_t<i_t, f_t>& cones,
 // boundary. Small/medium/large partitions come from segmented_sum_t.
 // =============================================================================
 
-template <std::floating_point f_t>
+template <typename f_t>
 HD f_t cone_step_length_from_scalars(
   f_t u0, f_t du0, f_t du_tail_sq, f_t u_tail_du_tail, f_t u_tail_sq, f_t alpha_max)
 {
@@ -1318,7 +1401,7 @@ HD f_t cone_step_length_from_scalars(
  * One warp per small cone: accumulate the three tail scalars, then solve for
  * alpha[cone]. Tail-only (local index 0 is the SOC head).
  */
-template <std::integral i_t, std::floating_point f_t, int warps_per_cta = 8>
+template <typename i_t, typename f_t, int warps_per_cta = 8>
 __global__ void __launch_bounds__(warps_per_cta* raft::WarpSize)
   step_length_small_kernel(raft::device_span<const f_t> u,
                            raft::device_span<const f_t> du,
@@ -1366,7 +1449,7 @@ __global__ void __launch_bounds__(warps_per_cta* raft::WarpSize)
 /**
  * One block per medium cone: three-scalar tail reduction + step solve.
  */
-template <std::integral i_t, std::floating_point f_t, int block_dim = 256>
+template <typename i_t, typename f_t, int block_dim = 256>
 __global__ void __launch_bounds__(block_dim)
   step_length_medium_kernel(raft::device_span<const f_t> u,
                             raft::device_span<const f_t> du,
@@ -1409,9 +1492,9 @@ __global__ void __launch_bounds__(block_dim)
     }
   }
 
-  f_t du_sq = f_t{0};
-  f_t u_du  = f_t{0};
-  f_t u_sq  = f_t{0};
+  f_t du_sq = 0.0;
+  f_t u_du  = 0.0;
+  f_t u_sq  = 0.0;
 #pragma unroll
   for (int k = 0; k < items_per_thread; ++k) {
     du_sq += acc_du_sq[k];
@@ -1430,7 +1513,7 @@ __global__ void __launch_bounds__(block_dim)
   }
 }
 
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 __global__ void step_length_large_solve_kernel(raft::device_span<const f_t> u,
                                                raft::device_span<const f_t> du,
                                                raft::device_span<f_t> alpha,
@@ -1453,7 +1536,7 @@ __global__ void step_length_large_solve_kernel(raft::device_span<const f_t> u,
  * Size-aware step length: one pass over each cone's tail, then write
  * alpha[cone].
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void launch_cone_step_length(segmented_sum_t<i_t>& partitions,
                              raft::device_span<const f_t> u,
                              raft::device_span<const f_t> du,
@@ -1536,7 +1619,7 @@ void launch_cone_step_length(segmented_sum_t<i_t>& partitions,
  *
  *   x + alpha dx in Q,  z + alpha dz in Q,  alpha <= alpha_max.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 f_t compute_cone_step_length(cone_data_t<i_t, f_t>& cones,
                              raft::device_span<const f_t> dx,
                              raft::device_span<const f_t> dz,
@@ -1578,7 +1661,7 @@ f_t compute_cone_step_length(cone_data_t<i_t, f_t>& cones,
  * On return, `out` holds `q`. Internally, `out` is reused for `W^{-T} dz_aff` and
  * then `d`; `scratch.temp_cone` is reused for `W dx_aff`, then `-p`.
  */
-template <std::integral i_t, std::floating_point f_t>
+template <typename i_t, typename f_t>
 void compute_combined_cone_rhs_term(raft::device_span<const f_t> dx_aff,
                                     raft::device_span<const f_t> dz_aff,
                                     cone_data_t<i_t, f_t>& cones,
diff --git a/cpp/src/barrier/translate_soc.hpp b/cpp/src/barrier/translate_soc.hpp
index f89994b949..4f6f251fe7 100644
--- a/cpp/src/barrier/translate_soc.hpp
+++ b/cpp/src/barrier/translate_soc.hpp
@@ -69,10 +69,14 @@ void convert_quadratic_constraints_to_second_order_cones(
   // Use a practical tolerance for text-parsed MPS numeric values.
   const f_t tol = std::numeric_limits<f_t>::epsilon() * 2;
 
+  // Row number before second-order cone translation.
+  user_problem.original_num_rows = csr_A.m;
+
   // Derive implied lower bounds from singleton inequality rows.
   // Used to check if SOC head variables have implied non-negativity from the constraint system
   // without actually modifying the variable bounds (which would add barrier terms).
   std::vector<f_t> implied_lower(n, -std::numeric_limits<f_t>::infinity());
+  std::vector<std::vector<std::pair<i_t, f_t>>> implied_rows(n);
   for (i_t i = 0; i < csr_A.m; i++) {
     const i_t row_start = csr_A.row_start[i];
     const i_t row_end   = csr_A.row_start[i + 1];
@@ -85,11 +89,21 @@ void convert_quadratic_constraints_to_second_order_cones(
     const f_t bound = b / a;
     if (sense == 'G' && a > 0) {
       implied_lower[j] = std::max(implied_lower[j], bound);
+      implied_rows[j].emplace_back(i, a);
     } else if (sense == 'L' && a < 0) {
       implied_lower[j] = std::max(implied_lower[j], bound);
+      implied_rows[j].emplace_back(i, a);
     }
   }
 
+  auto save_cone_head_bound = [&](i_t head) {
+    if (user_problem.lower[head] >= 0) { return; }
+    simplex::cone_head_bound_t<i_t, f_t> bound;
+    bound.head_col = head;
+    bound.rows     = implied_rows[head];
+    user_problem.cone_head_bounds.push_back(std::move(bound));
+  };
+
   // SOC conversion routes each quadratic constraint as follows:
   //
   // Fast path (pattern-matched, rhs = 0, no linear part in COLUMNS):
@@ -376,6 +390,7 @@ void convert_quadratic_constraints_to_second_order_cones(
                       "non-negative lower bound for the constraint to be convex",
                       qc.constraint_row_name.c_str(),
                       static_cast<int>(head));
+        save_cone_head_bound(head);
         cone.reserve(q_nnz);
         cone.push_back(head);
         cone.insert(cone.end(), tail_vars.begin(), tail_vars.end());
@@ -417,6 +432,7 @@ void convert_quadratic_constraints_to_second_order_cones(
                       "non-negative lower bound for the constraint to be convex",
                       qc.constraint_row_name.c_str(),
                       static_cast<int>(a));
+        save_cone_head_bound(a);
         cuopt_expects(std::max(user_problem.lower[b], implied_lower[b]) >= 0,
                       error_type_t::ValidationError,
                       "Quadratic constraint '%s': rotated second-order cone head variable (index "
@@ -424,6 +440,7 @@ void convert_quadratic_constraints_to_second_order_cones(
                       "non-negative lower bound for the constraint to be convex",
                       qc.constraint_row_name.c_str(),
                       static_cast<int>(b));
+        save_cone_head_bound(b);
         rotated_cones.push_back(rotated_soc_t{a, b, tail_vars, false, head_lift_sqrt_ratio});
       }
 
diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp
index d2731999ac..b053645e74 100644
--- a/cpp/src/dual_simplex/presolve.cpp
+++ b/cpp/src/dual_simplex/presolve.cpp
@@ -1862,9 +1862,11 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,
   }
   problem.Q.check_matrix("Before free variable expansion");
 
-  // Free linear variables. We handle them directly in QP/SOCP or split them in LP.
+  // Free linear variables. QP/SOCP keep them. A sequence LP does too: the v-w split would
+  // change the cached columns. barrier_eliminate_free_variables is false only on that path.
   const bool direct_free_linear =
-    settings.barrier_presolve && free_variables > 0 && (problem.Q.n > 0 || has_cones);
+    settings.barrier_presolve && free_variables > 0 &&
+    (problem.Q.n > 0 || has_cones || !settings.barrier_eliminate_free_variables);
   if (direct_free_linear) {
     presolve_info.free_variable_pairs.clear();
     presolve_info.direct_free_variables.clear();
@@ -1963,8 +1965,9 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,
     settings.log.printf("Dependent row check in %.2fs\n", toc(dependent_row_start));
   }
 
-  // LP already goes through PSLP; this substitution is for QP/SOCP only.
-  if (settings.barrier_presolve && (has_cones || problem.Q.n > 0)) {
+  // LP already goes through PSLP; this substitution is for QP/SOCP only. A cache fill skips it.
+  if (settings.barrier_presolve && settings.barrier_eliminate_free_variables &&
+      (has_cones || problem.Q.n > 0)) {
     const i_t old_free_count         = static_cast<i_t>(presolve_info.direct_free_variables.size());
     const f_t free_elimination_start = tic();
     const i_t pivot_rejected         = eliminate_free_variables(problem, presolve_info);
diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp
index 63d69aa42b..123e87efc3 100644
--- a/cpp/src/dual_simplex/presolve.hpp
+++ b/cpp/src/dual_simplex/presolve.hpp
@@ -40,7 +40,7 @@ row_bounds_t<f_t> get_range_bounds_from_sense(char row_sense, f_t rhs, f_t range
 
 template <typename i_t, typename f_t>
 struct lp_problem_t {
-  lp_problem_t(raft::handle_t const* handle_ptr_, i_t m, i_t n, i_t nz)
+  lp_problem_t(raft::handle_t const* handle_ptr_, i_t m, i_t n, i_t nz, i_t cone_var_start_ = 0)
     : handle_ptr(handle_ptr_),
       num_rows(m),
       num_cols(n),
@@ -50,7 +50,8 @@ struct lp_problem_t {
       rhs(m),
       lower(n),
       upper(n),
-      obj_constant(0.0)
+      obj_constant(0.0),
+      cone_var_start(cone_var_start_)
   {
   }
   raft::handle_t const* handle_ptr;
diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp
index 79fd9176e3..f352a88189 100644
--- a/cpp/src/dual_simplex/simplex_solver_settings.hpp
+++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp
@@ -68,6 +68,7 @@ struct simplex_solver_settings_t {
       eliminate_singletons(true),
       print_presolve_stats(true),
       barrier_presolve(false),
+      barrier_eliminate_free_variables(true),
       cudss_deterministic(false),
       deterministic(false),
       barrier(false),
@@ -172,11 +173,12 @@ struct simplex_solver_settings_t {
   bool relaxation;                 // true to only solve the LP relaxation of a MIP
   bool
     use_left_looking_lu;  // true to use left looking LU factorization, false to use right looking
-  bool eliminate_singletons;  // true to eliminate singletons from the basis
-  bool print_presolve_stats;  // true to print presolve stats
-  bool barrier_presolve;      // true to use barrier presolve
-  bool cudss_deterministic;   // true to use cuDSS deterministic mode, false for non-deterministic
-  bool barrier;               // true to use barrier method, false to use dual simplex method
+  bool eliminate_singletons;              // true to eliminate singletons from the basis
+  bool print_presolve_stats;              // true to print presolve stats
+  bool barrier_presolve;                  // true to use barrier presolve
+  bool barrier_eliminate_free_variables;  // true to eliminate zero-cost free variables
+  bool cudss_deterministic;  // true to use cuDSS deterministic mode, false for non-deterministic
+  bool barrier;              // true to use barrier method, false to use dual simplex method
   bool deterministic;  // true to use B&B deterministic mode, false to use non-deterministic mode
   bool eliminate_dense_columns;         // true to eliminate dense columns from A*D*A^T
   bool barrier_iterative_refinement;    // true to use iterative refinement for barrier method
diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp
index aac7db8682..35b523ea10 100644
--- a/cpp/src/dual_simplex/solve.cpp
+++ b/cpp/src/dual_simplex/solve.cpp
@@ -49,6 +49,7 @@ void unscale_uncrush_barrier_to_user(const user_problem_t<i_t, f_t>& user_proble
                                      const raft::handle_t* handle_ptr,
                                      i_t original_num_rows,
                                      i_t original_num_cols,
+                                     i_t converted_cone_var_start,
                                      const lp_problem_t<i_t, f_t>& barrier_lp,
                                      const presolve_info_t<i_t, f_t>& presolve_info,
                                      const std::vector<f_t>& column_scales,
@@ -69,8 +70,10 @@ void unscale_uncrush_barrier_to_user(const user_problem_t<i_t, f_t>& user_proble
                              unscaled_y,
                              unscaled_z);
 
-  // Dummy converted LP: sizes only. Bound-free=0 so uncrush_solution never reads A.
-  lp_problem_t<i_t, f_t> converted(handle_ptr, original_num_rows, original_num_cols, 0);
+  // Dummy converted LP. Bound-free=0 so uncrush_solution never reads A. cone_var_start undoes
+  // the shift convert applied to the cone block.
+  lp_problem_t<i_t, f_t> converted(
+    handle_ptr, original_num_rows, original_num_cols, 0, converted_cone_var_start);
   lp_solution_t<i_t, f_t> lp_solution(original_num_rows, original_num_cols);
   uncrush_solution(presolve_info,
                    barrier_settings,
@@ -522,20 +525,29 @@ lp_status_t solve_linear_program_with_barrier(
   const simplex_solver_settings_t<i_t, f_t>& settings,
   f_t start_time,
   lp_solution_t<i_t, f_t>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache,
+  cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache,
   const raft::handle_t* handle_ptr)
 {
   lp_status_t status                                   = lp_status_t::UNSET;
   simplex_solver_settings_t<i_t, f_t> barrier_settings = settings;
 
-  auto const* xf = (cache != nullptr && cache->dirty()) ? cache->transform() : nullptr;
-  const bool reuse_cached_data = cuopt::mathematical_optimization::can_reuse_barrier_cache(
-    xf,
-    settings.barrier_presolve_bound_free_variables,
-    user_problem.num_cols,
-    user_problem.num_rows,
-    !user_problem.Q_values.empty(),
-    !user_problem.second_order_cone_dims.empty());
+  auto const* xf          = (cache != nullptr && cache->dirty()) ? cache->transform() : nullptr;
+  const bool user_has_soc = !user_problem.second_order_cone_dims.empty();
+  // run_barrier already resolved -1 to 0. can_reuse_barrier_cache also rejects a cache whose
+  // presolve bounded free variables, which this path cannot replay.
+  const bool reuse_cached_data =
+    cuopt::mathematical_optimization::can_reuse_barrier_cache(
+      xf,
+      settings.barrier_presolve_bound_free_variables,
+      user_problem.num_cols,
+      user_problem.num_rows,
+      !user_problem.Q_values.empty(),
+      user_has_soc) &&
+    // Cone reuse also requires the expanded layout and dimensions to match. These checks stay
+    // here because the other caller of can_reuse_barrier_cache sees the pre-expansion problem.
+    (xf != nullptr && cuopt::mathematical_optimization::cone_layout_matches(*xf, user_problem)) &&
+    (!user_has_soc ||
+     (user_problem.num_cols == xf->user_num_cols && user_problem.num_rows == xf->user_num_rows));
 
   if (reuse_cached_data) {
     if (cache->rhs_infeasible()) {
@@ -553,6 +565,7 @@ lp_status_t solve_linear_program_with_barrier(
                                       cache->handle_ptr(),
                                       xf->original_num_rows,
                                       xf->original_num_cols,
+                                      xf->converted_cone_var_start,
                                       *xf->barrier_lp,
                                       xf->presolve_info,
                                       xf->column_scales,
@@ -614,9 +627,9 @@ lp_status_t solve_linear_program_with_barrier(
   lp_problem_t<i_t, f_t> const* solver_lp = &barrier_lp;
   if (cache != nullptr) {
     cache->clear();
-    auto xf           = std::make_unique<cuopt::mathematical_optimization::barrier_transform_t>();
-    xf->user_num_cols = user_problem.num_cols;
-    xf->user_num_rows = user_problem.num_rows;
+    auto xf = std::make_unique<cuopt::mathematical_optimization::barrier_transform_t<i_t, f_t>>();
+    xf->user_num_cols                = user_problem.num_cols;
+    xf->user_num_rows                = user_problem.num_rows;
     xf->original_num_cols            = original_lp.num_cols;
     xf->original_num_rows            = original_lp.num_rows;
     xf->obj_scale                    = user_problem.obj_scale;
@@ -624,13 +637,23 @@ lp_status_t solve_linear_program_with_barrier(
     xf->row_sense                    = user_problem.row_sense;
     xf->cone_var_start               = user_problem.cone_var_start;
     xf->second_order_cone_dims       = user_problem.second_order_cone_dims;
-    xf->expanded_original_num_cols   = user_problem.original_num_cols;
+    xf->pre_expansion_num_cols       = user_problem.original_num_cols;
     xf->original_col_to_expanded_col = user_problem.original_col_to_expanded_col;
-    xf->presolve_info                = presolve_info;
-    xf->column_scales                = column_scales;
-    xf->row_scales                   = row_scales;
+    xf->pre_expansion_num_rows       = user_problem.original_num_rows;
+    xf->converted_cone_var_start     = original_lp.cone_var_start;
+    xf->cone_head_bounds             = user_problem.cone_head_bounds;
+    // Append the expanded rhs.
+    if (!user_problem.original_col_to_expanded_col.empty()) {
+      xf->cone_row_rhs.assign(user_problem.rhs.begin() + user_problem.original_num_rows,
+                              user_problem.rhs.end());
+    }
+    xf->presolve_info = presolve_info;
+    xf->column_scales = column_scales;
+    xf->row_scales    = row_scales;
     // convert_range_rows zeroes rhs[i] onto the slack bounds and folding aggregates rows, so
     // neither leaves the user RHS in barrier_lp->rhs. Plain inequality/equality slacks do.
+    // Alias rows (alias - x = 0) are appended after the original constraints and stay in the
+    // cached tail, so they do not block an update of the original RHS.
     xf->rhs_update_supported =
       user_problem.num_range_rows == 0 && !presolve_info.folding_info.is_folded;
     xf->barrier_lp = std::make_unique<lp_problem_t<i_t, f_t>>(barrier_lp);
@@ -659,10 +682,15 @@ lp_status_t solve_linear_program_with_barrier(
       auto* xf = cache->transform();
       // Crushing this solve's own data also checks the maps still describe it: presolve may
       // have dualized an LP, in which case the crush throws.
+      // Both crushes take model coordinates. The expansion only appends rows, so the RHS prefix
+      // is already model-sized, but it permutes columns, so gather the objective back.
       bool objective_shift_ok = false;
       try {
+        auto const problem_objective =
+          cuopt::mathematical_optimization::gather_problem_objective<i_t>(*xf,
+                                                                          user_problem.objective);
         auto const crushed = cuopt::mathematical_optimization::crush_user_linear_objective(
-          *xf, user_problem.objective.data(), user_problem.num_cols);
+          *xf, problem_objective.data(), static_cast<i_t>(problem_objective.size()));
         objective_shift_ok =
           shift_from<i_t, f_t>(solver_lp->objective, crushed, xf->linear_obj_shift) == 0;
       } catch (std::exception const&) {
@@ -675,11 +703,11 @@ lp_status_t solve_linear_program_with_barrier(
       if (xf->rhs_update_supported) {
         std::vector<f_t> crushed;
         std::string error;
-        bool const rhs_shift_ok =
-          cuopt::mathematical_optimization::crush_user_rhs(
-            *xf, user_problem.rhs.data(), user_problem.num_rows, crushed, error) ==
-            cuopt::mathematical_optimization::crush_rhs_status_t::success &&
-          shift_from<i_t, f_t>(solver_lp->rhs, crushed, xf->rhs_shift) == 0;
+        const i_t problem_m     = cuopt::mathematical_optimization::problem_num_rows(*xf);
+        bool const rhs_shift_ok = cuopt::mathematical_optimization::crush_user_rhs(
+                                    *xf, user_problem.rhs.data(), problem_m, crushed, error) ==
+                                    cuopt::mathematical_optimization::crush_rhs_status_t::success &&
+                                  shift_from<i_t, f_t>(solver_lp->rhs, crushed, xf->rhs_shift) == 0;
         if (!rhs_shift_ok) {
           // Maps cannot reproduce this RHS, so refuse later updates.
           xf->rhs_update_supported = false;
@@ -975,7 +1003,7 @@ lp_status_t solve_linear_program_with_barrier(
   const user_problem_t<i_t, f_t>& user_problem,
   const simplex_solver_settings_t<i_t, f_t>& settings,
   lp_solution_t<i_t, f_t>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache)
+  cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache)
 {
   f_t start_time = tic();
   return solve_linear_program_with_barrier(
@@ -988,7 +1016,7 @@ lp_status_t solve_linear_program_with_barrier(
   const simplex_solver_settings_t<i_t, f_t>& settings,
   f_t start_time,
   lp_solution_t<i_t, f_t>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache)
+  cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache)
 {
   return solve_linear_program_with_barrier(
     user_problem, settings, start_time, solution, cache, user_problem.handle_ptr);
@@ -1244,7 +1272,7 @@ template lp_status_t solve_linear_program_with_barrier(
   const user_problem_t<int, double>& user_problem,
   const simplex_solver_settings_t<int, double>& settings,
   lp_solution_t<int, double>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache);
+  cuopt::mathematical_optimization::barrier_cache_t<int, double>* cache);
 
 template lp_status_t solve_linear_program_with_primal(
   const user_problem_t<int, double>& user_problem,
@@ -1257,14 +1285,14 @@ template lp_status_t solve_linear_program_with_barrier(
   const simplex_solver_settings_t<int, double>& settings,
   double start_time,
   lp_solution_t<int, double>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache);
+  cuopt::mathematical_optimization::barrier_cache_t<int, double>* cache);
 
 template lp_status_t solve_linear_program_with_barrier(
   const user_problem_t<int, double>& user_problem,
   const simplex_solver_settings_t<int, double>& settings,
   double start_time,
   lp_solution_t<int, double>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache,
+  cuopt::mathematical_optimization::barrier_cache_t<int, double>* cache,
   const raft::handle_t* handle_ptr);
 
 template lp_status_t solve_linear_program(const user_problem_t<int, double>& user_problem,
diff --git a/cpp/src/dual_simplex/solve.hpp b/cpp/src/dual_simplex/solve.hpp
index 08ae11a9dc..18f8ed0624 100644
--- a/cpp/src/dual_simplex/solve.hpp
+++ b/cpp/src/dual_simplex/solve.hpp
@@ -18,6 +18,7 @@ struct work_limit_context_t;
 }
 
 namespace cuopt::mathematical_optimization {
+template <typename i_t, typename f_t>
 class barrier_cache_t;
 }  // namespace cuopt::mathematical_optimization
 
@@ -155,7 +156,7 @@ lp_status_t solve_linear_program_with_barrier(
   const user_problem_t<i_t, f_t>& user_problem,
   const simplex_solver_settings_t<i_t, f_t>& settings,
   lp_solution_t<i_t, f_t>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache = nullptr);
+  cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache = nullptr);
 
 template <typename i_t, typename f_t>
 lp_status_t solve_linear_program_with_barrier(
@@ -163,7 +164,7 @@ lp_status_t solve_linear_program_with_barrier(
   const simplex_solver_settings_t<i_t, f_t>& settings,
   f_t start_time,
   lp_solution_t<i_t, f_t>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache = nullptr);
+  cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache = nullptr);
 
 template <typename i_t, typename f_t>
 lp_status_t solve_linear_program_with_primal(const user_problem_t<i_t, f_t>& user_problem,
@@ -176,7 +177,7 @@ lp_status_t solve_linear_program_with_barrier(
   const simplex_solver_settings_t<i_t, f_t>& settings,
   f_t start_time,
   lp_solution_t<i_t, f_t>& solution,
-  cuopt::mathematical_optimization::barrier_cache_t* cache,
+  cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache,
   const raft::handle_t* handle_ptr);
 
 template <typename i_t, typename f_t>
diff --git a/cpp/src/dual_simplex/user_problem.hpp b/cpp/src/dual_simplex/user_problem.hpp
index 419b33e1e1..4195e4d37a 100644
--- a/cpp/src/dual_simplex/user_problem.hpp
+++ b/cpp/src/dual_simplex/user_problem.hpp
@@ -14,6 +14,7 @@
 #include <raft/core/handle.hpp>
 
 #include <string>
+#include <utility>
 
 namespace cuopt::mathematical_optimization::simplex {
 
@@ -33,6 +34,15 @@ struct objective_step_t {
   bool has_step() const { return step_size > 0; }
 };
 
+// Singleton rows that force a cone head nonnegative. Recorded while the expansion proves the
+// head, so an RHS update can re-check them.
+template <typename i_t, typename f_t>
+struct cone_head_bound_t {
+  i_t head_col{0};
+  // (row, coefficient) pairs, each implying head >= rhs[row] / coefficient.
+  std::vector<std::pair<i_t, f_t>> rows;
+};
+
 template <typename i_t, typename f_t>
 struct user_problem_t {
   user_problem_t(raft::handle_t const* handle_ptr_)
@@ -69,6 +79,10 @@ struct user_problem_t {
   // expanded layout (num_cols) and must be projected back via original_col_to_expanded_col.
   i_t original_num_cols{0};
   std::vector<i_t> original_col_to_expanded_col;
+  // Row count before QCMATRIX->SOC expansion. The expansion only appends rows, so rows
+  // [0, original_num_rows) still hold the model's own constraints at their original indices.
+  i_t original_num_rows{0};
+  std::vector<cone_head_bound_t<i_t, f_t>> cone_head_bounds;
 };
 
 }  // namespace cuopt::mathematical_optimization::simplex
diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu
index 33d26f4641..1ddf20bfe9 100644
--- a/cpp/src/pdlp/solve.cu
+++ b/cpp/src/pdlp/solve.cu
@@ -83,12 +83,13 @@ template <typename i_t, typename f_t>
 simplex::user_problem_t<i_t, f_t> user_problem_from_transform(
   raft::handle_t const* handle_ptr,
   optimization_problem_t<i_t, f_t>& model,
-  cuopt::mathematical_optimization::barrier_transform_t const& xf)
+  cuopt::mathematical_optimization::barrier_transform_t<i_t, f_t> const& xf)
 {
   simplex::user_problem_t<i_t, f_t> user_problem(handle_ptr);
-  user_problem.num_rows  = xf.user_num_rows;
-  user_problem.num_cols  = xf.user_num_cols;
-  user_problem.objective = model.get_objective_coefficients_host();
+  user_problem.num_rows = xf.user_num_rows;
+  user_problem.num_cols = xf.user_num_cols;
+  user_problem.objective =
+    scatter_problem_objective<i_t>(xf, model.get_objective_coefficients_host());
   user_problem.row_sense = xf.row_sense;
   user_problem.rhs.assign(static_cast<std::size_t>(xf.user_num_rows), f_t(0));
   user_problem.obj_scale    = static_cast<f_t>(xf.obj_scale);
@@ -97,8 +98,9 @@ simplex::user_problem_t<i_t, f_t> user_problem_from_transform(
   user_problem.Q_values.assign(1, f_t(1));
   user_problem.cone_var_start               = xf.cone_var_start;
   user_problem.second_order_cone_dims       = xf.second_order_cone_dims;
-  user_problem.original_num_cols            = xf.expanded_original_num_cols;
+  user_problem.original_num_cols            = xf.pre_expansion_num_cols;
   user_problem.original_col_to_expanded_col = xf.original_col_to_expanded_col;
+  user_problem.original_num_rows            = xf.pre_expansion_num_rows;
   return user_problem;
 }
 
@@ -530,48 +532,62 @@ i_t effective_bound_free_variables(pdlp_solver_settings_t<i_t, f_t> const& setti
            : settings.barrier_presolve_bound_free_variables;
 }
 
+// Equality substitution stores a pivot RHS the reuse path cannot refresh, so a sequence solve
+// leaves those free columns in the barrier system.
+template <typename i_t, typename f_t>
+bool effective_eliminate_free_variables(pdlp_solver_settings_t<i_t, f_t> const& settings)
+{
+  return !settings.sequence_solve;
+}
+
 template <typename i_t, typename f_t>
 std::tuple<simplex::lp_solution_t<i_t, f_t>, simplex::lp_status_t, f_t, f_t, f_t> run_barrier(
   const simplex::user_problem_t<i_t, f_t>& user_problem,
   pdlp_solver_settings_t<i_t, f_t> const& settings,
   const timer_t& timer,
   const raft::handle_t* handle_ptr,
-  cuopt::mathematical_optimization::barrier_cache_t* cache = nullptr)
+  cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache = nullptr)
 {
   f_t norm_user_objective = vector_norm2<i_t, f_t>(user_problem.objective);
   f_t norm_rhs            = vector_norm2<i_t, f_t>(user_problem.rhs);
 
   simplex::simplex_solver_settings_t<i_t, f_t> barrier_settings;
-  barrier_settings.num_gpus                              = settings.num_gpus;
-  barrier_settings.time_limit                            = settings.time_limit;
-  barrier_settings.iteration_limit                       = settings.iteration_limit;
-  barrier_settings.concurrent_halt                       = settings.concurrent_halt;
-  barrier_settings.initial_perturbation                  = settings.initial_perturbation;
-  barrier_settings.remove_perturbation                   = settings.remove_perturbation;
-  barrier_settings.primal_pricing                        = settings.primal_pricing;
-  barrier_settings.folding                               = settings.folding;
+  barrier_settings.num_gpus             = settings.num_gpus;
+  barrier_settings.time_limit           = settings.time_limit;
+  barrier_settings.iteration_limit      = settings.iteration_limit;
+  barrier_settings.concurrent_halt      = settings.concurrent_halt;
+  barrier_settings.initial_perturbation = settings.initial_perturbation;
+  barrier_settings.remove_perturbation  = settings.remove_perturbation;
+  barrier_settings.primal_pricing       = settings.primal_pricing;
+  barrier_settings.folding              = settings.folding;
+  // Folding rewrites rows, so a sequence LP cannot replay an RHS update through the cache.
+  if (settings.sequence_solve && user_problem.Q_values.empty() &&
+      user_problem.second_order_cone_dims.empty()) {
+    barrier_settings.folding = 0;
+  }
   barrier_settings.augmented                             = settings.augmented;
   barrier_settings.dualize                               = settings.dualize;
   barrier_settings.ordering                              = settings.ordering;
   barrier_settings.barrier_dual_initial_point            = settings.barrier_dual_initial_point;
   barrier_settings.postsolve_info                        = settings.postsolve_info;
   barrier_settings.barrier_presolve_bound_free_variables = effective_bound_free_variables(settings);
-  barrier_settings.barrier_initial_point_safeguard       = settings.barrier_initial_point_safeguard;
-  barrier_settings.barrier                               = true;
-  barrier_settings.barrier_presolve                      = true;
-  barrier_settings.crossover                             = settings.crossover;
-  barrier_settings.eliminate_dense_columns               = settings.eliminate_dense_columns;
-  barrier_settings.barrier_iterative_refinement          = settings.barrier_iterative_refinement;
-  barrier_settings.barrier_adaptive_regularization       = settings.barrier_adaptive_regularization;
-  barrier_settings.barrier_primal_regularization         = settings.barrier_primal_regularization;
-  barrier_settings.barrier_dual_regularization           = settings.barrier_dual_regularization;
-  barrier_settings.barrier_soc_threshold                 = settings.barrier_soc_threshold;
-  barrier_settings.barrier_step_scale                    = settings.barrier_step_scale;
-  barrier_settings.qcqp_ruiz_equilibration               = settings.qcqp_ruiz_equilibration;
-  barrier_settings.gpu_ruiz_nnz_threshold                = settings.gpu_ruiz_nnz_threshold;
-  barrier_settings.cudss_deterministic                   = settings.cudss_deterministic;
-  barrier_settings.barrier_relaxed_feasibility_tol = settings.tolerances.relative_primal_tolerance;
-  barrier_settings.barrier_relaxed_optimality_tol  = settings.tolerances.relative_dual_tolerance;
+  barrier_settings.barrier_eliminate_free_variables = effective_eliminate_free_variables(settings);
+  barrier_settings.barrier_initial_point_safeguard  = settings.barrier_initial_point_safeguard;
+  barrier_settings.barrier                          = true;
+  barrier_settings.barrier_presolve                 = true;
+  barrier_settings.crossover                        = settings.crossover;
+  barrier_settings.eliminate_dense_columns          = settings.eliminate_dense_columns;
+  barrier_settings.barrier_iterative_refinement     = settings.barrier_iterative_refinement;
+  barrier_settings.barrier_adaptive_regularization  = settings.barrier_adaptive_regularization;
+  barrier_settings.barrier_primal_regularization    = settings.barrier_primal_regularization;
+  barrier_settings.barrier_dual_regularization      = settings.barrier_dual_regularization;
+  barrier_settings.barrier_soc_threshold            = settings.barrier_soc_threshold;
+  barrier_settings.barrier_step_scale               = settings.barrier_step_scale;
+  barrier_settings.qcqp_ruiz_equilibration          = settings.qcqp_ruiz_equilibration;
+  barrier_settings.gpu_ruiz_nnz_threshold           = settings.gpu_ruiz_nnz_threshold;
+  barrier_settings.cudss_deterministic              = settings.cudss_deterministic;
+  barrier_settings.barrier_relaxed_feasibility_tol  = settings.tolerances.relative_primal_tolerance;
+  barrier_settings.barrier_relaxed_optimality_tol   = settings.tolerances.relative_dual_tolerance;
   barrier_settings.barrier_relaxed_complementarity_tol = settings.tolerances.relative_gap_tolerance;
   barrier_settings.barrier_relaxed_relative_objective_gap_tol =
     settings.tolerances.relative_gap_tolerance;
@@ -607,7 +623,7 @@ optimization_problem_solution_t<i_t, f_t> run_barrier(
   mip::problem_t<i_t, f_t>& problem,
   pdlp_solver_settings_t<i_t, f_t> const& settings,
   const timer_t& timer,
-  cuopt::mathematical_optimization::barrier_cache_t* cache = nullptr)
+  cuopt::mathematical_optimization::barrier_cache_t<i_t, f_t>* cache = nullptr)
 {
   // Convert data structures to dual simplex format and back
   simplex::user_problem_t<i_t, f_t> dual_simplex_problem =
@@ -1994,14 +2010,23 @@ optimization_problem_solution_t<i_t, f_t> solve_qcqp(
 
     auto* cache    = settings.barrier_cache;
     auto const* xf = (cache != nullptr && cache->dirty()) ? cache->transform() : nullptr;
-    const bool reuse_from_cache = settings.user_problem_file.empty() &&
-                                  !op_problem.has_quadratic_constraints() &&
-                                  can_reuse_barrier_cache(xf,
-                                                          effective_bound_free_variables(settings),
-                                                          op_problem.get_n_variables(),
-                                                          op_problem.get_n_constraints(),
-                                                          op_problem.has_quadratic_objective(),
-                                                          false);
+    // Same decision as the reuse check in solve_linear_program_with_barrier. A hit copies the
+    // cached problem and skips conversion. Cone sizes are compared here, before expansion.
+    // can_reuse_barrier_cache skips that comparison: the cache stores the expanded size.
+    const bool user_has_soc = op_problem.has_quadratic_constraints();
+    const bool reuse_from_cache =
+      settings.user_problem_file.empty() &&
+      can_reuse_barrier_cache(xf,
+                              effective_bound_free_variables(settings),
+                              op_problem.get_n_variables(),
+                              op_problem.get_n_constraints(),
+                              op_problem.has_quadratic_objective(),
+                              user_has_soc) &&
+      (!user_has_soc || (static_cast<i_t>(op_problem.get_quadratic_constraints().size()) ==
+                           xf->num_quadratic_constraints &&
+                         static_cast<i_t>(xf->row_sense.size()) == xf->user_num_rows &&
+                         op_problem.get_n_variables() == problem_num_cols(*xf) &&
+                         op_problem.get_n_constraints() == problem_num_rows(*xf)));
 
     if (problem_checking && !reuse_from_cache) {
       problem_checking_t<i_t, f_t>::check_problem_representation(op_problem);
@@ -2049,6 +2074,10 @@ optimization_problem_solution_t<i_t, f_t> solve_qcqp(
     // space. Reuse keeps the sense of the workspace (and Q) built by that full solve.
     if (!reuse_from_cache && cache != nullptr && cache->transform() != nullptr) {
       cache->transform()->maximize = op_problem.get_sense();
+      // Reuse never re-runs the cone expansion, so the gate above rejects a model that has
+      // gained or lost a quadratic constraint since the cache was built.
+      cache->transform()->num_quadratic_constraints =
+        static_cast<i_t>(op_problem.get_quadratic_constraints().size());
     }
     auto solution = convert_dual_simplex_sol(op_problem,
                                              std::get<0>(sol_dual_simplex),
@@ -2224,7 +2253,7 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(
         CUOPT_LOG_INFO("Skipping presolve for small problem (nnz=%d < %d)",
                        op_problem.get_nnz(),
                        presolve_nnz_threshold);
-      } else {
+      } else if (!(settings.sequence_solve && settings.method == method_t::Barrier)) {
         settings.presolver = presolver_t::PSLP;
         CUOPT_LOG_INFO("Using PSLP presolver");
       }
@@ -2234,6 +2263,8 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(
     std::unique_ptr<mip::third_party_presolve_t<i_t, f_t>> presolver;
     auto run_presolve = settings.presolver != presolver_t::None;
     run_presolve = run_presolve && settings.get_pdlp_warm_start_data().total_pdlp_iterations_ == -1;
+    // The barrier cache stores the user's LP. PSLP would reduce a different problem each update.
+    if (settings.sequence_solve && settings.method == method_t::Barrier) { run_presolve = false; }
 
     // Declare result at outer scope so that result.reduced_problem (which may be
     // referenced by problem.original_problem_ptr) remains alive through the solve.
diff --git a/cpp/src/pdlp/utilities/cython_solve.cu b/cpp/src/pdlp/utilities/cython_solve.cu
index 48646d089c..e21ab027fd 100644
--- a/cpp/src/pdlp/utilities/cython_solve.cu
+++ b/cpp/src/pdlp/utilities/cython_solve.cu
@@ -42,7 +42,7 @@
 namespace cuopt {
 namespace cython {
 
-using mathematical_optimization::barrier_cache_t;
+using barrier_cache_t = mathematical_optimization::barrier_cache_t<int, double>;
 
 /**
  * @brief Wrapper for linear_programming to expose the API to cython
diff --git a/python/cuopt/cuopt/linear_programming/data_model/data_model.py b/python/cuopt/cuopt/linear_programming/data_model/data_model.py
index 8e256dcb5b..3be764b53c 100644
--- a/python/cuopt/cuopt/linear_programming/data_model/data_model.py
+++ b/python/cuopt/cuopt/linear_programming/data_model/data_model.py
@@ -230,12 +230,11 @@ def set_objective_coefficients(self, c):
     @catch_cuopt_exception
     def update_linear_objective(self, coefficients):
         """
-        Cache reuse is QP-only: quadratic constraints take a full solve.
-
         Update the linear objective coefficients for a sequence re-solve.
         Writes ``coefficients`` onto this DataModel. If a barrier cache is
         present, also maps them into the presolved space and marks
-        it dirty (quadratic ``Q``, ``A``, and bounds must stay unchanged).
+        it dirty (quadratic ``Q``, ``A``, bounds, and the quadratic
+        constraints must stay unchanged).
 
         Parameters
         ----------
@@ -252,8 +251,8 @@ def update_rhs(self, b):
 
         Writes ``b`` onto this DataModel. If a barrier cache is present, also
         maps ``b`` into the presolved space and marks it dirty
-        (quadratic ``Q``, ``A``, row senses, and bounds must stay unchanged).
-        Cache reuse is QP-only: quadratic constraints take a full solve.
+        (quadratic ``Q``, ``A``, row senses, bounds, and the quadratic
+        constraints must stay unchanged).
 
         Range rows and folding in the first solve are not supported and raise;
         run a full solve for those models. Rows that presolve dropped as empty
diff --git a/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx
index 95c3453e57..19680ea8cf 100644
--- a/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx
+++ b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx
@@ -26,9 +26,9 @@ cdef extern from "Python.h":
     void* PyCapsule_GetPointer(object cap, const char* name)
 
 cdef extern from "cuopt/mathematical_optimization/utilities/barrier_cache.hpp" namespace "cuopt::mathematical_optimization":  # noqa
-    cdef cppclass barrier_cache_t:
-        void update_linear_objective(const double* c, int n) except +
-        void update_rhs(const double* b, int m) except +
+    cdef cppclass barrier_cache_t[i_t, f_t]:
+        void update_linear_objective(const f_t* c, i_t n) except +
+        void update_rhs(const f_t* b, i_t m) except +
 
 
 def type_cast(np_obj, np_type, name):
@@ -177,7 +177,7 @@ cdef class DataModel:
         so a later reuse can skip convert/presolve. Crush runs first so a
         length error leaves the DataModel coefficients unchanged.
         """
-        cdef barrier_cache_t* cache
+        cdef barrier_cache_t[int, double]* cache
         cdef double[::1] c_view
         new_c = type_cast(coefficients, np.float64, "coefficients")
         if self.barrier_cache_capsule is not None:
@@ -185,13 +185,13 @@ cdef class DataModel:
                 self.barrier_cache_capsule, b"cuopt.barrier_cache"
             ):
                 raise ValueError("Invalid barrier cache stored on DataModel.")
-            cache = <barrier_cache_t*>PyCapsule_GetPointer(
+            cache = <barrier_cache_t[int, double]*>PyCapsule_GetPointer(
                 self.barrier_cache_capsule,
                 b"cuopt.barrier_cache",
             )
             c_view = np.ascontiguousarray(new_c, dtype=np.float64)
             if c_view.shape[0] == 0:
-                cache.update_linear_objective(NULL, 0)
+                cache.update_linear_objective(<const double*>NULL, 0)
             else:
                 cache.update_linear_objective(&c_view[0], <int>c_view.shape[0])
         self.c = new_c
@@ -203,7 +203,7 @@ cdef class DataModel:
         the barrier cache when this model owns one. Crush runs first so a
         length error leaves the DataModel RHS unchanged.
         """
-        cdef barrier_cache_t* cache
+        cdef barrier_cache_t[int, double]* cache
         cdef double[::1] b_view
         new_b = type_cast(b, np.float64, "b")
         if self.barrier_cache_capsule is not None:
@@ -211,13 +211,13 @@ cdef class DataModel:
                 self.barrier_cache_capsule, b"cuopt.barrier_cache"
             ):
                 raise ValueError("Invalid barrier cache stored on DataModel.")
-            cache = <barrier_cache_t*>PyCapsule_GetPointer(
+            cache = <barrier_cache_t[int, double]*>PyCapsule_GetPointer(
                 self.barrier_cache_capsule,
                 b"cuopt.barrier_cache",
             )
             b_view = np.ascontiguousarray(new_b, dtype=np.float64)
             if b_view.shape[0] == 0:
-                cache.update_rhs(NULL, 0)
+                cache.update_rhs(<const double*>NULL, 0)
             else:
                 cache.update_rhs(&b_view[0], <int>b_view.shape[0])
         self.b = new_b
diff --git a/python/cuopt/cuopt/linear_programming/solver/solver.pxd b/python/cuopt/cuopt/linear_programming/solver/solver.pxd
index 84db2565de..0dd6cb17c8 100644
--- a/python/cuopt/cuopt/linear_programming/solver/solver.pxd
+++ b/python/cuopt/cuopt/linear_programming/solver/solver.pxd
@@ -95,7 +95,7 @@ cdef extern from "cuopt/mathematical_optimization/utilities/cython_types.hpp" na
         vector[double] last_restart_duality_gap_dual_solution_
 
 cdef extern from "cuopt/mathematical_optimization/utilities/barrier_cache.hpp" namespace "cuopt::mathematical_optimization": # noqa
-    cdef cppclass barrier_cache_t:
+    cdef cppclass barrier_cache_t[i_t, f_t]:
         pass
 
 cdef extern from "cuopt/mathematical_optimization/utilities/cython_solve.hpp" namespace "cuopt::cython": # noqa
@@ -122,7 +122,7 @@ cdef extern from "cuopt/mathematical_optimization/utilities/cython_solve.hpp" na
         int nb_iterations_
         double solve_time_
         method_t solved_by_
-        barrier_cache_t* barrier_cache
+        barrier_cache_t[int, double]* barrier_cache
         bool is_gpu()
 
     # Unified MIP solution struct — solution_ variant accessed via helpers
@@ -152,7 +152,7 @@ cdef extern from "cuopt/mathematical_optimization/utilities/cython_solve.hpp" na
         solver_settings_t[int, double]* solver_settings,
         unsigned int flags,
         bool is_batch_mode,
-        barrier_cache_t* cache_in,
+        barrier_cache_t[int, double]* cache_in,
     ) except + nogil
 
     cdef pair[vector[unique_ptr[solver_ret_t]], double] call_batch_solve( # noqa
diff --git a/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx
index d74a66d35a..b5bc42a543 100644
--- a/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx
+++ b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx
@@ -93,7 +93,7 @@ cdef extern from "cuopt/mathematical_optimization/utilities/barrier_cache.hpp":
     {
       void *p = PyCapsule_GetPointer(cap, "cuopt.barrier_cache");
       if (p != nullptr) {
-        delete reinterpret_cast<cuopt::mathematical_optimization::barrier_cache_t *>(p);
+        delete reinterpret_cast<cuopt::mathematical_optimization::barrier_cache_t<int, double> *>(p);
       }
     }
     """
@@ -551,8 +551,8 @@ def prepare_solver_settings(SolverSettings settings, data_model=None, mip=False)
 def Solve(py_data_model_obj, SolverSettings settings, mip=False):
 
     cdef DataModel data_model_obj = <DataModel>py_data_model_obj
-    cdef barrier_cache_t* cache_in = NULL
-    cdef barrier_cache_t* cache_out = NULL
+    cdef barrier_cache_t[int, double]* cache_in = NULL
+    cdef barrier_cache_t[int, double]* cache_out = NULL
     cdef solver_ret_t* sol_ret
 
     if (
@@ -564,7 +564,7 @@ def Solve(py_data_model_obj, SolverSettings settings, mip=False):
             b"cuopt.barrier_cache",
         ):
             raise ValueError("Invalid barrier cache stored on DataModel.")
-        cache_in = <barrier_cache_t*>PyCapsule_GetPointer(
+        cache_in = <barrier_cache_t[int, double]*>PyCapsule_GetPointer(
             data_model_obj.barrier_cache_capsule,
             b"cuopt.barrier_cache",
         )
diff --git a/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py b/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py
index 5d705e054c..bad2f1c4b8 100644
--- a/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py
+++ b/python/cuopt/cuopt/tests/linear_programming/test_barrier_sequence_solve.py
@@ -1,7 +1,7 @@
 # SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
 # SPDX-License-Identifier: Apache-2.0
 
-"""Barrier cache reuse (``sequence_solve``) for ``DataModel.update_rhs``.
+"""Barrier cache reuse (``sequence_solve``) for the DataModel update APIs.
 
 Every re-solve through the cache is compared against a fresh full solve of the
 same model, so no assertion depends on a hand-derived optimum. The models are
@@ -25,7 +25,7 @@
 
 REUSE_LOG = "reusing cache"
 RHS_INFEASIBLE_LOG = "update_rhs made an empty constraint row infeasible"
-IPM_LOG = "Optimal solution found"
+BARRIER_LOG = "Optimal solution found"
 
 
 def _sequence_settings():
@@ -37,9 +37,18 @@ def _sequence_settings():
 
 
 def _build(
-    values, indices, offsets, rhs, senses, lower, upper, objective=None
+    values,
+    indices,
+    offsets,
+    rhs,
+    senses,
+    lower,
+    upper,
+    objective=None,
+    *,
+    quadratic=True,
 ):
-    """A QP with quadratic term x^T x, so the model is barrier-eligible."""
+    """Build a model. With ``quadratic``, add ``x^T x`` so it is barrier-eligible."""
     n = len(lower)
     model = data_model.DataModel()
     model.set_csr_constraint_matrix(
@@ -54,11 +63,12 @@ def _build(
         if objective is None
         else np.asarray(objective, dtype=np.float64)
     )
-    model.set_quadratic_objective_matrix(
-        np.ones(n),
-        np.arange(n, dtype=np.int32),
-        np.arange(n + 1, dtype=np.int32),
-    )
+    if quadratic:
+        model.set_quadratic_objective_matrix(
+            np.ones(n),
+            np.arange(n, dtype=np.int32),
+            np.arange(n + 1, dtype=np.int32),
+        )
     model.set_variable_lower_bounds(np.asarray(lower, dtype=np.float64))
     model.set_variable_upper_bounds(np.asarray(upper, dtype=np.float64))
     return model
@@ -183,7 +193,7 @@ def test_update_rhs_infeasible_empty_row_short_circuits(capfd):
 
     The row has no variables, so it is either satisfied for every x or for
     none. That makes the verdict exact and lets the next Solve answer without
-    running IPM. The cache is kept, so a later feasible RHS still reuses it.
+    running barrier. The cache is kept, so a later feasible RHS still reuses it.
     """
     settings = _sequence_settings()
     model = _build(**dict(EMPTY_ROW, rhs=[0.0, 2.0]))
@@ -195,7 +205,9 @@ def test_update_rhs_infeasible_empty_row_short_circuits(capfd):
     infeasible, log = _solve(model, settings, capfd)
     assert infeasible.get_termination_reason() == "PrimalInfeasible"
     assert RHS_INFEASIBLE_LOG in log
-    assert IPM_LOG not in log, "IPM ran despite a provably infeasible row"
+    assert BARRIER_LOG not in log, (
+        "barrier ran despite a provably infeasible row"
+    )
 
     # Same cache, feasible RHS again.
     rhs = [0.0, 4.0]
@@ -290,11 +302,271 @@ def test_bounded_free_variables_block_reuse(
     )
 
 
+# min x + 2y, x + y >= b, 0 <= x <= 2, y free. The optimum is unique: x = 2,
+# y = b - 2. y is free, so a normal LP barrier solve would split it into v - w.
+PURE_LP = dict(
+    values=[1.0, 1.0],
+    indices=[0, 1],
+    offsets=[0, 2],
+    senses="G",
+    lower=[0.0, -np.inf],
+    upper=[2.0, np.inf],
+    objective=[1.0, 2.0],
+)
+
+
+def test_lp_barrier_update_rhs_reuses_cache(capfd):
+    """A pure LP barrier sequence solve reuses the cache after an RHS update.
+
+    There is no quadratic term, so this is the LP path. PSLP, folding, and the
+    v-w split would each change the matrix the cache has to replay.
+    """
+    settings = _sequence_settings()
+    settings.set_parameter("method", solver_settings.SolverMethod.Barrier)
+    model = _build(**dict(PURE_LP, rhs=[3.0]), quadratic=False)
+
+    first, log = _solve(model, settings, capfd)
+    assert first.get_termination_reason() == "Optimal"
+    assert REUSE_LOG not in log
+    assert "Using PSLP presolver" not in log
+    assert "Handling 1 free variables directly" in log
+    assert "Folding:" not in log
+
+    rhs = [6.0]
+    model.update_rhs(np.asarray(rhs, dtype=np.float64))
+    reused, log = _solve(model, settings, capfd)
+    assert REUSE_LOG in log, (
+        "LP barrier sequence solve fell back to a full solve"
+    )
+
+    oracle_settings = solver_settings.SolverSettings()
+    oracle_settings.set_parameter(
+        "method", solver_settings.SolverMethod.Barrier
+    )
+    oracle = solver.Solve(
+        _build(**dict(PURE_LP, rhs=rhs), quadratic=False), oracle_settings
+    )
+    _assert_matches_oracle(reused, oracle)
+
+
+# Quadratic constraints reach the barrier as second-order cones. The conversion appends rows and
+# permutes columns, so an update has to be mapped from model coordinates into the expanded layout
+# the cache holds, and updates stay model-sized: the appended rows belong to the conversion.
+# The conversion also rejects cone variables that carry an explicit upper bound or a nonzero
+# lower bound, so the builders below leave the bounds of anything a cone touches open.
+
+
+def _build_lorentz(rhs, objective, cone_head_free=False):
+    """``||(x1, x2)|| <= t`` plus two linear rows, no quadratic objective.
+
+    Written as ``-t^2 + x1^2 + x2^2 <= 0``, which the conversion recognizes
+    and lifts by permuting ``(t, x1, x2)`` into a cone block, so all three are
+    conic variables. With ``cone_head_free``, ``t`` loses its lower bound and
+    row 1 becomes the singleton ``t >= b1`` that proves the head nonnegative.
+    """
+    n = 3
+    model = data_model.DataModel()
+    # row 0: x1 + x2 >= b0 ; row 1: t <= b1, or t >= b1 when the head is free
+    model.set_csr_constraint_matrix(
+        np.array([1.0, 1.0, 1.0], dtype=np.float64),
+        np.array([1, 2, 0], dtype=np.int32),
+        np.array([0, 2, 3], dtype=np.int32),
+    )
+    model.set_constraint_bounds(np.asarray(rhs, dtype=np.float64))
+    model.set_row_types("GG" if cone_head_free else "GL")
+    model.set_objective_coefficients(np.asarray(objective, dtype=np.float64))
+    lower = np.zeros(n)
+    if cone_head_free:
+        lower[0] = -np.inf
+    model.set_variable_lower_bounds(lower)
+    model.set_variable_upper_bounds(np.full(n, np.inf))
+    model.add_quadratic_constraint(
+        vals=np.array([-1.0, 1.0, 1.0]),
+        rows=np.array([0, 1, 2], dtype=np.int32),
+        cols=np.array([0, 1, 2], dtype=np.int32),
+        rhs_value=0.0,
+        sense="L",
+    )
+    return model
+
+
+def _build_general_cone(rhs, objective):
+    """``x1^2 + x2^2 + 2*x1 <= 8``, a shifted disk, plus two linear rows.
+
+    The linear part and the nonzero RHS keep this off the recognized Lorentz
+    pattern, so the conversion takes its general path: it factors Q, adds cone
+    variables of its own, and appends four rows instead of two, which is a
+    different expanded shape to map an update through. The disk bounds
+    ``min x1``, so the model stays bounded.
+    """
+    n = 2
+    model = data_model.DataModel()
+    # row 0: x1 + x2 >= b0 ; row 1: x1 - x2 <= b1
+    model.set_csr_constraint_matrix(
+        np.array([1.0, 1.0, 1.0, -1.0], dtype=np.float64),
+        np.array([0, 1, 0, 1], dtype=np.int32),
+        np.array([0, 2, 4], dtype=np.int32),
+    )
+    model.set_constraint_bounds(np.asarray(rhs, dtype=np.float64))
+    model.set_row_types("GL")
+    model.set_objective_coefficients(np.asarray(objective, dtype=np.float64))
+    model.set_variable_lower_bounds(np.full(n, -np.inf))
+    model.set_variable_upper_bounds(np.full(n, np.inf))
+    model.add_quadratic_constraint(
+        vals=np.array([1.0, 1.0]),
+        rows=np.array([0, 1], dtype=np.int32),
+        cols=np.array([0, 1], dtype=np.int32),
+        linear_values=np.array([2.0]),
+        linear_indices=np.array([0], dtype=np.int32),
+        rhs_value=8.0,
+        sense="L",
+    )
+    return model
+
+
+def _full_cone_solve(build, rhs, objective, **kwargs):
+    """Oracle: a fresh cone model with default settings, no cache in play."""
+    return solver.Solve(
+        build(rhs, objective, **kwargs), solver_settings.SolverSettings()
+    )
+
+
+# Builder, objective, first RHS, then the RHSs to update to. The updates stay
+# feasible: on the shifted disk, x1 + x2 tops out just above 2.16.
+CONE_CASES = {
+    "lorentz": (
+        _build_lorentz,
+        [1.0, 0.0, 0.0],
+        [2.0, 9.0],
+        [[3.0, 9.0], [1.5, 9.0]],
+    ),
+    "general": (
+        _build_general_cone,
+        [1.0, 0.0],
+        [1.0, 5.0],
+        [[1.5, 5.0], [0.5, 5.0]],
+    ),
+}
+
+
+@pytest.mark.parametrize("case", list(CONE_CASES))
+def test_cone_update_rhs_matches_full_solve(case, capfd):
+    """An RHS update on a cone model must agree with a fresh full solve.
+
+    The model RHS has two entries while the converted problem has more rows,
+    so this fails outright if the update is not mapped into the expanded
+    layout, and gives a wrong answer if the appended rows lose their own RHS.
+    """
+    build, objective, first_rhs, later_rhs = CONE_CASES[case]
+    settings = _sequence_settings()
+    model = build(first_rhs, objective)
+
+    first, _ = _solve(model, settings, capfd)
+    assert first.get_termination_reason() == "Optimal"
+
+    for rhs in later_rhs:
+        model.update_rhs(np.asarray(rhs, dtype=np.float64))
+        reused, log = _solve(model, settings, capfd)
+        assert REUSE_LOG in log, "update_rhs fell back to a full solve"
+        _assert_matches_oracle(reused, _full_cone_solve(build, rhs, objective))
+
+
+@pytest.mark.parametrize("case", list(CONE_CASES))
+def test_cone_update_rhs_moves_the_optimum(case):
+    """Guard the oracles: the RHS being updated has to matter.
+
+    Without this, a cached solve that ignored the new RHS would still match an
+    oracle that ignored it too, and every comparison above would pass.
+    """
+    build, objective, first_rhs, later_rhs = CONE_CASES[case]
+    objectives = {
+        _full_cone_solve(build, rhs, objective).get_primal_objective()
+        for rhs in [first_rhs] + later_rhs
+    }
+    assert len(objectives) == 1 + len(later_rhs)
+
+
+def test_lorentz_optimum_is_the_expected_one():
+    """Pin one closed form, so the oracles are not the only reference.
+
+    ``||(x1, x2)|| <= t`` with ``x1 + x2 >= b`` and ``x >= 0`` splits the sum
+    evenly, putting the optimum of ``min t`` at ``b / sqrt(2)``.
+    """
+    for b in (2.0, 3.0):
+        solution = _full_cone_solve(_build_lorentz, [b, 9.0], [1.0, 0.0, 0.0])
+        assert solution.get_termination_reason() == "Optimal"
+        assert solution.get_primal_objective() == pytest.approx(
+            b / np.sqrt(2.0), rel=1e-5, abs=1e-5
+        )
+
+
+def test_cone_update_linear_objective_matches_full_solve(capfd):
+    """The conversion permutes columns, so the objective needs remapping too."""
+    settings = _sequence_settings()
+    rhs = [2.0, 9.0]
+    model = _build_lorentz(rhs, [1.0, 0.0, 0.0])
+
+    first, _ = _solve(model, settings, capfd)
+    assert first.get_termination_reason() == "Optimal"
+
+    new_objective = [1.0, 0.25, 0.0]
+    model.update_linear_objective(np.asarray(new_objective, dtype=np.float64))
+    reused, log = _solve(model, settings, capfd)
+    assert REUSE_LOG in log, (
+        "update_linear_objective fell back to a full solve"
+    )
+    _assert_matches_oracle(
+        reused, _full_cone_solve(_build_lorentz, rhs, new_objective)
+    )
+
+
+def test_cone_update_rhs_rejects_lost_cone_head_bound(capfd):
+    """A cone head proved nonnegative by a row cannot lose that proof.
+
+    The head here is free, so the conversion only accepts the model because
+    row 1 forces t >= 0. An RHS that relaxes the row to t >= -1 makes this a
+    model a full solve refuses, and reuse has to refuse it the same way rather
+    than solving a stale cone formulation.
+    """
+    settings = _sequence_settings()
+    objective = [1.0, 0.0, 0.0]
+    model = _build_lorentz([2.0, 0.0], objective, cone_head_free=True)
+
+    first, _ = _solve(model, settings, capfd)
+    assert first.get_termination_reason() == "Optimal"
+
+    with pytest.raises(Exception, match="nonnegative"):
+        model.update_rhs(np.array([2.0, -1.0]))
+
+    # A full solve of the same model does not return an optimum either,
+    # whether it reports the rejection as an exception or a status.
+    try:
+        oracle = _full_cone_solve(
+            _build_lorentz, [2.0, -1.0], objective, cone_head_free=True
+        )
+    except Exception:
+        pass
+    else:
+        assert oracle.get_termination_reason() != "Optimal"
+
+
+def test_cone_update_rhs_rejects_wrong_length(capfd):
+    """Length is validated against the problem rows, not the converted rows."""
+    settings = _sequence_settings()
+    model = _build_lorentz([2.0, 9.0], [1.0, 0.0, 0.0])
+    first, _ = _solve(model, settings, capfd)
+    assert first.get_termination_reason() == "Optimal"
+
+    # Two model rows. Passing the converted row count must not be accepted.
+    with pytest.raises(Exception, match="match the cached problem row count"):
+        model.update_rhs(np.array([2.0, 9.0, 0.0, 0.0]))
+
+
 def test_update_rhs_rejects_wrong_length():
-    """Length is validated against the cached user row count."""
+    """Length is validated against the cached problem row count."""
     settings = _sequence_settings()
     model = _build(**dict(MIXED_SENSES, rhs=[5.0, 8.0, 3.0]))
     assert solver.Solve(model, settings).get_termination_reason() == "Optimal"
 
-    with pytest.raises(Exception, match="match the cached user row count"):
+    with pytest.raises(Exception, match="match the cached problem row count"):
         model.update_rhs(np.array([1.0, 2.0]))