diff --git a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp
index bd07aaac32..d97231c9f6 100644
--- a/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp
+++ b/cpp/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp
@@ -25,6 +25,8 @@
 namespace cuopt {
 namespace CUOPT_EXPORT mathematical_optimization {
 
+class barrier_cache_t;
+
 // Forward declare solver_settings_t for friend class
 template <typename i_t, typename f_t>
 class solver_settings_t;
@@ -365,6 +367,10 @@ class pdlp_solver_settings_t {
   // Used to force batch PDLP to solve a subbatch of the problems at a time
   // The 0 default value will make the solver use its heuristic to determine the subbatch size
   i_t fixed_batch_size{0};
+  /** When true, first GPU barrier/QCQP solve returns a ``barrier_cache_t`` capsule. */
+  bool sequence_solve{false};
+  /** Non-owning cache pointer set by ``call_solve`` for barrier symbolic reuse. */
+  barrier_cache_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
new file mode 100644
index 0000000000..cae544cb1a
--- /dev/null
+++ b/cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp
@@ -0,0 +1,84 @@
+/* clang-format off */
+/*
+ * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-License-Identifier: Apache-2.0
+ */
+/* clang-format on */
+
+#pragma once
+
+#include <cuopt/export.hpp>
+
+#include <memory>
+
+#include <raft/core/handle.hpp>
+#include <rmm/cuda_stream.hpp>
+
+namespace cuopt::mathematical_optimization::barrier {
+template <typename i_t, typename f_t>
+class iteration_data_t;
+
+void destroy_iteration_data(iteration_data_t<int, double>* data);
+
+void apply_barrier_linear_objective(iteration_data_t<int, double>& data,
+                                    double const* barrier_c,
+                                    int n);
+}  // namespace cuopt::mathematical_optimization::barrier
+
+namespace cuopt {
+namespace CUOPT_EXPORT mathematical_optimization {
+
+struct barrier_transform_t;
+
+/**
+ * @brief GPU solve cache owned by DataModel when sequence_solve is on.
+ *
+ * After an Optimal full solve, holds iteration_data_t and the user-barrier transform.
+ * update_linear_objective crushes the new linear objective and sets c_dirty so the next Solve
+ * reuses that workspace (skip convert/presolve/scaling).
+ */
+class barrier_cache_t {
+ public:
+  static std::unique_ptr<barrier_cache_t> create(unsigned stream_flags);
+
+  barrier_cache_t(barrier_cache_t&&) noexcept;
+  barrier_cache_t& operator=(barrier_cache_t&&) noexcept;
+  ~barrier_cache_t();
+
+  [[nodiscard]] raft::handle_t* handle_ptr();
+  [[nodiscard]] raft::handle_t const* handle_ptr() const;
+
+  /** Drop cached iteration workspace and transform (handle/stream stay). */
+  void clear();
+
+  /**
+   * @brief Take ownership of barrier iteration workspace. @p data may be null (clears).
+   */
+  void store_iteration_data(barrier::iteration_data_t<int, double>* data);
+
+  /**
+   * @brief Release ownership of cached iteration workspace; caller must delete or wrap it.
+   */
+  barrier::iteration_data_t<int, double>* 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 set_c_dirty(bool dirty);
+  [[nodiscard]] bool c_dirty() const;
+
+  /**
+   * Crush the input linear objective into cached iteration_data_t.c / d_c_ and set c_dirty.
+   * Requires a stored transform and iteration_data from an Optimal solve.
+   */
+  void update_linear_objective(double const* c, int n);
+
+ private:
+  barrier_cache_t(std::unique_ptr<rmm::cuda_stream> stream, std::unique_ptr<raft::handle_t> handle);
+
+  struct impl;
+  std::unique_ptr<impl> impl_;
+};
+
+}  // namespace CUOPT_EXPORT mathematical_optimization
+}  // namespace cuopt
diff --git a/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp b/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp
index f84119a8dc..0cee4e34f3 100644
--- a/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp
+++ b/cpp/include/cuopt/mathematical_optimization/utilities/cython_solve.hpp
@@ -8,12 +8,11 @@
 #pragma once
 
 #include <cuopt/export.hpp>
+#include <cuopt/mathematical_optimization/io/data_model_view.hpp>
 #include <cuopt/mathematical_optimization/optimization_problem_interface.hpp>
 #include <cuopt/mathematical_optimization/optimization_problem_solution_interface.hpp>
 #include <cuopt/mathematical_optimization/solver_settings.hpp>
 #include <cuopt/mathematical_optimization/utilities/cython_types.hpp>
-
-#include <cuopt/mathematical_optimization/io/data_model_view.hpp>
 #include <memory>
 #include <raft/core/handle.hpp>
 #include <string>
@@ -55,8 +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);
+  unsigned int flags                                   = cudaStreamNonBlocking,
+  bool is_batch_mode                                   = false,
+  mathematical_optimization::barrier_cache_t* 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 69d6f91604..e1f9c423f7 100644
--- a/cpp/include/cuopt/mathematical_optimization/utilities/cython_types.hpp
+++ b/cpp/include/cuopt/mathematical_optimization/utilities/cython_types.hpp
@@ -10,6 +10,7 @@
 #include <cuopt/export.hpp>
 #include <cuopt/mathematical_optimization/mip/solver_solution.hpp>
 #include <cuopt/mathematical_optimization/pdlp/solver_solution.hpp>
+#include <cuopt/mathematical_optimization/utilities/barrier_cache.hpp>
 #include <cuopt/mathematical_optimization/utilities/internals.hpp>
 
 #include <rmm/device_buffer.hpp>
@@ -86,6 +87,10 @@ struct linear_programming_ret_t {
   double solve_time_{};
   mathematical_optimization::method_t solved_by_{};
 
+  /** GPU barrier cache (stream + handle + iteration workspace); moved to Python capsule when set.
+   */
+  std::unique_ptr<mathematical_optimization::barrier_cache_t> barrier_cache;
+
   bool is_gpu() const { return std::holds_alternative<gpu_solutions_t>(solutions_); }
 };
 
diff --git a/cpp/src/barrier/CMakeLists.txt b/cpp/src/barrier/CMakeLists.txt
index 650bc733e9..9556f3992a 100644
--- a/cpp/src/barrier/CMakeLists.txt
+++ b/cpp/src/barrier/CMakeLists.txt
@@ -6,6 +6,7 @@
 set(BARRIER_SRC_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/cusparse_view.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/barrier.cu
+  ${CMAKE_CURRENT_SOURCE_DIR}/barrier_cache.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/device_sparse_matrix.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/pinned_host_allocator.cu
   )
diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu
index dca523669c..2763e142cf 100644
--- a/cpp/src/barrier/barrier.cu
+++ b/cpp/src/barrier/barrier.cu
@@ -29,6 +29,12 @@
 
 #include <linear_algebra/vector_math.cuh>
 
+#include <cuopt/error.hpp>
+#include <cuopt/mathematical_optimization/utilities/barrier_cache.hpp>
+
+#include <barrier/barrier_transform.hpp>
+
+#include <rmm/device_buffer.hpp>
 #include <rmm/device_scalar.hpp>
 #include <rmm/device_uvector.hpp>
 
@@ -469,7 +475,7 @@ class iteration_data_t {
       Q(Qin),
       cusparse_Q_view_(lp.handle_ptr, Q),
       cusparse_view_(lp.handle_ptr, lp.A),
-      cusparse_info(lp.handle_ptr),
+      cusparse_info_(nullptr),
       device_AD(lp.num_cols, lp.num_rows, 0, lp.handle_ptr->get_stream()),
       device_A(lp.num_cols, lp.num_rows, 0, lp.handle_ptr->get_stream()),
       device_ADAT(lp.num_rows, lp.num_rows, 0, lp.handle_ptr->get_stream()),
@@ -479,7 +485,6 @@ class iteration_data_t {
       device_Q_csc_(lp.handle_ptr->get_stream()),
       device_AT_csc_(lp.handle_ptr->get_stream()),
       d_original_A_values(0, lp.handle_ptr->get_stream()),
-      device_A_x_values(0, lp.handle_ptr->get_stream()),
       d_inv_diag_prime(0, lp.handle_ptr->get_stream()),
       d_flag_buffer(0, lp.handle_ptr->get_stream()),
       d_num_flag(lp.handle_ptr->get_stream()),
@@ -723,7 +728,7 @@ class iteration_data_t {
         if (n_dense_rows > 0) {
           settings.log.printf("Dense rows                  : %d\n", n_dense_rows);
         }
-        settings.log.printf("Density estimator time      : %.2fs\n", column_density_time);
+        settings.log.printf("Density estimator time      : %.3fs\n", column_density_time);
         if ((settings.augmented != 0) &&
             (n_dense_columns > 50 || n_dense_rows > 10 ||
              lp.A.m == 0 /* handle case with no constraints */ ||
@@ -843,9 +848,6 @@ class iteration_data_t {
                  handle_ptr->get_stream());
       // For efficient scaling of AD col we form the col index array
       device_AD.form_col_index(handle_ptr->get_stream());
-      device_A_x_values.resize(device_AD.x.size(), handle_ptr->get_stream());
-      raft::copy(
-        device_A_x_values.data(), device_AD.x.data(), device_AD.x.size(), handle_ptr->get_stream());
       device_AD.to_compressed_row(device_A, handle_ptr->get_stream());
       RAFT_CHECK_CUDA(handle_ptr->get_stream());
     }
@@ -856,7 +858,7 @@ class iteration_data_t {
       i_t factorization_size =
         use_augmented ? augmented_system_size(lp.num_cols, lp.num_rows) : lp.num_rows;
       chol = std::make_unique<sparse_cholesky_cudss_t<i_t, f_t>>(
-        handle_ptr, settings, factorization_size);
+        handle_ptr, settings_, factorization_size);
       chol->set_positive_definite(false);
     }
     if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; }
@@ -883,6 +885,68 @@ class iteration_data_t {
     }
   }
 
+  // Attach this solve's settings and rewind iterate-dependent state so barrier can
+  // Mehrotra-start with the new c. 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.
+  bool reset_iterate_state(const simplex_solver_settings_t<i_t, f_t>& settings)
+  {
+    if (chol == nullptr || symbolic_status != 0) { return false; }
+    settings_ = settings;
+
+    {
+      raft::common::nvtx::range fun_scope("Barrier: reset diagonal scaling");
+      const bool has_Q        = Q.n > 0;
+      const bool has_soc      = has_cones();
+      const bool adaptive_reg = should_use_adaptive_regularization(settings, has_soc);
+      primal_perturb          = (settings.barrier_primal_regularization >= 0)
+                                  ? settings.barrier_primal_regularization
+                                  : (has_soc ? 1e-8 : 1e-6);
+      dual_perturb            = (settings.barrier_dual_regularization >= 0)
+                                  ? settings.barrier_dual_regularization
+                                  : (adaptive_reg ? 1e-8 : 0);
+
+      diag.set_scalar(1.0);
+      for (i_t k = 0; k < n_upper_bounds; k++) {
+        diag[upper_bounds[k]] = 2.0;
+      }
+      if (has_Q && !use_augmented) {
+        for (i_t j = 0; j < Q.n; j++) {
+          diag[j] += Qdiag[j];
+        }
+      }
+
+      inv_diag.set_scalar(1.0);
+      if (n_upper_bounds > 0 || (has_Q && !use_augmented)) { diag.inverse(inv_diag); }
+      raft::copy(d_diag_.data(), diag.data(), diag.size(), stream_view_);
+      raft::copy(d_inv_diag.data(), inv_diag.data(), inv_diag.size(), stream_view_);
+      inv_sqrt_diag.set_scalar(1.0);
+      if (n_upper_bounds > 0 || (has_Q && !use_augmented)) { inv_diag.sqrt(inv_sqrt_diag); }
+    }
+
+    if (use_augmented) {
+      form_augmented(false);
+    } else {
+      form_adat(false);
+      handle_ptr->sync_stream();
+      if (chol != nullptr) { chol->rebind_csr_matrix(device_ADAT); }
+    }
+    if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) { return false; }
+
+    // Drop the last iterate's numeric factor and residual history; keep symbolic analysis.
+    has_factorization                      = false;
+    has_solve_info                         = false;
+    relative_primal_residual_save          = inf;
+    relative_dual_residual_save            = inf;
+    relative_complementarity_residual_save = inf;
+    primal_residual_norm_save              = inf;
+    dual_residual_norm_save                = inf;
+    complementarity_residual_norm_save     = inf;
+    if (chol != nullptr) { chol->invalidate_numeric_factor(); }
+    handle_ptr->sync_stream();
+    return true;
+  }
+
   bool has_cones() const { return cones_.has_value(); }
 
   cone_data_t<i_t, f_t>& cones()
@@ -1130,8 +1194,11 @@ class iteration_data_t {
     if (first_call) {
       raft::common::nvtx::range scope("Barrier: Form ADAT: cusparse init");
       try {
+        if (!cusparse_info_) {
+          cusparse_info_ = std::make_unique<cusparse_info_t<i_t, f_t>>(handle_ptr);
+        }
         initialize_cusparse_data<i_t, f_t>(
-          handle_ptr, device_A, device_AD, device_ADAT, cusparse_info);
+          handle_ptr, device_A, device_AD, device_ADAT, spgemm_info());
       } catch (const raft::cuda_error& e) {
         settings_.log.printf("Error in initialize_cusparse_data: %s\n", e.what());
         return;
@@ -1141,7 +1208,7 @@ class iteration_data_t {
 
     {
       raft::common::nvtx::range scope("Barrier: Form ADAT: ADAT multiply");
-      multiply_kernels<i_t, f_t>(handle_ptr, device_A, device_AD, device_ADAT, cusparse_info);
+      multiply_kernels<i_t, f_t>(handle_ptr, device_A, device_AD, device_ADAT, spgemm_info());
       handle_ptr->sync_stream();
     }
 
@@ -1149,7 +1216,7 @@ class iteration_data_t {
     float64_t adat_time = toc(start_form_adat);
 
     if (num_factorizations == 0) {
-      settings_.log.printf("ADAT time                   : %.2fs\n", adat_time);
+      settings_.log.printf("ADAT time                   : %.3fs\n", adat_time);
       settings_.log.printf("ADAT nonzeros               : %.2e\n",
                            static_cast<float64_t>(adat_nnz));
       settings_.log.printf(
@@ -2148,7 +2215,6 @@ class iteration_data_t {
   device_csr_matrix_t<i_t, f_t> device_ADAT;
   device_csr_matrix_t<i_t, f_t> device_A;
   device_csc_matrix_t<i_t, f_t> device_AD;
-  rmm::device_uvector<f_t> device_A_x_values;
   // For GPU Form ADAT
   rmm::device_uvector<f_t> d_inv_diag_prime;
   rmm::device_buffer d_flag_buffer;
@@ -2171,6 +2237,13 @@ class iteration_data_t {
   std::vector<f_t> Qdiag;
   bool Q_diagonal;
   rmm::device_uvector<i_t> d_augmented_diagonal_indices_;
+
+  cusparse_info_t<i_t, f_t>& spgemm_info()
+  {
+    cuopt_assert(cusparse_info_ != nullptr, "spgemm_info: cusparse workspace unset");
+    return *cusparse_info_;
+  }
+
   cone_kkt_data_t<i_t, f_t> cone_kkt_data_;
   bool indefinite_Q;
   cusparse_view_t<i_t, f_t> cusparse_Q_view_;
@@ -2194,7 +2267,7 @@ class iteration_data_t {
   bool has_solve_info;
   i_t num_factorizations;
 
-  cusparse_info_t<i_t, f_t> cusparse_info;
+  std::unique_ptr<cusparse_info_t<i_t, f_t>> cusparse_info_;
   cusparse_view_t<i_t, f_t> cusparse_view_;
   pdlp::cusparse_dn_vec_descr_wrapper_t<f_t> cusparse_tmp4_;
   pdlp::cusparse_dn_vec_descr_wrapper_t<f_t> cusparse_h_;
@@ -2288,7 +2361,7 @@ class iteration_data_t {
 
   rmm::cuda_stream_view stream_view_;
 
-  const simplex_solver_settings_t<i_t, f_t>& settings_;
+  simplex_solver_settings_t<i_t, f_t> settings_;
 };
 
 // Move the Cholesky debug logic to a reusable function.
@@ -4201,7 +4274,7 @@ lp_status_t barrier_solver_t<i_t, f_t>::check_for_suboptimal_solution(
                      solution);
     settings.log.printf("\n");
     settings.log.printf(
-      "Suboptimal solution found in %d iterations and %.2f seconds\n", iter, toc(start_time));
+      "Suboptimal solution found in %d iterations and %.3f seconds\n", iter, toc(start_time));
     settings.log.printf("Objective %+.8e\n", compute_user_objective(lp, primal_objective));
     settings.log.printf("Primal infeasibility (abs/rel): %8.2e/%8.2e\n",
                         primal_residual_norm,
@@ -4240,7 +4313,7 @@ lp_status_t barrier_solver_t<i_t, f_t>::check_for_suboptimal_solution(
                      solution);
     settings.log.printf("\n");
     settings.log.printf(
-      "Suboptimal solution found in %d iterations and %.2f seconds\n", iter, toc(start_time));
+      "Suboptimal solution found in %d iterations and %.3f seconds\n", iter, toc(start_time));
     settings.log.printf("Objective %+.8e\n", compute_user_objective(lp, primal_objective_save));
     settings.log.printf("Primal infeasibility (abs/rel): %8.2e/%8.2e\n",
                         data.primal_residual_norm_save,
@@ -4265,56 +4338,11 @@ lp_status_t barrier_solver_t<i_t, f_t>::check_for_suboptimal_solution(
 }
 
 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)
+lp_status_t barrier_solver_t<i_t, f_t>::barrier_advanced_solve(f_t start_time,
+                                                               lp_solution_t<i_t, f_t>& solution,
+                                                               iteration_data_t<i_t, f_t>& data)
 {
-  settings.log.printf("Barrier solver started at %.2f seconds\n", toc(start_time));
-  try {
-    raft::common::nvtx::range fun_scope("Barrier: solve");
-
-    i_t n = lp.num_cols;
-    i_t m = lp.num_rows;
-
-    solution.resize(m, n);
-    settings.log.printf(
-      "Barrier solver: %d constraints, %d variables, %ld nonzeros\n", m, n, lp.A.col_start[n]);
-
-    settings.log.printf("\n");
-
-    if (lp.Q.n > 0) {
-      settings.log.printf("Quadratic objective matrix  : %d nonzeros\n", lp.Q.row_start[lp.Q.n]);
-    }
-    if (lp.second_order_cone_dims.size() > 0) {
-      settings.log.printf("Second-order cones          : %d\n",
-                          static_cast<int>(lp.second_order_cone_dims.size()));
-    }
-
-    // Compute the number of free variables
-    i_t num_free_variables = presolve_info.free_variable_pairs.size() / 2;
-    if (num_free_variables > 0) {
-      settings.log.printf("Free variables              : %d\n", num_free_variables);
-    }
-
-    // Compute the number of upper bounds
-    i_t num_upper_bounds = 0;
-    for (i_t j = 0; j < n; j++) {
-      if (lp.upper[j] < inf) { num_upper_bounds++; }
-    }
-
-    csc_matrix_t<i_t, f_t> Q(lp.num_cols, 0, 0);
-    if (lp.Q.n > 0) { create_Q(lp, Q); }
-
-    iteration_data_t<i_t, f_t> data(
-      lp, num_upper_bounds, presolve_info.direct_free_variables, Q, settings);
-    if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) {
-      settings.log.printf("Barrier solver halted\n");
-      return lp_status_t::CONCURRENT_LIMIT;
-    }
-    if (data.indefinite_Q) { return lp_status_t::NUMERICAL_ISSUES; }
-    if (data.symbolic_status != 0) {
-      settings.log.printf("Error in symbolic analysis\n");
-      return lp_status_t::NUMERICAL_ISSUES;
-    }
-
+  {
     data.cusparse_dual_residual_ = data.cusparse_view_.create_vector(data.d_dual_residual_);
     data.cusparse_r1_            = data.cusparse_view_.create_vector(data.d_r1_);
     data.cusparse_tmp4_          = data.cusparse_view_.create_vector(data.d_tmp4_);
@@ -4322,9 +4350,9 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
     data.cusparse_dx_residual_   = data.cusparse_view_.create_vector(data.d_dx_residual_);
     data.cusparse_u_             = data.cusparse_view_.create_vector(data.d_u_);
     data.cusparse_y_residual_    = data.cusparse_view_.create_vector(data.d_y_residual_);
-    data.restrict_u_.resize(num_upper_bounds);
+    data.restrict_u_.resize(data.n_upper_bounds);
 
-    settings.log.printf("Elapsed time                : %.2fs\n", toc(start_time));
+    settings.log.printf("Elapsed time                : %.3fs\n", toc(start_time));
 
     if (toc(start_time) > settings.time_limit) {
       settings.log.printf("Barrier time limit exceeded\n");
@@ -4353,7 +4381,6 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
       settings.log.printf("Unable to compute initial point\n");
       return lp_status_t::NUMERICAL_ISSUES;
     }
-
     // Upload initial point to device and compute initial residuals/norms on GPU
     data.d_complementarity_wv_residual_.resize(data.n_upper_bounds, stream_view_);
     data.d_complementarity_wv_rhs_.resize(data.n_upper_bounds, stream_view_);
@@ -4421,7 +4448,7 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
     settings.log.printf(
       "Iter   Primal              Dual                Primal   Dual    Compl.   Elapsed\n");
     float64_t elapsed_time = toc(start_time);
-    settings.log.printf("%3d   %+.12e %+.12e %.2e %.2e %.2e %.1f\n",
+    settings.log.printf("%3d   %+.12e %+.12e %.2e %.2e %.2e %.3f\n",
                         iter,
                         compute_user_objective(lp, primal_objective),
                         compute_user_objective(lp, dual_objective),
@@ -4615,7 +4642,7 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
                                              solution);
       }
 
-      settings.log.printf("%3d   %+.12e %+.12e %.2e %.2e %.2e %.1f\n",
+      settings.log.printf("%3d   %+.12e %+.12e %.2e %.2e %.2e %.3f\n",
                           iter,
                           compute_user_objective(lp, primal_objective),
                           compute_user_objective(lp, dual_objective),
@@ -4701,6 +4728,153 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
                      data.cusparse_view_,
                      solution);
     return lp_status_t::ITERATION_LIMIT;
+  }
+}
+
+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)
+{
+  settings.log.printf("Barrier solver started at %.3f seconds\n", toc(start_time));
+  try {
+    raft::common::nvtx::range fun_scope("Barrier: barrier_advanced_solve");
+
+    i_t n = lp.num_cols;
+    i_t m = lp.num_rows;
+    solution.resize(m, n);
+    settings.log.printf(
+      "Barrier solver: %d constraints, %d variables, %ld nonzeros\n", m, n, lp.A.col_start[n]);
+    if (lp.Q.n > 0) {
+      settings.log.printf("Quadratic objective matrix  : %d nonzeros\n", lp.Q.row_start[lp.Q.n]);
+    }
+
+    std::unique_ptr<iteration_data_t<i_t, f_t>> owned_data;
+    if (cache != nullptr) {
+      if (auto* cached = cache->release_iteration_data()) { owned_data.reset(cached); }
+    }
+    if (!owned_data) {
+      if (cache != nullptr) { cache->clear(); }
+      settings.log.printf(
+        "Barrier: cache reuse failed; cached iteration_data is missing or invalid\n");
+      return lp_status_t::NUMERICAL_ISSUES;
+    }
+    try {
+      if (!owned_data->reset_iterate_state(settings)) {
+        owned_data.reset();
+        if (cache != nullptr) { cache->clear(); }
+        settings.log.printf(
+          "Barrier: cache reuse failed; cached iteration_data is missing or invalid\n");
+        return lp_status_t::NUMERICAL_ISSUES;
+      }
+    } catch (const raft::cuda_error&) {
+      owned_data.reset();
+      if (cache != nullptr) { cache->clear(); }
+      settings.log.printf(
+        "Barrier: cache reuse failed; cached iteration_data is missing or invalid\n");
+      return lp_status_t::NUMERICAL_ISSUES;
+    }
+
+    iteration_data_t<i_t, f_t>& data = *owned_data;
+    if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) {
+      settings.log.printf("Barrier solver halted\n");
+      if (cache != nullptr) { cache->clear(); }
+      return lp_status_t::CONCURRENT_LIMIT;
+    }
+    if (data.indefinite_Q) {
+      if (cache != nullptr) { cache->clear(); }
+      return lp_status_t::NUMERICAL_ISSUES;
+    }
+    if (data.symbolic_status != 0) {
+      settings.log.printf("Error in symbolic analysis\n");
+      if (cache != nullptr) { cache->clear(); }
+      return lp_status_t::NUMERICAL_ISSUES;
+    }
+    settings.log.printf("Barrier setup complete at %.3f seconds\n", toc(start_time));
+    lp_status_t status = barrier_advanced_solve(start_time, solution, *owned_data);
+    return store_or_clear_cache(cache, owned_data, status);
+  } catch (const raft::cuda_error& e) {
+    settings.log.printf("Error in barrier_solver_t: %s\n", e.what());
+    return lp_status_t::NUMERICAL_ISSUES;
+  } catch (const std::bad_alloc& e) {
+    settings.log.printf("Out of memory in barrier_solver_t: %s\n", e.what());
+    return lp_status_t::NUMERICAL_ISSUES;
+  }
+}
+
+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)
+{
+  settings.log.printf("Barrier solver started at %.3f seconds\n", toc(start_time));
+  try {
+    raft::common::nvtx::range fun_scope("Barrier: solve");
+
+    i_t n = lp.num_cols;
+    i_t m = lp.num_rows;
+
+    solution.resize(m, n);
+    settings.log.printf(
+      "Barrier solver: %d constraints, %d variables, %ld nonzeros\n", m, n, lp.A.col_start[n]);
+
+    settings.log.printf("\n");
+
+    if (lp.Q.n > 0) {
+      settings.log.printf("Quadratic objective matrix  : %d nonzeros\n", lp.Q.row_start[lp.Q.n]);
+    }
+    if (lp.second_order_cone_dims.size() > 0) {
+      settings.log.printf("Second-order cones          : %d\n",
+                          static_cast<int>(lp.second_order_cone_dims.size()));
+    }
+
+    // Compute the number of free variables
+    i_t num_free_variables = presolve_info.free_variable_pairs.size() / 2;
+    if (num_free_variables > 0) {
+      settings.log.printf("Free variables              : %d\n", num_free_variables);
+    }
+
+    // Compute the number of upper bounds
+    i_t num_upper_bounds = 0;
+    for (i_t j = 0; j < n; j++) {
+      if (lp.upper[j] < inf) { num_upper_bounds++; }
+    }
+
+    std::unique_ptr<iteration_data_t<i_t, f_t>> owned_data;
+    if (cache != nullptr) { cache->store_iteration_data(nullptr); }
+
+    csc_matrix_t<i_t, f_t> local_Q(lp.num_cols, 0, 0);
+    csc_matrix_t<i_t, f_t>* Qin = &local_Q;
+    if (cache != nullptr && cache->transform() != nullptr) {
+      auto* xf      = cache->transform();
+      xf->barrier_Q = std::make_unique<csc_matrix_t<i_t, f_t>>(lp.num_cols, 0, 0);
+      Qin           = xf->barrier_Q.get();
+    }
+    if (lp.Q.n > 0) { create_Q(lp, *Qin); }
+    owned_data = std::make_unique<iteration_data_t<i_t, f_t>>(
+      lp, num_upper_bounds, presolve_info.direct_free_variables, *Qin, settings);
+    iteration_data_t<i_t, f_t>& data = *owned_data;
+
+    if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) {
+      settings.log.printf("Barrier solver halted\n");
+      if (cache != nullptr) { cache->clear(); }
+      return lp_status_t::CONCURRENT_LIMIT;
+    }
+    if (data.indefinite_Q) {
+      if (cache != nullptr) { cache->clear(); }
+      return lp_status_t::NUMERICAL_ISSUES;
+    }
+    if (data.symbolic_status != 0) {
+      settings.log.printf("Error in symbolic analysis\n");
+      if (cache != nullptr) { cache->clear(); }
+      return lp_status_t::NUMERICAL_ISSUES;
+    }
+
+    settings.log.printf("Barrier setup complete at %.3f seconds\n", toc(start_time));
+    lp_status_t status = barrier_advanced_solve(start_time, solution, *owned_data);
+    return store_or_clear_cache(cache, owned_data, status);
   } catch (const raft::cuda_error& e) {
     settings.log.printf("Error in barrier_solver_t: %s\n", e.what());
     return lp_status_t::NUMERICAL_ISSUES;
@@ -4713,6 +4887,38 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
   }
 }
 
+// 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,
+                                 std::unique_ptr<iteration_data_t<i_t, f_t>>& owned_data,
+                                 lp_status_t status)
+{
+  if (cache != nullptr && owned_data) {
+    if (status == lp_status_t::OPTIMAL) {
+      cache->store_iteration_data(owned_data.release());
+    } else {
+      cache->clear();
+    }
+  }
+  return status;
+}
+
+void destroy_iteration_data(iteration_data_t<int, double>* data) { delete data; }
+
+void apply_barrier_linear_objective(iteration_data_t<int, double>& data,
+                                    double const* barrier_c,
+                                    int n)
+{
+  cuopt_expects(barrier_c != nullptr && static_cast<int>(data.c.size()) == n &&
+                  static_cast<int>(data.d_c_.size()) == n,
+                error_type_t::ValidationError,
+                "update_linear_objective: barrier linear objective size does not match cached "
+                "iteration_data_t.");
+  std::copy(barrier_c, barrier_c + n, data.c.data());
+  raft::copy(
+    data.d_c_.data(), data.c.data(), static_cast<std::size_t>(n), data.handle_ptr->get_stream());
+}
+
 #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE
 template bool validate_barrier_cone_layout<int, double>(
   const lp_problem_t<int, double>& problem, const simplex_solver_settings_t<int, double>& settings);
diff --git a/cpp/src/barrier/barrier.hpp b/cpp/src/barrier/barrier.hpp
index 9865df693a..3d15803115 100644
--- a/cpp/src/barrier/barrier.hpp
+++ b/cpp/src/barrier/barrier.hpp
@@ -17,9 +17,16 @@
 #include <linear_algebra/sparse_matrix.hpp>
 #include <math_optimization/tic_toc.hpp>
 
+#include <memory>
+
 #include <rmm/device_uvector.hpp>
 
 #include <utility>
+
+namespace cuopt::mathematical_optimization {
+class barrier_cache_t;
+}
+
 namespace cuopt::mathematical_optimization::barrier {
 
 /** Validates SOC layout on an simplex::lp_problem_t before barrier presolve/solve. */
@@ -36,9 +43,20 @@ class barrier_solver_t {
   barrier_solver_t(const simplex::lp_problem_t<i_t, f_t>& lp,
                    const simplex::presolve_info_t<i_t, f_t>& presolve,
                    const simplex::simplex_solver_settings_t<i_t, f_t>& settings);
-  simplex::lp_status_t solve(f_t start_time, simplex::lp_solution_t<i_t, f_t>& solution);
+  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);
+  // 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);
 
  private:
+  simplex::lp_status_t barrier_advanced_solve(f_t start_time,
+                                              simplex::lp_solution_t<i_t, f_t>& solution,
+                                              iteration_data_t<i_t, f_t>& data);
   void my_pop_range(bool debug) const;
   void create_Q(const simplex::lp_problem_t<i_t, f_t>& lp, csc_matrix_t<i_t, f_t>& Q);
   int initial_point(iteration_data_t<i_t, f_t>& data);
diff --git a/cpp/src/barrier/barrier_cache.cu b/cpp/src/barrier/barrier_cache.cu
new file mode 100644
index 0000000000..4013eb690f
--- /dev/null
+++ b/cpp/src/barrier/barrier_cache.cu
@@ -0,0 +1,128 @@
+/* clang-format off */
+/*
+ * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-License-Identifier: Apache-2.0
+ */
+/* clang-format on */
+
+#include <cuopt/error.hpp>
+#include <cuopt/mathematical_optimization/utilities/barrier_cache.hpp>
+
+#include <barrier/barrier_transform.hpp>
+
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+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*)>;
+
+struct barrier_cache_t::impl {
+  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)),
+      iteration_data(nullptr, &barrier::destroy_iteration_data)
+  {
+  }
+
+  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;
+  bool c_dirty{false};
+};
+
+barrier_cache_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;
+
+barrier_cache_t::barrier_cache_t(barrier_cache_t&&) noexcept            = default;
+barrier_cache_t& barrier_cache_t::operator=(barrier_cache_t&&) noexcept = default;
+
+std::unique_ptr<barrier_cache_t> barrier_cache_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>(
+    new barrier_cache_t(std::move(stream), std::move(handle)));
+}
+
+raft::handle_t* barrier_cache_t::handle_ptr() { return impl_->handle.get(); }
+
+raft::handle_t const* barrier_cache_t::handle_ptr() const { return impl_->handle.get(); }
+
+void barrier_cache_t::clear()
+{
+  impl_->iteration_data.reset();
+  impl_->transform.reset();
+  impl_->c_dirty = false;
+}
+
+void barrier_cache_t::store_iteration_data(barrier_iteration_data_t* data)
+{
+  impl_->iteration_data.reset(data);
+}
+
+barrier_iteration_data_t* barrier_cache_t::release_iteration_data()
+{
+  return impl_->iteration_data.release();
+}
+
+void barrier_cache_t::store_transform(std::unique_ptr<barrier_transform_t> transform)
+{
+  impl_->transform = std::move(transform);
+}
+
+barrier_transform_t* barrier_cache_t::transform() { return impl_->transform.get(); }
+
+barrier_transform_t const* barrier_cache_t::transform() const { return impl_->transform.get(); }
+
+void barrier_cache_t::set_c_dirty(bool dirty) { impl_->c_dirty = dirty; }
+
+bool barrier_cache_t::c_dirty() const
+{
+  return impl_->c_dirty && impl_->transform != nullptr && impl_->iteration_data.get() != nullptr;
+}
+
+void barrier_cache_t::update_linear_objective(double const* c, int n)
+{
+  cuopt_expects(impl_->transform != nullptr,
+                error_type_t::ValidationError,
+                "update_linear_objective: no barrier transform; Solve with sequence_solve first.");
+  cuopt_expects(impl_->iteration_data.get() != nullptr,
+                error_type_t::ValidationError,
+                "update_linear_objective: no cached iteration_data; Solve a QP to Optimal first.");
+  std::vector<double> crushed;
+  try {
+    crushed = crush_user_linear_objective(*impl_->transform, c, n);
+  } catch (std::invalid_argument const& e) {
+    cuopt_expects(false, error_type_t::ValidationError, "%s", e.what());
+  }
+  if (impl_->transform->linear_obj_shift.size() == crushed.size()) {
+    for (std::size_t j = 0; j < crushed.size(); ++j) {
+      crushed[j] += impl_->transform->linear_obj_shift[j];
+    }
+  }
+  // The next solve builds its solver from barrier_lp, so keep its objective and the cached
+  // iteration workspace on the same c.
+  auto& barrier_objective = impl_->transform->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_->c_dirty = true;
+}
+
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/barrier/barrier_transform.hpp b/cpp/src/barrier/barrier_transform.hpp
new file mode 100644
index 0000000000..2f0fa953da
--- /dev/null
+++ b/cpp/src/barrier/barrier_transform.hpp
@@ -0,0 +1,111 @@
+/* clang-format off */
+/*
+ * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-License-Identifier: Apache-2.0
+ */
+/* clang-format on */
+
+#pragma once
+
+#include <dual_simplex/presolve.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+
+#include <memory>
+#include <stdexcept>
+#include <vector>
+
+namespace cuopt::mathematical_optimization {
+
+/**
+ * User-to-barrier transform retained on barrier_cache_t after Optimal:
+ * convert / presolve / scaling, plus the scaled LP.
+ * Enough to crush a new linear objective from the original problem into barrier
+ * coordinates and to uncrush a solution without rerunning those algorithms.
+ */
+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};
+
+  // 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;
+  // Barrier linear objective minus crush(user c) from the first solve (Q*ell shift, etc.).
+  std::vector<double> linear_obj_shift;
+  std::unique_ptr<cuopt::mathematical_optimization::simplex::lp_problem_t<int, double>> 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;
+};
+
+inline std::vector<double> crush_user_linear_objective(barrier_transform_t const& xf,
+                                                       double const* c,
+                                                       int n)
+{
+  if (c == nullptr || n != xf.user_num_cols) {
+    throw std::invalid_argument(
+      "update_linear_objective: linear objective length must match the cached user column count.");
+  }
+  if (xf.original_num_cols < xf.user_num_cols) {
+    throw std::invalid_argument(
+      "update_linear_objective: cached original column count is smaller than user n.");
+  }
+  if (xf.barrier_lp == nullptr) {
+    throw std::invalid_argument("update_linear_objective: cached barrier LP is missing.");
+  }
+
+  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];
+  }
+  for (int j : xf.presolve_info.negated_variables) {
+    orig[static_cast<std::size_t>(j)] *= -1.0;
+  }
+
+  std::vector<double> 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])];
+    }
+  } else {
+    presolved = std::move(orig);
+  }
+
+  auto const& pairs = xf.presolve_info.free_variable_pairs;
+  if (!pairs.empty()) {
+    if (pairs.size() % 2 != 0) {
+      throw std::invalid_argument("update_linear_objective: free_variable_pairs size is not even.");
+    }
+    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)];
+    }
+  }
+
+  if (static_cast<int>(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 / "
+      "column_scales.");
+  }
+  for (std::size_t j = 0; j < presolved.size(); ++j) {
+    presolved[j] /= xf.column_scales[j];
+  }
+  return presolved;
+}
+
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/barrier/device_sparse_matrix.cuh b/cpp/src/barrier/device_sparse_matrix.cuh
index 974e2b0f4a..8194bda6fc 100644
--- a/cpp/src/barrier/device_sparse_matrix.cuh
+++ b/cpp/src/barrier/device_sparse_matrix.cuh
@@ -179,6 +179,10 @@ class device_csc_matrix_t {
   {
   }
 
+  device_csc_matrix_t(device_csc_matrix_t&&)                 = default;
+  device_csc_matrix_t& operator=(device_csc_matrix_t&&)      = default;
+  device_csc_matrix_t& operator=(const device_csc_matrix_t&) = delete;
+
   device_csc_matrix_t(const csc_matrix_t<i_t, f_t>& A, rmm::cuda_stream_view stream)
     : m(A.m),
       n(A.n),
@@ -318,6 +322,10 @@ class device_csr_matrix_t {
   {
   }
 
+  device_csr_matrix_t(device_csr_matrix_t&&)                 = default;
+  device_csr_matrix_t& operator=(device_csr_matrix_t&&)      = default;
+  device_csr_matrix_t& operator=(const device_csr_matrix_t&) = delete;
+
   device_csr_matrix_t(const csr_matrix_t<i_t, f_t>& A, rmm::cuda_stream_view stream)
     : m(A.m),
       n(A.n),
diff --git a/cpp/src/barrier/sparse_cholesky.cuh b/cpp/src/barrier/sparse_cholesky.cuh
index 01045847d1..ecfbe95e09 100644
--- a/cpp/src/barrier/sparse_cholesky.cuh
+++ b/cpp/src/barrier/sparse_cholesky.cuh
@@ -33,6 +33,8 @@ class sparse_cholesky_base_t {
   virtual i_t solve(const dense_vector_t<i_t, f_t>& b, dense_vector_t<i_t, f_t>& x) = 0;
   virtual i_t solve(rmm::device_uvector<f_t>& b, rmm::device_uvector<f_t>& x)       = 0;
   virtual void set_positive_definite(bool positive_definite)                        = 0;
+  virtual void invalidate_numeric_factor() {}
+  virtual void rebind_csr_matrix(device_csr_matrix_t<i_t, f_t>& Arow) {}
 };
 
 #define CUDSS_EXAMPLE_FREE \
@@ -144,6 +146,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
       positive_definite(true),
       A_created(false),
       settings_(settings),
+      symbolic_done_(false),
       stream(handle_ptr->get_stream())
   {
     int major, minor, patch;
@@ -158,6 +161,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     if (CUDART_VERSION >= 13000 && settings_.concurrent_halt != nullptr &&
         settings_.num_gpus == 1) {
       cuGetErrorString_func = cuopt::get_driver_entry_point("cuGetErrorString");
+
       // 1. Set up the GPU resources
       CUdevResource initial_device_GPU_resources = {};
       auto cuDeviceGetDevResource_func = cuopt::get_driver_entry_point("cuDeviceGetDevResource");
@@ -493,7 +497,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
         return -1;
       }
       f_t reordering_time = toc(start_symbolic);
-      settings_.log.printf("Reordering time             : %.2fs\n", reordering_time);
+      settings_.log.printf("Reordering time             : %.3fs\n", reordering_time);
       start_symbolic_factor = tic();
 
       status = cudssExecute(
@@ -511,7 +515,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     }
     RAFT_CUDA_TRY(cudaStreamSynchronize(stream));
     f_t symbolic_factorization_time = toc(start_symbolic_factor);
-    settings_.log.printf("Symbolic factorization time : %.2fs\n", symbolic_factorization_time);
+    settings_.log.printf("Symbolic factorization time : %.3fs\n", symbolic_factorization_time);
     int64_t lu_nz       = 0;
     size_t size_written = 0;
     CUDSS_CALL_AND_CHECK(
@@ -524,12 +528,22 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     RAFT_CUDA_TRY(cudaStreamSynchronize(stream));
     handle_ptr_->get_stream().synchronize();
 
+    symbolic_done_ = true;
     return 0;
   }
   i_t factorize(device_csr_matrix_t<i_t, f_t>& Arow) override
   {
     raft::common::nvtx::range fun_scope("Factorize: cuDSS");
 
+    if (!symbolic_done_ || !A_created) {
+      settings_.log.printf(
+        "Error: cuDSS factorize(device_csr) called before analyze (symbolic_done=%d "
+        "A_created=%d)\n",
+        static_cast<int>(symbolic_done_),
+        static_cast<int>(A_created));
+      return -1;
+    }
+
 // #define PRINT_MATRIX_NORM
 #ifdef PRINT_MATRIX_NORM
     cudaStreamSynchronize(stream);
@@ -590,7 +604,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     }
 
     if (first_factor) {
-      settings_.log.debug("Factorization time          : %.2fs\n", numeric_time);
+      settings_.log.debug("Factorization time          : %.3fs\n", numeric_time);
       first_factor = false;
     }
     if (status != CUDSS_STATUS_SUCCESS) {
@@ -714,7 +728,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
 
     f_t symbolic_time = toc(start_symbolic);
     f_t analysis_time = toc(start_analysis);
-    settings_.log.printf("Symbolic factorization time : %.2fs\n", symbolic_time);
+    settings_.log.printf("Symbolic factorization time : %.3fs\n", symbolic_time);
     if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) {
       RAFT_CUDA_TRY(cudaStreamSynchronize(stream));
       handle_ptr_->get_stream().synchronize();
@@ -784,7 +798,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     }
 
     if (first_factor) {
-      settings_.log.debug("Factorization time          : %.2fs\n", numeric_time);
+      settings_.log.debug("Factorization time          : %.3fs\n", numeric_time);
       first_factor = false;
     }
     if (status != CUDSS_STATUS_SUCCESS) {
@@ -866,6 +880,55 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     this->positive_definite = positive_definite;
   }
 
+  /// Re-point cuDSS CSR wrapper at current device buffers after in-place value refresh.
+  void rebind_csr_matrix(device_csr_matrix_t<i_t, f_t>& Arow) override
+  {
+    if (!symbolic_done_ || !A_created) { return; }
+    auto d_nnz = Arow.row_start.element(Arow.m, Arow.row_start.stream());
+    if (d_nnz != nnz) { return; }
+    status = cudssMatrixDestroy(A);
+    if (status != CUDSS_STATUS_SUCCESS) {
+      settings_.log.printf("cudssMatrixDestroy for A rebind failed: %d\n", status);
+      return;
+    }
+#if CUDSS_VERSION_MAJOR > 0 || (CUDSS_VERSION_MAJOR == 0 && CUDSS_VERSION_MINOR >= 8)
+    status = cudssMatrixCreateCsr(&A,
+                                  n,
+                                  n,
+                                  nnz,
+                                  Arow.row_start.data(),
+                                  nullptr,
+                                  Arow.j.data(),
+                                  Arow.x.data(),
+                                  CUDSS_R_32I,
+                                  CUDSS_R_32I,
+                                  CUDSS_R_64F,
+                                  positive_definite ? CUDSS_MTYPE_SPD : CUDSS_MTYPE_SYMMETRIC,
+                                  CUDSS_MVIEW_FULL,
+                                  CUDSS_BASE_ZERO);
+#else
+    status = cudssMatrixCreateCsr(&A,
+                                  n,
+                                  n,
+                                  nnz,
+                                  Arow.row_start.data(),
+                                  nullptr,
+                                  Arow.j.data(),
+                                  Arow.x.data(),
+                                  CUDA_R_32I,
+                                  CUDA_R_64F,
+                                  positive_definite ? CUDSS_MTYPE_SPD : CUDSS_MTYPE_SYMMETRIC,
+                                  CUDSS_MVIEW_FULL,
+                                  CUDSS_BASE_ZERO);
+#endif
+    if (status != CUDSS_STATUS_SUCCESS) {
+      settings_.log.printf("cudssMatrixCreateCsr rebind failed: %d\n", status);
+      A_created = false;
+      return;
+    }
+    A_created = true;
+  }
+
  private:
   raft::handle_t const* handle_ptr_;
   i_t n;
@@ -889,7 +952,9 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
   f_t* x_values_d;
   f_t* b_values_d;
 
+  bool symbolic_done_;
   const simplex::simplex_solver_settings_t<i_t, f_t>& settings_;
+
   CUgreenCtx barrier_green_ctx;
   CUstream stream;
   void* cuGetErrorString_func;
diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp
index 388bb43b35..ae26585a0c 100644
--- a/cpp/src/dual_simplex/solve.cpp
+++ b/cpp/src/dual_simplex/solve.cpp
@@ -26,17 +26,75 @@
 #include <math_optimization/tic_toc.hpp>
 #include <math_optimization/types.hpp>
 
+#include <barrier/barrier_transform.hpp>
+#include <cuopt/mathematical_optimization/utilities/barrier_cache.hpp>
+
 #include <raft/core/nvtx.hpp>
 
 #include <cstdio>
 #include <cstdlib>
+#include <memory>
 #include <queue>
+#include <stdexcept>
 #include <string>
 
 namespace cuopt::mathematical_optimization::simplex {
 
 namespace {
 
+template <typename i_t, typename f_t>
+void unscale_uncrush_barrier_to_user(const user_problem_t<i_t, f_t>& user_problem,
+                                     const raft::handle_t* handle_ptr,
+                                     i_t original_num_rows,
+                                     i_t original_num_cols,
+                                     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,
+                                     const std::vector<f_t>& row_scales,
+                                     const simplex_solver_settings_t<i_t, f_t>& barrier_settings,
+                                     const lp_solution_t<i_t, f_t>& barrier_solution,
+                                     lp_solution_t<i_t, f_t>& solution)
+{
+  std::vector<f_t> unscaled_x(barrier_lp.num_cols);
+  std::vector<f_t> unscaled_y(barrier_lp.num_rows);
+  std::vector<f_t> unscaled_z(barrier_lp.num_cols);
+  unscale_solution<i_t, f_t>(column_scales,
+                             row_scales,
+                             barrier_solution.x,
+                             barrier_solution.y,
+                             barrier_solution.z,
+                             unscaled_x,
+                             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);
+  lp_solution_t<i_t, f_t> lp_solution(original_num_rows, original_num_cols);
+  uncrush_solution(presolve_info,
+                   barrier_settings,
+                   converted,
+                   unscaled_x,
+                   unscaled_y,
+                   unscaled_z,
+                   lp_solution.x,
+                   lp_solution.y,
+                   lp_solution.z);
+
+  uncrush_primal_solution(user_problem, converted, lp_solution.x, solution.x);
+  uncrush_dual_solution(
+    user_problem, converted, lp_solution.y, lp_solution.z, solution.y, solution.z);
+  solution.objective =
+    barrier_solution.user_objective / user_problem.obj_scale - user_problem.obj_constant;
+  solution.user_objective     = barrier_solution.user_objective;
+  solution.l2_primal_residual = barrier_solution.l2_primal_residual;
+  solution.l2_dual_residual   = barrier_solution.l2_dual_residual;
+  solution.iterations         = barrier_solution.iterations;
+}
+
+}  // namespace
+
+namespace {
+
 template <typename i_t, typename f_t>
 void write_matlab(const std::string& filename, const simplex::lp_problem_t<i_t, f_t>& lp)
 {
@@ -354,18 +412,55 @@ lp_status_t solve_linear_program_with_advanced_basis(
 }
 
 template <typename i_t, typename f_t>
-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,
-                                              f_t start_time,
-                                              lp_solution_t<i_t, f_t>& solution,
-                                              const raft::handle_t* handle_ptr)
+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,
+  f_t start_time,
+  lp_solution_t<i_t, f_t>& solution,
+  cuopt::mathematical_optimization::barrier_cache_t* cache,
+  const raft::handle_t* handle_ptr)
 {
-  lp_status_t status = lp_status_t::UNSET;
+  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->c_dirty()) ? cache->transform() : nullptr;
+  const bool reuse_c_only =
+    xf != nullptr && xf->barrier_lp != nullptr && !user_problem.Q_values.empty() &&
+    user_problem.second_order_cone_dims.empty() && xf->second_order_cone_dims.empty() &&
+    xf->barrier_lp->second_order_cone_dims.empty() &&
+    settings.barrier_presolve_bound_free_variables == 0 &&
+    user_problem.num_cols == xf->user_num_cols && user_problem.num_rows == xf->user_num_rows;
+
+  if (reuse_c_only) {
+    settings.log.printf("Barrier: reusing cache (skip convert/presolve/scaling)\n");
+    lp_solution_t<i_t, f_t> barrier_solution(xf->barrier_lp->num_rows, xf->barrier_lp->num_cols);
+    barrier::barrier_solver_t<i_t, f_t> barrier_solver(
+      *xf->barrier_lp, xf->presolve_info, barrier_settings);
+    lp_status_t barrier_status =
+      barrier_solver.solve_with_cache(start_time, barrier_solution, cache);
+    if (barrier_status == lp_status_t::OPTIMAL) {
+      unscale_uncrush_barrier_to_user(user_problem,
+                                      cache->handle_ptr(),
+                                      xf->original_num_rows,
+                                      xf->original_num_cols,
+                                      *xf->barrier_lp,
+                                      xf->presolve_info,
+                                      xf->column_scales,
+                                      xf->row_scales,
+                                      barrier_settings,
+                                      barrier_solution,
+                                      solution);
+      cache->set_c_dirty(false);
+    } else {
+      cache->clear();
+    }
+    return barrier_status;
+  }
+
   lp_problem_t<i_t, f_t> original_lp(handle_ptr, 1, 1, 1);
 
   // Convert the user problem to a linear program with only equality constraints
   std::vector<i_t> new_slacks;
-  simplex_solver_settings_t<i_t, f_t> barrier_settings = settings;
   dualize_info_t<i_t, f_t> dualize_info;
   convert_user_problem(user_problem, barrier_settings, original_lp, new_slacks, dualize_info);
   if (!barrier::validate_barrier_cone_layout(original_lp, barrier_settings)) {
@@ -377,7 +472,7 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& us
   // Presolve the linear program
   presolve_info_t<i_t, f_t> presolve_info;
   lp_problem_t<i_t, f_t> presolved_lp(handle_ptr, 1, 1, 1);
-  const i_t ok = presolve(original_lp, barrier_settings, presolved_lp, presolve_info);
+  i_t ok = presolve(original_lp, barrier_settings, presolved_lp, presolve_info);
   if (ok == CONCURRENT_HALT_RETURN) { return lp_status_t::CONCURRENT_LIMIT; }
   if (ok == TIME_LIMIT_RETURN) { return lp_status_t::TIME_LIMIT; }
   if (ok == -1) { return lp_status_t::INFEASIBLE; }
@@ -394,8 +489,56 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& us
   // Solve using barrier
   lp_solution_t<i_t, f_t> barrier_solution(barrier_lp.num_rows, barrier_lp.num_cols);
 
-  barrier::barrier_solver_t<i_t, f_t> barrier_solver(barrier_lp, presolve_info, barrier_settings);
-  lp_status_t barrier_status = barrier_solver.solve(start_time, barrier_solution);
+  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;
+    xf->original_num_cols            = original_lp.num_cols;
+    xf->original_num_rows            = original_lp.num_rows;
+    xf->obj_scale                    = user_problem.obj_scale;
+    xf->obj_constant                 = user_problem.obj_constant;
+    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->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->barrier_lp                   = std::make_unique<lp_problem_t<i_t, f_t>>(barrier_lp);
+    solver_lp                        = xf->barrier_lp.get();
+    cache->store_transform(std::move(xf));
+  }
+
+  barrier::barrier_solver_t<i_t, f_t> barrier_solver(*solver_lp, presolve_info, barrier_settings);
+  lp_status_t barrier_status = barrier_solver.solve(start_time, barrier_solution, cache);
+
+  if (cache != nullptr) {
+    if (barrier_status == lp_status_t::OPTIMAL) {
+      auto* xf = cache->transform();
+      {
+        try {
+          auto crushed = cuopt::mathematical_optimization::crush_user_linear_objective(
+            *xf, user_problem.objective.data(), user_problem.num_cols);
+          xf->linear_obj_shift.resize(static_cast<std::size_t>(solver_lp->num_cols), 0.0);
+          if (static_cast<int>(crushed.size()) == solver_lp->num_cols) {
+            for (int j = 0; j < solver_lp->num_cols; ++j) {
+              xf->linear_obj_shift[static_cast<std::size_t>(j)] =
+                solver_lp->objective[static_cast<std::size_t>(j)] -
+                crushed[static_cast<std::size_t>(j)];
+            }
+          }
+        } catch (std::exception const&) {
+          xf->linear_obj_shift.assign(static_cast<std::size_t>(solver_lp->num_cols), 0.0);
+        }
+      }
+    } else {
+      cache->clear();
+    }
+  }
+
   if (barrier_status == lp_status_t::OPTIMAL) {
 #ifdef COMPUTE_SCALED_RESIDUALS
     std::vector<f_t> scaled_residual = barrier_lp.rhs;
@@ -686,22 +829,27 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& us
 }
 
 template <typename i_t, typename f_t>
-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,
-                                              f_t start_time,
-                                              lp_solution_t<i_t, f_t>& solution)
+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)
 {
+  f_t start_time = tic();
   return solve_linear_program_with_barrier(
-    user_problem, settings, start_time, solution, user_problem.handle_ptr);
+    user_problem, settings, start_time, solution, cache, user_problem.handle_ptr);
 }
 
 template <typename i_t, typename f_t>
-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)
+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,
+  f_t start_time,
+  lp_solution_t<i_t, f_t>& solution,
+  cuopt::mathematical_optimization::barrier_cache_t* cache)
 {
-  f_t start_time = tic();
-  return solve_linear_program_with_barrier(user_problem, settings, start_time, solution);
+  return solve_linear_program_with_barrier(
+    user_problem, settings, start_time, solution, cache, user_problem.handle_ptr);
 }
 
 template <typename i_t, typename f_t>
