Add update apis - RHS

cpp/include/cuopt/mathematical_optimization/utilities/barrier_cache.hpp

@@ -23,6 +23,8 @@ 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);
+
+void apply_barrier_rhs(iteration_data_t<int, double>& data, double const* barrier_b, int m);
 }  // namespace cuopt::mathematical_optimization::barrier
 
 namespace cuopt {
@@ -34,8 +36,8 @@ struct barrier_transform_t;
  * @brief GPU solve cache owned by DataModel when CUOPT_SEQUENCE_SOLVE is enabled.
  *
  * 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).
+ * The update APIs crush new user data into that workspace and mark the cache dirty so the
+ * next Solve reuses it (skip convert/presolve/scaling).
  */
 class barrier_cache_t {
  public:
@@ -64,15 +66,25 @@ class barrier_cache_t {
   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;
+  /** True when an update API has staged new data that the next Solve should reuse. */
+  [[nodiscard]] bool dirty() const;
+  void mark_clean();
+
+  /** True when the last update_rhs made a row presolve dropped as empty infeasible. */
+  [[nodiscard]] bool rhs_infeasible() const;
 
   /**
-   * Crush the input linear objective into cached iteration_data_t.c / d_c_ and set c_dirty.
+   * Crush the input linear objective into cached iteration_data_t.c / d_c_ and mark dirty.
    * Requires a stored transform and iteration_data from an Optimal solve.
    */
   void update_linear_objective(double const* c, int n);
 
+  /**
+   * Crush the input constraint RHS into cached iteration_data_t.b / d_b_ and mark dirty.
+   * Requires a stored transform and iteration_data from an Optimal solve.
+   */
+  void update_rhs(double const* b, int m);
+
  private:
   barrier_cache_t(std::unique_ptr<rmm::cuda_stream> stream, std::unique_ptr<raft::handle_t> handle);
 

cpp/src/barrier/barrier.cu

@@ -895,7 +895,7 @@ class iteration_data_t {
   }
 
   // Attach this solve's settings and rewind iterate-dependent state so barrier can
-  // start with the new c. A and Q are unchanged; the previous solve
+  // start with the new c / b. A and Q are unchanged; the previous solve
   // left D and the KKT values at its last iterate. Reuse is QP-only (no cones),
   // so form_*(false) updates values in the existing CSR; no symbolic rebuild.
   bool reset_iterate_state(const simplex_solver_settings_t<i_t, f_t>& settings)
@@ -4911,6 +4911,18 @@ void apply_barrier_linear_objective(iteration_data_t<int, double>& data,
     data.d_c_.data(), data.c.data(), static_cast<std::size_t>(n), data.handle_ptr->get_stream());
 }
 
+void apply_barrier_rhs(iteration_data_t<int, double>& data, double const* barrier_b, int m)
+{
+  cuopt_expects(
+    barrier_b != nullptr && static_cast<int>(data.b.size()) == m &&
+      static_cast<int>(data.d_b_.size()) == m,
+    error_type_t::ValidationError,
+    "update_rhs: barrier RHS size does not match cached iteration_data_t.");
+  std::copy(barrier_b, barrier_b + m, data.b.data());
+  raft::copy(
+    data.d_b_.data(), data.b.data(), static_cast<std::size_t>(m), 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);

cpp/src/barrier/barrier_cache.cu

@@ -22,6 +22,29 @@ using barrier_iteration_data_t = barrier::iteration_data_t<int, double>;
 using barrier_iteration_data_ptr =
   std::unique_ptr<barrier_iteration_data_t, void (*)(barrier_iteration_data_t*)>;
 
+static void require_warm_cache(barrier_transform_t const* transform,
+                               barrier_iteration_data_t const* data,
+                               char const* api)
+{
+  cuopt_expects(transform != nullptr,
+                error_type_t::ValidationError,
+                "%s: no barrier transform; Solve with CUOPT_SEQUENCE_SOLVE enabled first.",
+                api);
+  cuopt_expects(data != nullptr,
+                error_type_t::ValidationError,
+                "%s: no cached iteration_data; Solve a QP to Optimal first.",
+                api);
+}
+
+// Re-adds the first solve's barrier-minus-crush shift so the update lands in cached coordinates.
+static void add_shift(std::vector<double>& crushed, std::vector<double> const& shift)
+{
+  if (shift.size() != crushed.size()) { return; }
+  for (std::size_t i = 0; i < crushed.size(); ++i) {
+    crushed[i] += shift[i];
+  }
+}
+
 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)),