@@ -845,19 +993,22 @@ template lp_status_t solve_linear_program_with_advanced_basis(
 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);
+  lp_solution_t<int, double>& solution,
+  cuopt::mathematical_optimization::barrier_cache_t* 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);
+  lp_solution_t<int, double>& solution,
+  cuopt::mathematical_optimization::barrier_cache_t* 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,
   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 308c462de5..f8966e29c1 100644
--- a/cpp/src/dual_simplex/solve.hpp
+++ b/cpp/src/dual_simplex/solve.hpp
@@ -17,6 +17,10 @@ namespace cuopt {
 struct work_limit_context_t;
 }
 
+namespace cuopt::mathematical_optimization {
+class barrier_cache_t;
+}  // namespace cuopt::mathematical_optimization
+
 namespace cuopt::mathematical_optimization::simplex {
 
 template <typename i_t, typename f_t>
@@ -91,22 +95,28 @@ lp_status_t solve_linear_program_with_advanced_basis(
   work_limit_context_t* work_unit_context = nullptr);
 
 template <typename i_t, typename f_t>
-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);
+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);
 
 template <typename i_t, typename f_t>
-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,
-                                              f_t start_time,
-                                              lp_solution_t<i_t, f_t>& solution);
+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,
+  f_t start_time,
+  lp_solution_t<i_t, f_t>& solution,
+  cuopt::mathematical_optimization::barrier_cache_t* cache = nullptr);
 
 template <typename i_t, typename f_t>
-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,
-                                              f_t start_time,
-                                              lp_solution_t<i_t, f_t>& solution,
-                                              const raft::handle_t* handle_ptr);
+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,
+  f_t start_time,
+  lp_solution_t<i_t, f_t>& solution,
+  cuopt::mathematical_optimization::barrier_cache_t* cache,
+  const raft::handle_t* handle_ptr);
 
 template <typename i_t, typename f_t>
 lp_status_t solve_linear_program(const user_problem_t<i_t, f_t>& user_problem,
diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu
index d08a36d178..cf92e27447 100644
--- a/cpp/src/pdlp/solve.cu
+++ b/cpp/src/pdlp/solve.cu
@@ -26,6 +26,7 @@
 #include <mip_heuristics/presolve/trivial_presolve.cuh>
 #include <mip_heuristics/solver.cuh>
 
+#include <barrier/barrier_transform.hpp>
 #include <cuopt/mathematical_optimization/backend_selection.hpp>
 #include <cuopt/mathematical_optimization/cpu_optimization_problem.hpp>
 #include <cuopt/mathematical_optimization/cpu_optimization_problem_solution.hpp>
@@ -35,6 +36,7 @@
 #include <cuopt/mathematical_optimization/pdlp/pdlp_hyper_params.cuh>
 #include <cuopt/mathematical_optimization/pdlp/solver_settings.hpp>
 #include <cuopt/mathematical_optimization/solve.hpp>
+#include <cuopt/mathematical_optimization/utilities/barrier_cache.hpp>
 
 #include <cuopt/mathematical_optimization/io/mps_data_model.hpp>
 #include <utilities/copy_helpers.hpp>
@@ -62,6 +64,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <cstdint>
 #include <exception>
 #include <optional>
 #include <set>
@@ -72,6 +75,33 @@
 
 namespace cuopt::mathematical_optimization {
 
+namespace {
+
+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)
+{
+  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.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);
+  user_problem.obj_constant = static_cast<f_t>(xf.obj_constant);
+  // Nonempty Q so the cache-reuse path accepts this as a QP (it rejects empty Q).
+  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_col_to_expanded_col = xf.original_col_to_expanded_col;
+  return user_problem;
+}
+
+}  // namespace
+
 template <typename From, typename To>
 extern rmm::device_uvector<To> gpu_cast(const rmm::device_uvector<From>& src,
                                         rmm::cuda_stream_view stream);