@@ -36,6 +59,8 @@ struct barrier_cache_t::impl {
   std::unique_ptr<barrier_transform_t> transform;
   barrier_iteration_data_ptr iteration_data;
   bool c_dirty{false};
+  bool b_dirty{false};
+  bool rhs_infeasible{false};
 };
 
 barrier_cache_t::barrier_cache_t(std::unique_ptr<rmm::cuda_stream> stream,
@@ -66,7 +91,7 @@ void barrier_cache_t::clear()
 {
   impl_->iteration_data.reset();
   impl_->transform.reset();
-  impl_->c_dirty = false;
+  mark_clean();
 }
 
 void barrier_cache_t::store_iteration_data(barrier_iteration_data_t* data)
@@ -88,22 +113,25 @@ 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::dirty() const
+{
+  return (impl_->c_dirty || impl_->b_dirty) && impl_->transform != nullptr &&
+         impl_->iteration_data.get() != nullptr;
+}
 
-bool barrier_cache_t::c_dirty() const
+void barrier_cache_t::mark_clean()
 {
-  return impl_->c_dirty && impl_->transform != nullptr && impl_->iteration_data.get() != nullptr;
+  impl_->c_dirty        = false;
+  impl_->b_dirty        = false;
+  impl_->rhs_infeasible = false;
 }
 
+bool barrier_cache_t::rhs_infeasible() const { return impl_->rhs_infeasible; }
+
 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 CUOPT_SEQUENCE_SOLVE "
-                "enabled 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.");
+  require_warm_cache(
+    impl_->transform.get(), impl_->iteration_data.get(), "update_linear_objective");
   // Cached Q and c are in minimization space.
   std::vector<double> user_objective;
   if (impl_->transform->maximize && c != nullptr && n > 0) {
@@ -135,11 +163,7 @@ void barrier_cache_t::update_linear_objective(double const* c, int n)
     }
     barrier_lp.obj_constant += obj_constant_delta;
   }
-  if (linear_obj_shift.size() == crushed.size()) {
-    for (std::size_t j = 0; j < crushed.size(); ++j) {
-      crushed[j] += linear_obj_shift[j];
-    }
-  }
+  add_shift(crushed, linear_obj_shift);
   // The next solve builds its solver from barrier_lp, so keep its objective and the cached
   // iteration workspace on the same c.
   auto& barrier_objective = barrier_lp.objective;
@@ -153,4 +177,33 @@ void barrier_cache_t::update_linear_objective(double const* c, int n)
   impl_->c_dirty = true;
 }
 
+void barrier_cache_t::update_rhs(double const* b, int m)
+{
+  require_warm_cache(impl_->transform.get(), impl_->iteration_data.get(), "update_rhs");
+  std::vector<double> crushed;
+  try {
+    crushed = crush_user_rhs(*impl_->transform, b, m);
+  } catch (update_rhs_infeasible_error const&) {
+    // Cache stays usable for a later feasible update; the next Solve reports INFEASIBLE from
+    // this flag without running IPM.
+    impl_->rhs_infeasible = true;
+    impl_->b_dirty        = true;
+    return;
+  } catch (std::invalid_argument const& e) {
+    cuopt_expects(false, error_type_t::ValidationError, "%s", e.what());
+  }
+  impl_->rhs_infeasible = false;
+  add_shift(crushed, impl_->transform->rhs_shift);
+  // barrier_lp->rhs also seeds the next solve's Mehrotra start, so keep it and the cached
+  // workspace on the same b.
+  auto& barrier_rhs = impl_->transform->barrier_lp->rhs;
+  cuopt_expects(barrier_rhs.size() == crushed.size(),
+                error_type_t::ValidationError,
+                "update_rhs: crushed RHS size does not match the cached barrier LP.");
+  barrier_rhs = crushed;
+  barrier::apply_barrier_rhs(
+    *impl_->iteration_data, crushed.data(), static_cast<int>(crushed.size()));
+  impl_->b_dirty = true;
+}
+
 }  // namespace cuopt::mathematical_optimization

cpp/src/barrier/barrier_transform.hpp

@@ -10,17 +10,19 @@
 #include <dual_simplex/presolve.hpp>
 #include <linear_algebra/sparse_matrix.hpp>
 