@@ -492,7 +522,8 @@ std::tuple<simplex::lp_solution_t<i_t, f_t>, simplex::lp_status_t, f_t, f_t, f_t
   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)
+  const raft::handle_t* handle_ptr,
+  cuopt::mathematical_optimization::barrier_cache_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);
@@ -533,7 +564,7 @@ std::tuple<simplex::lp_solution_t<i_t, f_t>, simplex::lp_status_t, f_t, f_t, f_t
 
   simplex::lp_solution_t<i_t, f_t> solution(user_problem.num_rows, user_problem.num_cols);
   auto status = simplex::solve_linear_program_with_barrier<i_t, f_t>(
-    user_problem, barrier_settings, timer.get_tic_start(), solution, handle_ptr);
+    user_problem, barrier_settings, timer.get_tic_start(), solution, cache, handle_ptr);
 
   if (status == simplex::lp_status_t::OPTIMAL) {
     barrier::project_barrier_solution_to_model_variables(user_problem, solution);
@@ -557,12 +588,14 @@ template <typename i_t, typename f_t>
 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)
+  const timer_t& timer,
+  cuopt::mathematical_optimization::barrier_cache_t* cache = nullptr)
 {
   // Convert data structures to dual simplex format and back
   simplex::user_problem_t<i_t, f_t> dual_simplex_problem =
     cuopt_problem_to_user_problem<i_t, f_t>(problem.handle_ptr, problem, false);
-  auto sol_dual_simplex = run_barrier(dual_simplex_problem, settings, timer, problem.handle_ptr);
+  auto sol_dual_simplex =
+    run_barrier(dual_simplex_problem, settings, timer, problem.handle_ptr, cache);
   return convert_dual_simplex_sol(problem,
                                   std::get<0>(sol_dual_simplex),
                                   std::get<1>(sol_dual_simplex),
@@ -1813,7 +1846,7 @@ optimization_problem_solution_t<i_t, f_t> solve_lp_with_method(
     if (settings.method == method_t::DualSimplex) {
       return run_dual_simplex(problem, settings, timer);
     } else if (settings.method == method_t::Barrier) {
-      return run_barrier(problem, settings, timer);
+      return run_barrier(problem, settings, timer, settings.barrier_cache);
     } else if (settings.method == method_t::Concurrent) {
       return run_concurrent(problem, settings, timer, is_batch_mode);
     } else {
@@ -1845,7 +1878,17 @@ optimization_problem_solution_t<i_t, f_t> solve_qcqp(
 
     auto qcqp_timer = cuopt::timer_t(settings.time_limit);
 
-    if (problem_checking) {
+    auto* cache    = settings.barrier_cache;
+    auto const* xf = (cache != nullptr && cache->c_dirty()) ? cache->transform() : nullptr;
+    const bool reuse_from_cache =
+      settings.user_problem_file.empty() && xf != nullptr && xf->barrier_lp != nullptr &&
+      settings.barrier_presolve_bound_free_variables == 0 && op_problem.has_quadratic_objective() &&
+      !op_problem.has_quadratic_constraints() && xf->second_order_cone_dims.empty() &&
+      static_cast<int>(xf->row_sense.size()) == xf->user_num_rows &&
+      op_problem.get_n_variables() == xf->user_num_cols &&
+      op_problem.get_n_constraints() == xf->user_num_rows;
+
+    if (problem_checking && !reuse_from_cache) {
       problem_checking_t<i_t, f_t>::check_problem_representation(op_problem);
       if (problem_checking_t<i_t, f_t>::has_crossing_bounds(op_problem)) {
         return optimization_problem_solution_t<i_t, f_t>(
@@ -1873,12 +1916,20 @@ optimization_problem_solution_t<i_t, f_t> solve_qcqp(
       CUOPT_LOG_INFO("Writing user problem to file: %s", settings.user_problem_file.c_str());
       op_problem.write_to_mps(settings.user_problem_file);
     }
-    // Convert data structures to dual simplex format and back
-    simplex::user_problem_t<i_t, f_t> dual_simplex_problem =
-      cuopt_optimization_problem_to_user_problem<i_t, f_t>(op_problem.get_handle_ptr(), op_problem);
-    auto sol_dual_simplex =
-      run_barrier(dual_simplex_problem, settings, qcqp_timer, op_problem.get_handle_ptr());
-    auto solution = convert_dual_simplex_sol(op_problem,
+    simplex::user_problem_t<i_t, f_t> dual_simplex_problem(op_problem.get_handle_ptr());
+    if (reuse_from_cache) {
+      dual_simplex_problem = user_problem_from_transform(
+        op_problem.get_handle_ptr(), op_problem, *settings.barrier_cache->transform());
+    } else {
+      dual_simplex_problem = cuopt_optimization_problem_to_user_problem<i_t, f_t>(
+        op_problem.get_handle_ptr(), op_problem);
+    }
+    auto sol_dual_simplex = run_barrier(dual_simplex_problem,
+                                        settings,
+                                        qcqp_timer,
+                                        op_problem.get_handle_ptr(),
+                                        settings.barrier_cache);
+    auto solution         = convert_dual_simplex_sol(op_problem,
                                              std::get<0>(sol_dual_simplex),
                                              std::get<1>(sol_dual_simplex),
                                              std::get<2>(sol_dual_simplex),
diff --git a/cpp/src/pdlp/utilities/cython_solve.cu b/cpp/src/pdlp/utilities/cython_solve.cu
index ed77c4f722..b106eb5e0e 100644
--- a/cpp/src/pdlp/utilities/cython_solve.cu
+++ b/cpp/src/pdlp/utilities/cython_solve.cu
@@ -6,6 +6,7 @@
 /* clang-format on */
 
 #include <cuopt/error.hpp>
+
 #include <cuopt/mathematical_optimization/backend_selection.hpp>
 #include <cuopt/mathematical_optimization/cpu_optimization_problem.hpp>
 #include <cuopt/mathematical_optimization/cpu_optimization_problem_solution.hpp>
@@ -17,7 +18,9 @@
 #include <cuopt/mathematical_optimization/optimization_problem_utils.hpp>
 #include <cuopt/mathematical_optimization/solve.hpp>
 #include <cuopt/mathematical_optimization/solver_settings.hpp>
+#include <cuopt/mathematical_optimization/utilities/barrier_cache.hpp>
 #include <cuopt/mathematical_optimization/utilities/cython_solve.hpp>
+
 #include <mip_heuristics/logger.hpp>
 #include <utilities/copy_helpers.hpp>
 #include <utilities/logger.hpp>
@@ -30,11 +33,15 @@
 #include <utility>
 #include <vector>
 
+#include <chrono>
+
 #include <unistd.h>
 
 namespace cuopt {
 namespace cython {
 
+using mathematical_optimization::barrier_cache_t;
+
 /**
  * @brief Wrapper for linear_programming to expose the API to cython
  *
@@ -96,24 +103,57 @@ std::unique_ptr<solver_ret_t> call_solve(
   cuopt::mathematical_optimization::io::data_model_view_t<int, double>* data_model,
   cuopt::mathematical_optimization::solver_settings_t<int, double>* solver_settings,
   unsigned int flags,
-  bool is_batch_mode)
+  bool is_batch_mode,
+  barrier_cache_t* cache_in)
 {
   raft::common::nvtx::range fun_scope("Call Solve");
 
+  cuopt_expects(
+    data_model != nullptr, error_type_t::ValidationError, "call_solve: data_model is null.");
+  cuopt_expects(solver_settings != nullptr,
+                error_type_t::ValidationError,
+                "call_solve: solver_settings is null.");
+
   // Determine memory backend based on execution mode
   auto memory_backend = cuopt::mathematical_optimization::get_memory_backend_type();
 
   solver_ret_t response;
 
-  // Create problem instance and CUDA resources based on memory backend
+  auto& pdlp_settings       = solver_settings->get_pdlp_settings();
+  const bool sequence_solve = pdlp_settings.sequence_solve;
+  const bool barrier_path =
+    data_model->has_quadratic_objective() || data_model->has_quadratic_constraints() ||
+    pdlp_settings.method == cuopt::mathematical_optimization::method_t::Barrier;
+  const bool want_cache =
+    (cache_in != nullptr || sequence_solve) && barrier_path &&
+    memory_backend == cuopt::mathematical_optimization::memory_backend_t::GPU && !is_batch_mode;
+
+  std::unique_ptr<barrier_cache_t> owned_cache;
+  barrier_cache_t* active_cache = cache_in;
+  pdlp_settings.barrier_cache   = nullptr;
+
+  // Create problem instance and CUDA resources based on memory backend.
+  // Do not construct rmm::cuda_stream until we know we are on GPU: CPU-only /
+  // remote-gRPC hosts have no device (CUDA_VISIBLE_DEVICES="") and stream
+  // construction would throw cudaErrorNoDevice.
   if (memory_backend == cuopt::mathematical_optimization::memory_backend_t::GPU) {
-    // GPU memory backend: Create CUDA resources and GPU problem
-    rmm::cuda_stream stream(static_cast<rmm::cuda_stream::flags>(flags));
-    const raft::handle_t handle_{stream};
+    rmm::cuda_stream ephemeral_stream(static_cast<rmm::cuda_stream::flags>(flags));
+    raft::handle_t ephemeral_handle(ephemeral_stream);
+    raft::handle_t* solve_handle = &ephemeral_handle;
+
+    if (want_cache) {
+      if (active_cache == nullptr) {
+        owned_cache  = barrier_cache_t::create(flags);
+        active_cache = owned_cache.get();
+      }
+      solve_handle                = active_cache->handle_ptr();
+      pdlp_settings.barrier_cache = active_cache;
+    }
 
-    auto problem = cuopt::mathematical_optimization::optimization_problem_t<int, double>(&handle_);
+    auto problem =
+      cuopt::mathematical_optimization::optimization_problem_t<int, double>(solve_handle);
     cuopt::mathematical_optimization::populate_from_data_model_view(
-      &problem, data_model, solver_settings, &handle_);
+      &problem, data_model, solver_settings, solve_handle);
 
     // Call appropriate solve function and convert to ret struct
     if (problem.get_problem_category() == mathematical_optimization::problem_category_t::LP) {
@@ -142,6 +182,8 @@ std::unique_ptr<solver_ret_t> call_solve(
       gpu_sols.last_restart_duality_gap_primal_solution_->set_stream(rmm::cuda_stream_per_thread);
       gpu_sols.last_restart_duality_gap_dual_solution_->set_stream(rmm::cuda_stream_per_thread);
 
+      if (owned_cache) { response.lp_ret.barrier_cache = std::move(owned_cache); }
+
     } else {
       // MIP solve
       auto mip_solution_ptr =
@@ -200,6 +242,8 @@ std::unique_ptr<solver_ret_t> call_solve(
     }
   }
 
+  pdlp_settings.barrier_cache = nullptr;
+
   return std::make_unique<solver_ret_t>(std::move(response));
 }
 
@@ -288,7 +332,8 @@ std::pair<std::vector<std::unique_ptr<solver_ret_t>>, double> call_batch_solve(
 
 #pragma omp parallel for num_threads(max_thread)
   for (std::size_t i = 0; i < size; ++i)
-    list[i] = call_solve(data_models[i], solver_settings, cudaStreamNonBlocking, is_batch_mode);
+    list[i] =
+      call_solve(data_models[i], solver_settings, cudaStreamNonBlocking, is_batch_mode, nullptr);
 
   auto end      = std::chrono::high_resolution_clock::now();
   auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start_solver);
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 7bcdfaea9b..b8712c4bf0 100644
--- a/python/cuopt/cuopt/linear_programming/data_model/data_model.py
+++ b/python/cuopt/cuopt/linear_programming/data_model/data_model.py
@@ -228,6 +228,24 @@ def set_objective_coefficients(self, c):
         """
         super().set_objective_coefficients(c)
 
+    @catch_cuopt_exception
+    def update_linear_objective(self, c):
+        """
+        Update the linear objective coefficients (c) for a sequence re-solve.
+
+        Writes ``c`` onto this DataModel. If a barrier cache is
+        present, also maps ``c`` into the cached barrier workspace and marks
+        it dirty (quadratic ``Q``, ``A``, and bounds must stay unchanged).
+        Cache reuse is QP-only: quadratic constraints take a full solve.
+
+        Parameters
+        ----------
+        c : array-like of float64
+            Linear objective coefficients, length equal to the number of
+            variables on the first ``sequence_solve``.
+        """
+        super().update_linear_objective(c)
+
     @catch_cuopt_exception
     def set_objective_scaling_factor(self, objective_scaling_factor):
         """
diff --git a/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pxd b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pxd
index 6c401b59f5..7633fb69f7 100644
--- a/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pxd
+++ b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pxd
@@ -1,4 +1,4 @@
-# SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # noqa
+# SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
 # SPDX-License-Identifier: Apache-2.0
 
 from .data_model cimport *
@@ -12,6 +12,7 @@ from libcpp.memory cimport unique_ptr
 
 cdef class DataModel:
     cdef unique_ptr[data_model_view_t[int, double]] c_data_model_view
+    cdef object barrier_cache_capsule
 
     cdef void _set_cpp_quadratic_constraints(
         self, data_model_view_t[int, double]* c_data_model_view
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 b298cf57a5..8fed4296f6 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
@@ -21,6 +21,14 @@ from libcpp.string cimport string
 from libcpp.utility cimport move
 from libcpp.vector cimport vector
 
+cdef extern from "Python.h":
+    bint PyCapsule_IsValid(object cap, const char* name)
+    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 +
+
 
 def type_cast(np_obj, np_type, name):
     if not isinstance(np_obj, np.ndarray):
@@ -41,6 +49,7 @@ cdef class DataModel:
 
     def __init__(self):
         self.c_data_model_view.reset(new data_model_view_t[int, double]())
+        self.barrier_cache_capsule = None
 
         self.maximize = False
         self.A_values = np.array([])
@@ -158,6 +167,34 @@ cdef class DataModel:
     def set_objective_coefficients(self, c):
         self.c = type_cast(c, np.float64, "c")
 
+    def update_linear_objective(self, c):
+        """Update linear objective coefficients (``c`` on this DataModel).
+
+        Always writes the DataModel objective. If this model owns a solver
+        cache from a prior Barrier solve, also crushes ``c`` into the cached
+        ``iteration_data_t`` and sets ``c_dirty`` 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 double[::1] c_view
+        new_c = type_cast(c, np.float64, "c")
+        if self.barrier_cache_capsule is not None:
+            if not PyCapsule_IsValid(
+                self.barrier_cache_capsule, b"cuopt.barrier_cache"
+            ):
+                raise ValueError("Invalid barrier cache stored on DataModel.")
+            cache = <barrier_cache_t*>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)
+            else:
+                cache.update_linear_objective(&c_view[0], <int>c_view.shape[0])
+        self.c = new_c
+
     def set_objective_scaling_factor(self, objective_scaling_factor):
         self.objective_scaling_factor = objective_scaling_factor
 
diff --git a/python/cuopt/cuopt/linear_programming/solver/solver.pxd b/python/cuopt/cuopt/linear_programming/solver/solver.pxd
index 04b75ce4e5..1b673a363d 100644
--- a/python/cuopt/cuopt/linear_programming/solver/solver.pxd
+++ b/python/cuopt/cuopt/linear_programming/solver/solver.pxd
@@ -1,4 +1,4 @@
-# SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # noqa
+# SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
 # SPDX-License-Identifier: Apache-2.0
 
 
@@ -8,6 +8,7 @@
 # cython: language_level = 3
 
 from libcpp cimport bool
+from libcpp.memory cimport unique_ptr
 from libcpp.pair cimport pair
 from libcpp.string cimport string
 from libcpp.vector cimport vector
@@ -93,6 +94,10 @@ cdef extern from "cuopt/mathematical_optimization/utilities/cython_types.hpp" na
         vector[double] last_restart_duality_gap_primal_solution_
         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:
+        pass
+
 cdef extern from "cuopt/mathematical_optimization/utilities/cython_solve.hpp" namespace "cuopt::cython": # noqa
     # Unified LP solution struct — solutions_ variant accessed via helpers
     cdef cppclass linear_programming_ret_t:
@@ -117,6 +122,7 @@ cdef extern from "cuopt/mathematical_optimization/utilities/cython_solve.hpp" na
         int nb_iterations_
         double solve_time_
         method_t solved_by_
+        unique_ptr[barrier_cache_t] barrier_cache
         bool is_gpu()
 
     # Unified MIP solution struct — solution_ variant accessed via helpers
@@ -144,6 +150,9 @@ cdef extern from "cuopt/mathematical_optimization/utilities/cython_solve.hpp" na
     cdef unique_ptr[solver_ret_t] call_solve(
         data_model_view_t[int, double]* data_model,
         solver_settings_t[int, double]* solver_settings,
+        unsigned int flags,
+        bool is_batch_mode,
+        barrier_cache_t* 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 7a27a87140..4246d7ee6c 100644
--- a/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx
+++ b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx
@@ -18,7 +18,7 @@ from dateutil.relativedelta import relativedelta
 
 from cuopt.utilities import type_cast
 
-from libc.stdint cimport uintptr_t
+from libc.stdint cimport uintptr_t, uint32_t
 from libc.stdlib cimport free, malloc
 from libc.string cimport memcpy, strcpy, strlen
 from libcpp cimport bool
@@ -28,6 +28,13 @@ from libcpp.string cimport string
 from libcpp.utility cimport move
 from libcpp.vector cimport vector
 
+from cpython.pycapsule cimport (
+    PyCapsule_Destructor,
+    PyCapsule_GetPointer,
+    PyCapsule_IsValid,
+    PyCapsule_New,
+)
+
 from rmm.pylibrmm.device_buffer cimport DeviceBuffer
 
 from cuopt.linear_programming.data_model.data_model cimport data_model_view_t
@@ -43,6 +50,7 @@ from cuopt.linear_programming.solver.solver cimport (
     linear_programming_ret_t,
     lp_cpu_solutions_t,
     lp_gpu_solutions_t,
+    barrier_cache_t,
     mip_ret_t,
     mip_termination_status_t,
     pdlp_solver_mode_t,
@@ -79,6 +87,25 @@ cdef extern from "cuopt/mathematical_optimization/utilities/internals.hpp" names
     cdef cppclass base_solution_callback_t
 
 
+cdef extern from *:
+    """
+    #include <cuopt/mathematical_optimization/utilities/barrier_cache.hpp>
+
+    static void cuopt_barrier_cache_capsule_dtor(PyObject *cap) noexcept
+    {
+      void *p = PyCapsule_GetPointer(cap, "cuopt.barrier_cache");
+      if (p != nullptr) {
+        delete reinterpret_cast<cuopt::mathematical_optimization::barrier_cache_t *>(p);
+      }
+    }
+    """
+    void cuopt_barrier_cache_capsule_dtor(object cap) noexcept
+
+
+cdef extern from "driver_types.h":
+    cdef uint32_t cudaStreamNonBlocking
+
+
 class MILPTerminationStatus(IntEnum):
     NoTermination = mip_termination_status_t.NoTermination
     Optimal = mip_termination_status_t.Optimal
@@ -341,7 +368,7 @@ cdef create_solution_with_names(unique_ptr[solver_ret_t] sol_ret_ptr,
     if sol_ret.problem_type == ProblemCategory.MIP or sol_ret.problem_type == ProblemCategory.IP: # noqa
         mip_ptr = &sol_ret.mip_ret
 
-        # Extract solution vector — branch only for the buffer type
+        # Extract solution vector -- branch only for the buffer type
         if mip_ptr.is_gpu():
             solution_buf = DeviceBuffer.c_from_unique_ptr(move(get_gpu_mip_solution(mip_ptr[0]))) # noqa
             solution = series_from_buf(solution_buf, pa.float64()).to_numpy()
@@ -371,7 +398,7 @@ cdef create_solution_with_names(unique_ptr[solver_ret_t] sol_ret_ptr,
     else:
         lp_ptr = &sol_ret.lp_ret
 
-        # Extract solution vectors — branch only for the buffer type
+        # Extract solution vectors -- branch only for the buffer type
         if lp_ptr.is_gpu():
             gpu_sols = &get_gpu_lp_solutions(lp_ptr[0])
 
@@ -452,7 +479,7 @@ cdef create_solution_with_names(unique_ptr[solver_ret_t] sol_ret_ptr,
                     last_restart_primal = None
                     last_restart_dual = None
 
-        # Shared scalar access — written once regardless of GPU/CPU backend
+        # Shared scalar access -- written once regardless of GPU/CPU backend
         if not is_batch:
             return Solution(
                 ProblemCategory(sol_ret.problem_type),
@@ -526,6 +553,19 @@ 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 solver_ret_t* sol_ret
+
+    if settings.sequence_solve and data_model_obj.barrier_cache_capsule is not None:
+        if not PyCapsule_IsValid(
+            data_model_obj.barrier_cache_capsule,
+            b"cuopt.barrier_cache",
+        ):
+            raise ValueError("Invalid barrier cache stored on DataModel.")
+        cache_in = <barrier_cache_t*>PyCapsule_GetPointer(
+            data_model_obj.barrier_cache_capsule,
+            b"cuopt.barrier_cache",
+        )
 
     data_model_obj.variable_types = type_cast(
         data_model_obj.variable_types, "S1", "variable_types"
@@ -541,7 +581,22 @@ def Solve(py_data_model_obj, SolverSettings settings, mip=False):
         sol_ret_ptr = move(call_solve(
             data_model_obj.c_data_model_view.get(),
             settings.c_solver_settings.get(),
+            cudaStreamNonBlocking,
+            False,
+            cache_in,
         ))
+
+    sol_ret = sol_ret_ptr.get()
+    if (
+        sol_ret.problem_type == ProblemCategory.LP
+        and sol_ret.lp_ret.barrier_cache.get() != NULL
+    ):
+        data_model_obj.barrier_cache_capsule = PyCapsule_New(
+            <void*>sol_ret.lp_ret.barrier_cache.release(),
+            b"cuopt.barrier_cache",
+            <PyCapsule_Destructor>cuopt_barrier_cache_capsule_dtor,
+        )
+
     return create_solution(move(sol_ret_ptr), data_model_obj)
 
 
diff --git a/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pxd b/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pxd
index 03958d2286..51b113f0f6 100644
--- a/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pxd
+++ b/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pxd
@@ -1,4 +1,4 @@
-# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # noqa
+# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
 # SPDX-License-Identifier: Apache-2.0
 
 
@@ -15,6 +15,10 @@ from libcpp.vector cimport vector
 cdef extern from "cuopt/mathematical_optimization/utilities/internals.hpp" namespace "cuopt::internals": # noqa
     cdef cppclass base_solution_callback_t
 
+cdef extern from "cuopt/mathematical_optimization/utilities/barrier_cache.hpp" namespace "cuopt::mathematical_optimization": # noqa
+    cdef cppclass barrier_cache_t:
+        pass
+
 cdef extern from "cuopt/mathematical_optimization/pdlp/solver_settings.hpp" namespace "cuopt::mathematical_optimization": # noqa
     ctypedef enum pdlp_solver_mode_t "cuopt::mathematical_optimization::pdlp_solver_mode_t": # noqa
         Stable1 "cuopt::mathematical_optimization::pdlp_solver_mode_t::Stable1" # noqa
@@ -30,6 +34,10 @@ cdef extern from "cuopt/mathematical_optimization/pdlp/solver_settings.hpp" name
         Barrier "cuopt::mathematical_optimization::method_t::Barrier" # noqa
         Unset "cuopt::mathematical_optimization::method_t::Unset" # noqa
 
+    cdef cppclass pdlp_solver_settings_t[i_t, f_t]:
+        bool sequence_solve
+        barrier_cache_t* barrier_cache
+
 cdef extern from "cuopt/mathematical_optimization/solver_settings.hpp" namespace "cuopt::mathematical_optimization": # noqa
 
     cdef cppclass solver_settings_t[i_t, f_t]:
@@ -90,9 +98,12 @@ cdef extern from "cuopt/mathematical_optimization/solver_settings.hpp" namespace
 
         void load_parameters_from_file(const string& path) except +
 
+        pdlp_solver_settings_t[i_t, f_t]& get_pdlp_settings()
+
 
 cdef class SolverSettings:
     cdef unique_ptr[solver_settings_t[int, double]] c_solver_settings
     cdef public dict settings_dict
     cdef public object pdlp_warm_start_data
     cdef public list mip_callbacks
+    cdef public bint sequence_solve
diff --git a/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx b/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx
index a5dcc78d18..d228df88e1 100644
--- a/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx
+++ b/python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.pyx
@@ -1,4 +1,4 @@
-# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # noqa
+# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
 # SPDX-License-Identifier: Apache-2.0
 
 # cython: profile=False
@@ -117,6 +117,7 @@ cdef class SolverSettings:
         self.settings_dict = {}
         self.pdlp_warm_start_data = None
         self.mip_callbacks = []
+        self.sequence_solve = False
 
     def to_base_type(self, value):
         """Convert a string to a base type.
@@ -459,6 +460,8 @@ cdef class SolverSettings:
                 warm_start_data.iterations_since_last_restart # noqa
             )
 
+        c_solver_settings.get_pdlp_settings().sequence_solve = self.sequence_solve
+
     def dump_parameters_to_file(self, path, hyperparameters_only=True):
         """Apply ``settings_dict`` / warm start to C++, then dump parameters to *path*.