+#include <cmath>
 #include <memory>
 #include <stdexcept>
+#include <string>
 #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.
+ * Enough to crush new linear objective or RHS data 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};
@@ -43,6 +45,12 @@ struct barrier_transform_t {
   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;
+  // Barrier RHS minus crush(user b) from the first solve (fixed/lower-bound shifts).
+  std::vector<double> rhs_shift;
+  // False when range rows or folding put the user RHS somewhere other than barrier_lp->rhs.
+  bool rhs_update_supported{false};
+  // Absolute primal tolerance of the first solve, used to test rows presolve dropped as empty.
+  double primal_tol{1e-6};
   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.
@@ -109,4 +117,77 @@ inline std::vector<double> crush_user_linear_objective(barrier_transform_t const
   return presolved;
 }
 
+// Distinct from the invalid_argument cases so the caller can report INFEASIBLE rather than a
+// validation failure.
+struct update_rhs_infeasible_error : std::runtime_error {
+  explicit update_rhs_infeasible_error(std::string const& message) : std::runtime_error(message) {}
+};
+
+inline std::vector<double> crush_user_rhs(barrier_transform_t const& xf, double const* b, int m)
+{
+  if (b == nullptr || m != xf.user_num_rows) {
+    throw std::invalid_argument("update_rhs: RHS length must match the cached user row count.");
+  }
+  if (!xf.rhs_update_supported) {
+    throw std::invalid_argument(
+      "update_rhs: cached convert used range rows or folding; run a full Solve.");
+  }
+  if (xf.original_num_rows != xf.user_num_rows) {
+    throw std::invalid_argument(
+      "update_rhs: cached original row count does not match the user row count.");
+  }
+  if (static_cast<int>(xf.row_sense.size()) != xf.user_num_rows) {
+    throw std::invalid_argument(
+      "update_rhs: cached row-sense count does not match the user row count.");
+  }
+  if (xf.barrier_lp == nullptr) {
+    throw std::invalid_argument("update_rhs: cached barrier LP is missing.");
+  }
+
+  // convert turns 'G' rows into 'L' rows by negating the row and its RHS.
+  std::vector<double> original(static_cast<std::size_t>(xf.original_num_rows));
+  for (int i = 0; i < m; ++i) {
+    original[static_cast<std::size_t>(i)] =
+      xf.row_sense[static_cast<std::size_t>(i)] == 'G' ? -b[i] : b[i];
+  }
+
+  // Dropped rows were empty, so the new RHS never reaches the barrier: 'E' needs 0 == b_i and
+  // the rest need 0 <= b_i.
+  for (int i : xf.presolve_info.removed_constraints) {
+    if (i < 0 || i >= m) {
+      throw std::invalid_argument("update_rhs: removed constraint index is out of range.");
+    }
+    double const converted_rhs = original[static_cast<std::size_t>(i)];
+    bool const infeasible      = xf.row_sense[static_cast<std::size_t>(i)] == 'E'
+                                   ? std::abs(converted_rhs) > xf.primal_tol
+                                   : converted_rhs < -xf.primal_tol;
+    if (infeasible) {
+      throw update_rhs_infeasible_error("update_rhs: empty constraint row " + std::to_string(i) +
+                                        " is infeasible with the new RHS.");
+    }
+  }
+
+  // Empty remaining_constraints means either no empty-row pass ran, or every row was dropped
+  // and accepted above.
+  std::vector<double> presolved;
+  if (!xf.presolve_info.remaining_constraints.empty()) {
+    presolved.resize(xf.presolve_info.remaining_constraints.size());
+    for (std::size_t k = 0; k < xf.presolve_info.remaining_constraints.size(); ++k) {
+      presolved[k] = original[static_cast<std::size_t>(xf.presolve_info.remaining_constraints[k])];
+    }
+  } else if (xf.presolve_info.removed_constraints.empty()) {
+    presolved = std::move(original);
+  }
+
+  if (static_cast<int>(presolved.size()) != xf.barrier_lp->num_rows ||
+      xf.row_scales.size() != presolved.size()) {
+    throw std::invalid_argument(
+      "update_rhs: crushed RHS size does not match barrier rows / row_scales.");
+  }
+  for (std::size_t i = 0; i < presolved.size(); ++i) {
+    presolved[i] /= xf.row_scales[i];
+  }
+  return presolved;
+}
+
 }  // namespace cuopt::mathematical_optimization