diff --git a/ortools/base/parse_text_proto.h b/ortools/base/parse_text_proto.h
index 625ad29ad1..5e78d3dfc3 100644
--- a/ortools/base/parse_text_proto.h
+++ b/ortools/base/parse_text_proto.h
@@ -14,6 +14,8 @@
 #ifndef OR_TOOLS_BASE_PARSE_TEXT_PROTO_H_
 #define OR_TOOLS_BASE_PARSE_TEXT_PROTO_H_
 
+#include <string_view>
+
 #include "absl/log/absl_check.h"
 #include "google/protobuf/message.h"
 #include "google/protobuf/text_format.h"
@@ -32,6 +34,31 @@ T ParseTextOrDie(const std::string& input) {
   return result;
 }
 
+namespace text_proto_internal {
+
+class ParseProtoHelper {
+ public:
+  explicit ParseProtoHelper(std::string_view asciipb) : asciipb_(asciipb) {}
+  template <class T>
+  operator T() {  // NOLINT(runtime/explicit)
+    T result;
+    const bool ok = ::google::protobuf::TextFormat::TextFormat::ParseFromString(
+        asciipb_, &result);
+    CHECK(ok) << "Failed to parse text proto: " << asciipb_;
+    return result;
+  }
+
+ private:
+  const std::string asciipb_;
+};
+
+}  // namespace text_proto_internal
+
+text_proto_internal::ParseProtoHelper ParseTextProtoOrDie(
+    std::string_view input) {
+  return text_proto_internal::ParseProtoHelper(input);
+}
+
 }  // namespace google::protobuf::contrib::parse_proto
 
 #endif  // OR_TOOLS_BASE_PARSE_TEXT_PROTO_H_
diff --git a/ortools/glop/revised_simplex.cc b/ortools/glop/revised_simplex.cc
index 88549afa30..2a7e20dc8a 100644
--- a/ortools/glop/revised_simplex.cc
+++ b/ortools/glop/revised_simplex.cc
@@ -154,29 +154,72 @@ void RevisedSimplex::SetStartingVariableValuesForNextSolve(
   variable_starting_values_ = values;
 }
 
-void RevisedSimplex::NotifyThatMatrixIsUnchangedForNextSolve() {
-  notify_that_matrix_is_unchanged_ = true;
-}
+Status RevisedSimplex::MinimizeFromTransposedMatrixWithSlack(
+    const DenseRow& objective, Fractional objective_scaling_factor,
+    Fractional objective_offset, TimeLimit* time_limit) {
+  const double start_time = time_limit->GetElapsedTime();
+  default_logger_.EnableLogging(parameters_.log_search_progress());
+  default_logger_.SetLogToStdOut(parameters_.log_to_stdout());
+  parameters_ = initial_parameters_;
+  PropagateParameters();
+
+  // The source of truth is the transposed matrix.
+  if (transpose_was_changed_) {
+    compact_matrix_.PopulateFromTranspose(transposed_matrix_);
+    num_rows_ = compact_matrix_.num_rows();
+    num_cols_ = compact_matrix_.num_cols();
+    first_slack_col_ = num_cols_ - RowToColIndex(num_rows_);
+  }
+
+  DCHECK_EQ(num_cols_, objective.size());
+
+  // Copy objective
+  objective_scaling_factor_ = objective_scaling_factor;
+  objective_offset_ = objective_offset;
+  const bool objective_is_unchanged = objective_ == objective;
+  objective_ = objective;
+  InitializeObjectiveLimit();
+
+  // Initialize variable infos from the mutated bounds.
+  variables_info_.InitializeFromMutatedState();
+
+  if (objective_is_unchanged && parameters_.use_dual_simplex() &&
+      !transpose_was_changed_ && !solution_state_has_been_set_externally_ &&
+      !solution_state_.IsEmpty()) {
+    // Fast track if we just changed variable bounds.
+    primal_edge_norms_.Clear();
+    variables_info_.InitializeFromBasisState(first_slack_col_, ColIndex(0),
+                                             solution_state_);
+    variable_values_.ResetAllNonBasicVariableValues(variable_starting_values_);
+    variable_values_.RecomputeBasicVariableValues();
+    return SolveInternal(start_time, false, objective, time_limit);
+  } else {
+    GLOP_RETURN_IF_ERROR(FinishInitialization(true));
+  }
 
-void RevisedSimplex::NotifyThatMatrixIsChangedForNextSolve() {
-  notify_that_matrix_is_unchanged_ = false;
+  return SolveInternal(start_time, false, objective, time_limit);
 }
 
 Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) {
-  SCOPED_TIME_STAT(&function_stats_);
+  const double start_time = time_limit->GetElapsedTime();
+  default_logger_.EnableLogging(parameters_.log_search_progress());
+  default_logger_.SetLogToStdOut(parameters_.log_to_stdout());
+
   DCHECK(lp.IsCleanedUp());
+  GLOP_RETURN_IF_ERROR(Initialize(lp));
+  return SolveInternal(start_time, lp.IsMaximizationProblem(),
+                       lp.objective_coefficients(), time_limit);
+}
+
+ABSL_MUST_USE_RESULT Status RevisedSimplex::SolveInternal(
+    double start_time, bool is_maximization_problem,
+    const DenseRow& objective_coefficients, TimeLimit* time_limit) {
+  SCOPED_TIME_STAT(&function_stats_);
   GLOP_RETURN_ERROR_IF_NULL(time_limit);
   Cleanup update_deterministic_time_on_return(
       [this, time_limit]() { AdvanceDeterministicTime(time_limit); });
 
-  default_logger_.EnableLogging(parameters_.log_search_progress());
-  default_logger_.SetLogToStdOut(parameters_.log_to_stdout());
   SOLVER_LOG(logger_, "");
-
-  // Initialization. Note That Initialize() must be called first since it
-  // analyzes the current solver state.
-  const double start_time = time_limit->GetElapsedTime();
-  GLOP_RETURN_IF_ERROR(Initialize(lp));
   if (logger_->LoggingIsEnabled()) {
     DisplayBasicVariableStatistics();
   }
@@ -310,7 +353,13 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) {
 
     // After the primal phase I, we need to restore the objective.
     if (problem_status_ != ProblemStatus::PRIMAL_INFEASIBLE) {
-      InitializeObjectiveAndTestIfUnchanged(lp);
+      objective_ = objective_coefficients;
+      if (is_maximization_problem) {
+        for (Fractional& value : objective_) {
+          value = -value;
+        }
+      }
+      objective_.resize(num_cols_, 0.0);  // For the slack.
       reduced_costs_.ResetForNewObjective();
     }
   }
@@ -639,7 +688,7 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) {
   solution_reduced_costs_ = reduced_costs_.GetReducedCosts();
   SaveState();
 
-  if (lp.IsMaximizationProblem()) {
+  if (is_maximization_problem) {
     ChangeSign(&solution_dual_values_);
     ChangeSign(&solution_reduced_costs_);
   }
@@ -650,7 +699,7 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) {
     solution_objective_value_ =
         (problem_status_ == ProblemStatus::DUAL_UNBOUNDED) ? kInfinity
                                                            : -kInfinity;
-    if (lp.IsMaximizationProblem()) {
+    if (is_maximization_problem) {
       solution_objective_value_ = -solution_objective_value_;
     }
   }
@@ -1379,21 +1428,13 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) {
   ColIndex num_new_cols(0);
   bool only_change_is_new_rows = false;
   bool only_change_is_new_cols = false;
-  bool matrix_is_unchanged = true;
-  bool only_new_bounds = false;
-  if (solution_state_.IsEmpty() || !notify_that_matrix_is_unchanged_) {
-    matrix_is_unchanged = InitializeMatrixAndTestIfUnchanged(
-        lp, lp_is_in_equation_form, &only_change_is_new_rows,
-        &only_change_is_new_cols, &num_new_cols);
-    only_new_bounds = only_change_is_new_cols && num_new_cols > 0 &&
-                      OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero(
-                          lp, lp_is_in_equation_form, num_new_cols);
-  } else if (DEBUG_MODE) {
-    CHECK(InitializeMatrixAndTestIfUnchanged(
-        lp, lp_is_in_equation_form, &only_change_is_new_rows,
-        &only_change_is_new_cols, &num_new_cols));
-  }
-  notify_that_matrix_is_unchanged_ = false;
+  const bool matrix_is_unchanged = InitializeMatrixAndTestIfUnchanged(
+      lp, lp_is_in_equation_form, &only_change_is_new_rows,
+      &only_change_is_new_cols, &num_new_cols);
+  const bool only_new_bounds =
+      only_change_is_new_cols && num_new_cols > 0 &&
+      OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero(
+          lp, lp_is_in_equation_form, num_new_cols);
 
   // TODO(user): move objective with ReducedCosts class.
   const bool objective_is_unchanged = InitializeObjectiveAndTestIfUnchanged(lp);
@@ -1509,6 +1550,10 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) {
     }
   }
 
+  return FinishInitialization(solve_from_scratch);
+}
+
+Status RevisedSimplex::FinishInitialization(bool solve_from_scratch) {
   // If we couldn't perform a "quick" warm start above, we can at least try to
   // reuse the variable statuses.
   if (solve_from_scratch && !solution_state_.IsEmpty()) {
@@ -1589,6 +1634,8 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) {
     SOLVER_LOG(logger_, "Starting basis: incremental solve.");
   }
   DCHECK(BasisIsConsistent());
+
+  transpose_was_changed_ = false;
   return Status::OK();
 }
 
diff --git a/ortools/glop/revised_simplex.h b/ortools/glop/revised_simplex.h
index a4fa4b51ad..90e55698e4 100644
--- a/ortools/glop/revised_simplex.h
+++ b/ortools/glop/revised_simplex.h
@@ -170,14 +170,6 @@ class RevisedSimplex {
   // variables.
   void SetStartingVariableValuesForNextSolve(const DenseRow& values);
 
-  // Advanced usage. Tells the next Solve() that the matrix inside the linear
-  // program will not change compared to the one used the last time Solve() was
-  // called. This allows to bypass the somewhat costly check of comparing both
-  // matrices. Note that this call will be ignored if Solve() was never called
-  // or if ClearStateForNextSolve() was called.
-  void NotifyThatMatrixIsUnchangedForNextSolve();
-  void NotifyThatMatrixIsChangedForNextSolve();
-
   // Getters to retrieve all the information computed by the last Solve().
   RowIndex GetProblemNumRows() const;
   ColIndex GetProblemNumCols() const;
@@ -252,6 +244,24 @@ class RevisedSimplex {
 
   void SetLogger(SolverLogger* logger) { logger_ = logger; }
 
+  // Advanced usage. For fast incremental call to the solver, it is better not
+  // to use LinearProgram at all. This api allows to directly modify the
+  // internal data of glop and then call solve.
+  const CompactSparseMatrix& MatrixWithSlack() const { return compact_matrix_; }
+  CompactSparseMatrix* MutableTransposedMatrixWithSlack() {
+    transpose_was_changed_ = true;
+    return &transposed_matrix_;
+  }
+  DenseRow* MutableLowerBounds() {
+    return variables_info_.MutableLowerBounds();
+  }
+  DenseRow* MutableUpperBounds() {
+    return variables_info_.MutableUpperBounds();
+  }
+  ABSL_MUST_USE_RESULT Status MinimizeFromTransposedMatrixWithSlack(
+      const DenseRow& objective, Fractional objective_scaling_factor,
+      Fractional objective_offset, TimeLimit* time_limit);
+
  private:
   struct IterationStats : public StatsGroup {
     IterationStats()
@@ -303,6 +313,10 @@ class RevisedSimplex {
     FINAL_CHECK
   };
 
+  ABSL_MUST_USE_RESULT Status SolveInternal(double start_time, bool maximize,
+                                            const DenseRow& objective,
+                                            TimeLimit* time_limit);
+
   // Propagates parameters_ to all the other classes that need it.
   //
   // TODO(user): Maybe a better design is for them to have a reference to a
@@ -427,6 +441,7 @@ class RevisedSimplex {
 
   // Entry point for the solver initialization.
   ABSL_MUST_USE_RESULT Status Initialize(const LinearProgram& lp);
+  ABSL_MUST_USE_RESULT Status FinishInitialization(bool solve_from_scratch);
 
   // Saves the current variable statuses in solution_state_.
   void SaveState();
@@ -715,9 +730,8 @@ class RevisedSimplex {
   // If this is cleared, we assume they are none.
   DenseRow variable_starting_values_;
 
-  // Flag used by NotifyThatMatrixIsUnchangedForNextSolve() and changing
-  // the behavior of Initialize().
-  bool notify_that_matrix_is_unchanged_ = false;
+  // See MutableTransposedMatrixWithSlack().
+  bool transpose_was_changed_ = false;
 
   // This is known as 'd' in the literature and is set during each pivot to the
   // right inverse of the basic entering column of A by ComputeDirection().
diff --git a/ortools/glop/variables_info.cc b/ortools/glop/variables_info.cc
index d100f15eea..dd5b0d8d71 100644
--- a/ortools/glop/variables_info.cc
+++ b/ortools/glop/variables_info.cc
@@ -46,6 +46,16 @@ bool VariablesInfo::LoadBoundsAndReturnTrueIfUnchanged(
   return false;
 }
 
+void VariablesInfo::InitializeFromMutatedState() {
+  const ColIndex num_cols = matrix_.num_cols();
+  DCHECK_EQ(num_cols, lower_bounds_.size());
+  DCHECK_EQ(num_cols, upper_bounds_.size());
+  variable_type_.resize(num_cols, VariableType::UNCONSTRAINED);
+  for (ColIndex col(0); col < num_cols; ++col) {
+    variable_type_[col] = ComputeVariableType(col);
+  }
+}
+
 bool VariablesInfo::LoadBoundsAndReturnTrueIfUnchanged(
     const DenseRow& variable_lower_bounds,
     const DenseRow& variable_upper_bounds,
diff --git a/ortools/glop/variables_info.h b/ortools/glop/variables_info.h
index 2a8fa4f7e6..944d39cfc4 100644
--- a/ortools/glop/variables_info.h
+++ b/ortools/glop/variables_info.h
@@ -174,6 +174,12 @@ class VariablesInfo {
   void EndDualPhaseI(Fractional dual_feasibility_tolerance,
                      DenseRow::ConstView reduced_costs);
 
+  // Advanced incremental API to reuse directly the internal storage.
+  // This saves two copy per solves, and only matter on large easy problems.
+  void InitializeFromMutatedState();
+  DenseRow* MutableLowerBounds() { return &lower_bounds_; }
+  DenseRow* MutableUpperBounds() { return &upper_bounds_; }
+
  private:
   // Computes the initial/default variable status from its type. A constrained
   // variable is set to the lowest of its 2 bounds in absolute value.
diff --git a/ortools/lp_data/lp_data_utils.cc b/ortools/lp_data/lp_data_utils.cc
index 25d03a176a..126121e65b 100644
--- a/ortools/lp_data/lp_data_utils.cc
+++ b/ortools/lp_data/lp_data_utils.cc
@@ -13,7 +13,12 @@
 
 #include "ortools/lp_data/lp_data_utils.h"
 
+#include <algorithm>
+#include <cstdlib>
+#include <limits>
+
 #include "absl/log/check.h"
+#include "absl/types/span.h"
 #include "ortools/glop/parameters.pb.h"
 #include "ortools/lp_data/lp_data.h"
 #include "ortools/lp_data/lp_types.h"
@@ -84,14 +89,40 @@ void Scale(LinearProgram* lp, SparseMatrixScaler* scaler,
 void LpScalingHelper::Scale(LinearProgram* lp) { Scale(GlopParameters(), lp); }
 
 void LpScalingHelper::Scale(const GlopParameters& params, LinearProgram* lp) {
-  scaler_.Clear();
-  ::operations_research::glop::Scale(lp, &scaler_, params.scaling_method());
+  SparseMatrixScaler scaler;
+  ::operations_research::glop::Scale(lp, &scaler, params.scaling_method());
   bound_scaling_factor_ = 1.0 / lp->ScaleBounds();
   objective_scaling_factor_ = 1.0 / lp->ScaleObjective(params.cost_scaling());
+
+  matrix_is_scaled_ = true;
+  row_unscaling_factors_ = scaler.row_scales();
+  col_unscaling_factors_ = scaler.col_scales();
+
+  // It is possible the scaler didn't do anything.
+  // we still allocate the vector though since we don't test that below.
+  row_unscaling_factors_.resize(lp->num_constraints(), 1.0);
+  col_unscaling_factors_.resize(lp->num_variables(), 1.0);
+}
+
+void LpScalingHelper::ConfigureFromFactors(
+    absl::Span<const double> row_factors,
+    absl::Span<const double> col_factors) {
+  matrix_is_scaled_ = true;
+  const RowIndex num_rows(row_factors.size());
+  row_unscaling_factors_.resize(num_rows, 1.0);
+  for (RowIndex row(0); row < num_rows; ++row) {
+    row_unscaling_factors_[row] = 1.0 / row_factors[row.value()];
+  }
+
+  const ColIndex num_cols(col_factors.size());
+  col_unscaling_factors_.resize(num_cols, 1.0);
+  for (ColIndex col(0); col < num_cols; ++col) {
+    col_unscaling_factors_[col] = 1.0 / col_factors[col.value()];
+  }
 }
 
 void LpScalingHelper::Clear() {
-  scaler_.Clear();
+  matrix_is_scaled_ = false;
   bound_scaling_factor_ = 1.0;
   objective_scaling_factor_ = 1.0;
 }
@@ -99,71 +130,81 @@ void LpScalingHelper::Clear() {
 Fractional LpScalingHelper::VariableScalingFactor(ColIndex col) const {
   // During scaling a col was multiplied by ColScalingFactor() and the variable
   // bounds divided by it.
-  return scaler_.ColUnscalingFactor(col) * bound_scaling_factor_;
+  return ColUnscalingFactor(col) * bound_scaling_factor_;
+}
+
+Fractional LpScalingHelper::VariableScalingFactorWithSlack(ColIndex col) const {
+  if (!matrix_is_scaled_) return bound_scaling_factor_;
+  const ColIndex num_cols = col_unscaling_factors_.size();
+  if (col < num_cols) {
+    return col_unscaling_factors_[col] * bound_scaling_factor_;
+  }
+  return row_unscaling_factors_[ColToRowIndex(col - num_cols)] *
+         bound_scaling_factor_;
 }
 
 Fractional LpScalingHelper::ScaleVariableValue(ColIndex col,
                                                Fractional value) const {
-  return value * scaler_.ColUnscalingFactor(col) * bound_scaling_factor_;
+  return value * ColUnscalingFactor(col) * bound_scaling_factor_;
 }
 
 Fractional LpScalingHelper::ScaleReducedCost(ColIndex col,
                                              Fractional value) const {
   // The reduced cost move like the objective and the col scale.
-  return value / scaler_.ColUnscalingFactor(col) * objective_scaling_factor_;
+  return value / ColUnscalingFactor(col) * objective_scaling_factor_;
 }
 
 Fractional LpScalingHelper::ScaleDualValue(RowIndex row,
                                            Fractional value) const {
   // The dual value move like the objective and the inverse of the row scale.
-  return value * (scaler_.RowUnscalingFactor(row) * objective_scaling_factor_);
+  return value * (RowUnscalingFactor(row) * objective_scaling_factor_);
 }
 
 Fractional LpScalingHelper::ScaleConstraintActivity(RowIndex row,
                                                     Fractional value) const {
   // The activity move with the row_scale and the bound_scaling_factor.
-  return value / scaler_.RowUnscalingFactor(row) * bound_scaling_factor_;
+  return value / RowUnscalingFactor(row) * bound_scaling_factor_;
 }
 
 Fractional LpScalingHelper::UnscaleVariableValue(ColIndex col,
                                                  Fractional value) const {
   // Just the opposite of ScaleVariableValue().
-  return value / (scaler_.ColUnscalingFactor(col) * bound_scaling_factor_);
+  return value / (ColUnscalingFactor(col) * bound_scaling_factor_);
 }
 
 Fractional LpScalingHelper::UnscaleReducedCost(ColIndex col,
                                                Fractional value) const {
   // The reduced cost move like the objective and the col scale.
-  return value * scaler_.ColUnscalingFactor(col) / objective_scaling_factor_;
+  return value * ColUnscalingFactor(col) / objective_scaling_factor_;
 }
 
 Fractional LpScalingHelper::UnscaleDualValue(RowIndex row,
                                              Fractional value) const {
   // The dual value move like the objective and the inverse of the row scale.
-  return value / (scaler_.RowUnscalingFactor(row) * objective_scaling_factor_);
+  return value / (RowUnscalingFactor(row) * objective_scaling_factor_);
 }
 
 Fractional LpScalingHelper::UnscaleConstraintActivity(RowIndex row,
                                                       Fractional value) const {
   // The activity move with the row_scale and the bound_scaling_factor.
-  return value * scaler_.RowUnscalingFactor(row) / bound_scaling_factor_;
+  return value * RowUnscalingFactor(row) / bound_scaling_factor_;
 }
 
 void LpScalingHelper::UnscaleUnitRowLeftSolve(
     ColIndex basis_col, ScatteredRow* left_inverse) const {
-  const Fractional global_factor = scaler_.ColUnscalingFactor(basis_col);
+  const Fractional global_factor = ColUnscalingFactor(basis_col);
 
   // We have left_inverse * [RowScale * B * ColScale] = unit_row.
   if (left_inverse->non_zeros.empty()) {
     const ColIndex num_rows = left_inverse->values.size();
     for (ColIndex col(0); col < num_rows; ++col) {
       left_inverse->values[col] /=
-          scaler_.RowUnscalingFactor(ColToRowIndex(col)) * global_factor;
+          RowUnscalingFactor(ColToRowIndex(col)) * global_factor;
     }
   } else {
     for (const ColIndex col : left_inverse->non_zeros) {
       left_inverse->values[col] /=
-          scaler_.RowUnscalingFactor(ColToRowIndex(col)) * global_factor;
+          RowUnscalingFactor(ColToRowIndex(col)) * global_factor;
     }
   }
 }
@@ -171,7 +212,7 @@ void LpScalingHelper::UnscaleUnitRowLeftSolve(
 void LpScalingHelper::UnscaleColumnRightSolve(
     const RowToColMapping& basis, ColIndex col,
     ScatteredColumn* right_inverse) const {
-  const Fractional global_factor = scaler_.ColScalingFactor(col);
+  const Fractional global_factor = 1.0 / ColUnscalingFactor(col);
 
   // [RowScale * B * BColScale] * inverse = RowScale * column * ColScale.
   // That is B * (BColScale * inverse) = column * ColScale[col].
@@ -179,14 +220,71 @@ void LpScalingHelper::UnscaleColumnRightSolve(
     const RowIndex num_rows = right_inverse->values.size();
     for (RowIndex row(0); row < num_rows; ++row) {
       right_inverse->values[row] /=
-          scaler_.ColUnscalingFactor(basis[row]) * global_factor;
+          ColUnscalingFactor(basis[row]) * global_factor;
     }
   } else {
     for (const RowIndex row : right_inverse->non_zeros) {
       right_inverse->values[row] /=
-          scaler_.ColUnscalingFactor(basis[row]) * global_factor;
+          ColUnscalingFactor(basis[row]) * global_factor;
+    }
+  }
+}
+
+void LpScalingHelper::AverageCostScaling(DenseRow* objective) {
+  Fractional sum = 0.0;
+  int num_terms = 0;
+  for (const Fractional f : *objective) {
+    if (f == 0) continue;
+    ++num_terms;
+    sum += std::abs(f);
+  }
+  if (num_terms == 0) {
+    objective_scaling_factor_ = 1.0;
+    return;
+  }
+
+  const Fractional average = sum / static_cast<double>(num_terms);
+  objective_scaling_factor_ = 1.0 / average;
+  for (Fractional& f : *objective) {
+    f *= objective_scaling_factor_;
+  }
+}
+
+void LpScalingHelper::ContainOneBoundScaling(DenseRow* upper_bounds,
+                                             DenseRow* lower_bounds) {
+  const double infinity = std::numeric_limits<double>::infinity();
+  Fractional min_magnitude = infinity;
+  Fractional max_magnitude = 0.0;
+  for (const Fractional f : *lower_bounds) {
+    const Fractional m = std::abs(f);
+    if (m == 0 || m == infinity) continue;
+    min_magnitude = std::min(min_magnitude, m);
+    max_magnitude = std::max(max_magnitude, m);
+  }
+  for (const Fractional f : *upper_bounds) {
+    const Fractional m = std::abs(f);
+    if (m == 0 || m == infinity) continue;
+    min_magnitude = std::min(min_magnitude, m);
+    max_magnitude = std::max(max_magnitude, m);
+  }
+
+  bound_scaling_factor_ = 1.0;
+  if (min_magnitude != infinity) {
+    CHECK_LE(min_magnitude, max_magnitude);
+    if (min_magnitude > 1.0) {
+      bound_scaling_factor_ = 1.0 / min_magnitude;
+    } else if (max_magnitude < 1.0) {
+      bound_scaling_factor_ = 1.0 / max_magnitude;
     }
   }
+
+  if (bound_scaling_factor_ == 1.0) return;
+  for (Fractional& f : *lower_bounds) {
+    f *= bound_scaling_factor_;
+  }
+  for (Fractional& f : *upper_bounds) {
+    f *= bound_scaling_factor_;
+  }
 }
 
 }  // namespace glop
diff --git a/ortools/lp_data/lp_data_utils.h b/ortools/lp_data/lp_data_utils.h
index d63eb52a3e..37e94249ce 100644
--- a/ortools/lp_data/lp_data_utils.h
+++ b/ortools/lp_data/lp_data_utils.h
@@ -51,12 +51,14 @@ void Scale(LinearProgram* lp, SparseMatrixScaler* scaler);
 // sense to have a single place where all the scaling formulas are kept.
 class LpScalingHelper {
  public:
+  // Clear all scaling coefficients.
+  void Clear();
+
   // Scale the given LP.
   void Scale(LinearProgram* lp);
   void Scale(const GlopParameters& params, LinearProgram* lp);
-
-  // Clear all scaling coefficients.
-  void Clear();
+  void ConfigureFromFactors(absl::Span<const double> row_factors,
+                            absl::Span<const double> col_factors);
 
   // Transforms value from unscaled domain to the scaled one.
   Fractional ScaleVariableValue(ColIndex col, Fractional value) const;
@@ -83,18 +85,37 @@ class LpScalingHelper {
   // to be in the scaled domain.
   Fractional VariableScalingFactor(ColIndex col) const;
 
-  // Visible for testing. All objective coefficients of the original LP where
-  // multiplied by this factor. Nothing else changed.
-  Fractional BoundsScalingFactor() const { return bound_scaling_factor_; }
+  // Same as VariableScalingFactor() except that ColIndex greater than the
+  // number of columns will be interpreted as "slack" variable whose scaling
+  // factor depends on the row.
+  Fractional VariableScalingFactorWithSlack(ColIndex col) const;
+
+  // Extra scaling function, to scale objective/bounds.
+  void AverageCostScaling(DenseRow* objective);
+  void ContainOneBoundScaling(DenseRow* upper_bounds, DenseRow* lower_bounds);
 
   // Visible for testing. All variable/constraint bounds of the original LP
   // where multiplied by this factor. Nothing else changed.
+  Fractional BoundsScalingFactor() const { return bound_scaling_factor_; }
+
+  // Visible for testing. All objective coefficients of the original LP where
+  // multiplied by this factor. Nothing else changed.
   Fractional ObjectiveScalingFactor() const {
     return objective_scaling_factor_;
   }
 
  private:
-  SparseMatrixScaler scaler_;
+  Fractional RowUnscalingFactor(RowIndex row) const {
+    return matrix_is_scaled_ ? row_unscaling_factors_[row] : 1.0;
+  }
+  Fractional ColUnscalingFactor(ColIndex col) const {
+    return matrix_is_scaled_ ? col_unscaling_factors_[col] : 1.0;
+  }
+
+  bool matrix_is_scaled_ = false;
+  DenseColumn row_unscaling_factors_;
+  DenseRow col_unscaling_factors_;
+
   Fractional bound_scaling_factor_ = 1.0;
   Fractional objective_scaling_factor_ = 1.0;
 };
diff --git a/ortools/lp_data/sparse.cc b/ortools/lp_data/sparse.cc
index 1928aa7d55..585744ecd9 100644
--- a/ortools/lp_data/sparse.cc
+++ b/ortools/lp_data/sparse.cc
@@ -577,6 +577,17 @@ void TriangularMatrix::Reset(RowIndex num_rows, ColIndex col_capacity) {
   starts_[ColIndex(0)] = 0;
 }
 
+void CompactSparseMatrix::AddEntryToCurrentColumn(RowIndex row,
+                                                  Fractional coeff) {
+  rows_.push_back(row);
+  coefficients_.push_back(coeff);
+}
+
+void CompactSparseMatrix::CloseCurrentColumn() {
+  starts_.push_back(rows_.size());
+  ++num_cols_;
+}
+
 ColIndex CompactSparseMatrix::AddDenseColumn(const DenseColumn& dense_column) {
   return AddDenseColumnPrefix(dense_column.const_view(), RowIndex(0));
 }
@@ -832,6 +843,7 @@ void TriangularMatrix::UpperSolveInternal(DenseColumn::View rhs) const {
   const auto entry_rows = rows_.view();
   const auto entry_coefficients = coefficients_.view();
   const auto diagonal_coefficients = diagonal_coefficients_.view();
+  const auto starts = starts_.view();
   for (ColIndex col(diagonal_coefficients.size() - 1); col >= end; --col) {
     const Fractional value = rhs[ColToRowIndex(col)];
     if (value == 0.0) continue;
@@ -844,8 +856,8 @@ void TriangularMatrix::UpperSolveInternal(DenseColumn::View rhs) const {
     // It is faster to iterate this way (instead of i : Column(col)) because of
     // cache locality. Note that the floating-point computations are exactly the
     // same in both cases.
-    const EntryIndex i_end = starts_[col];
-    for (EntryIndex i(starts_[col + 1] - 1); i >= i_end; --i) {
+    const EntryIndex i_end = starts[col];
+    for (EntryIndex i(starts[col + 1] - 1); i >= i_end; --i) {
       rhs[entry_rows[i]] -= coeff * entry_coefficients[i];
     }
   }
diff --git a/ortools/lp_data/sparse.h b/ortools/lp_data/sparse.h
index 96c4571444..19053dba32 100644
--- a/ortools/lp_data/sparse.h
+++ b/ortools/lp_data/sparse.h
@@ -369,6 +369,10 @@ class CompactSparseMatrix {
   // Add*() functions below.
   void Reset(RowIndex num_rows);
 
+  // Api to add columns one at the time.
+  void AddEntryToCurrentColumn(RowIndex row, Fractional coeff);
+  void CloseCurrentColumn();
+
   // Adds a dense column to the CompactSparseMatrix (only the non-zero will be
   // actually stored). This work in O(input.size()) and returns the index of the
   // added column.
diff --git a/ortools/sat/2d_rectangle_presolve_test.cc b/ortools/sat/2d_rectangle_presolve_test.cc
index 899dc9cd37..20700a826d 100644
--- a/ortools/sat/2d_rectangle_presolve_test.cc
+++ b/ortools/sat/2d_rectangle_presolve_test.cc
@@ -48,13 +48,14 @@ using ::testing::IsEmpty;
 std::vector<Rectangle> BuildFromAsciiArt(std::string_view input) {
   std::vector<Rectangle> rectangles;
   std::vector<std::string_view> lines = absl::StrSplit(input, '\n');
+  const int num_lines = lines.size();
   for (int i = 0; i < lines.size(); i++) {
     for (int j = 0; j < lines[i].size(); j++) {
       if (lines[i][j] != ' ') {
         rectangles.push_back({.x_min = j,
                               .x_max = j + 1,
-                              .y_min = 2 * lines.size() - 2 * i,
-                              .y_max = 2 * lines.size() - 2 * i + 2});
+                              .y_min = 2 * num_lines - 2 * i,
+                              .y_max = 2 * num_lines - 2 * i + 2});
       }
     }
   }
diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel
index 5e1b5ba3d7..e4d9c5c1f9 100644
--- a/ortools/sat/BUILD.bazel
+++ b/ortools/sat/BUILD.bazel
@@ -2146,6 +2146,7 @@ cc_library(
         "//ortools/lp_data:base",
         "//ortools/lp_data:lp_data_utils",
         "//ortools/lp_data:scattered_vector",
+        "//ortools/lp_data:sparse",
         "//ortools/lp_data:sparse_column",
         "//ortools/util:bitset",
         "//ortools/util:rev",
@@ -2156,6 +2157,7 @@ cc_library(
         "@com_google_absl//absl/log:check",
         "@com_google_absl//absl/numeric:int128",
         "@com_google_absl//absl/strings",
+        "@com_google_absl//absl/strings:str_format",
         "@com_google_absl//absl/types:span",
     ],
 )
@@ -3404,6 +3406,26 @@ cc_library(
     ],
 )
 
+cc_test(
+    name = "work_assignment_test",
+    srcs = ["work_assignment_test.cc"],
+    deps = [
+        ":cp_model",
+        ":cp_model_cc_proto",
+        ":cp_model_checker",
+        ":cp_model_loader",
+        ":cp_model_solver",
+        ":integer",
+        ":model",
+        ":sat_parameters_cc_proto",
+        ":synchronization",
+        ":work_assignment",
+        "//ortools/base:gmock_main",
+        "//ortools/base:parse_text_proto",
+        "@com_google_absl//absl/strings:string_view",
+    ],
+)
+
 cc_test(
     name = "inclusion_test",
     size = "small",
diff --git a/ortools/sat/cp_model_checker.cc b/ortools/sat/cp_model_checker.cc
index 926d58d04a..8645013cd0 100644
--- a/ortools/sat/cp_model_checker.cc
+++ b/ortools/sat/cp_model_checker.cc
@@ -854,9 +854,9 @@ std::string ValidateSolutionHint(const CpModelProto& model) {
   if (hint.vars().size() != hint.values().size()) {
     return "Invalid solution hint: vars and values do not have the same size.";
   }
-  for (const int ref : hint.vars()) {
-    if (!VariableReferenceIsValid(model, ref)) {
-      return absl::StrCat("Invalid variable reference in solution hint: ", ref);
+  for (const int var : hint.vars()) {
+    if (!VariableIndexIsValid(model, var)) {
+      return absl::StrCat("Invalid variable in solution hint: ", var);
     }
   }
 
diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc
index e9cceb9d4d..eaa97cdc80 100644
--- a/ortools/sat/cp_model_expand.cc
+++ b/ortools/sat/cp_model_expand.cc
@@ -1768,7 +1768,8 @@ void CompressAndExpandPositiveTable(ConstraintProto* ct,
     }
   }
 
-  VLOG(2) << "Table compression" << " var=" << vars.size()
+  VLOG(2) << "Table compression"
+          << " var=" << vars.size()
           << " cost=" << domain_sizes.size() - vars.size()
           << " tuples= " << num_tuples_before_compression << " -> "
           << num_tuples_after_first_compression << " -> "
diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc
index a3ccabe003..48404abdfd 100644
--- a/ortools/sat/cp_model_lns.cc
+++ b/ortools/sat/cp_model_lns.cc
@@ -448,7 +448,7 @@ std::vector<int> NeighborhoodGeneratorHelper::GetActiveIntervals(
                              initial_solution);
 }
 
-std::vector<std::pair<int, int>>
+std::vector<NeighborhoodGeneratorHelper::ActiveRectangle>
 NeighborhoodGeneratorHelper::GetActiveRectangles(
     const CpSolverResponse& initial_solution) const {
   const std::vector<int> active_intervals =
@@ -456,7 +456,7 @@ NeighborhoodGeneratorHelper::GetActiveRectangles(
   const absl::flat_hash_set<int> active_intervals_set(active_intervals.begin(),
                                                       active_intervals.end());
 
-  std::vector<std::pair<int, int>> active_rectangles;
+  absl::flat_hash_map<std::pair<int, int>, std::vector<int>> active_rectangles;
   for (const int ct_index : TypeToConstraints(ConstraintProto::kNoOverlap2D)) {
     const NoOverlap2DConstraintProto& ct =
         model_proto_.constraints(ct_index).no_overlap_2d();
@@ -465,12 +465,20 @@ NeighborhoodGeneratorHelper::GetActiveRectangles(
       const int y_i = ct.y_intervals(i);
       if (active_intervals_set.contains(x_i) ||
           active_intervals_set.contains(y_i)) {
-        active_rectangles.push_back({x_i, y_i});
+        active_rectangles[{x_i, y_i}].push_back(ct_index);
       }
     }
   }
 
-  return active_rectangles;
+  std::vector<ActiveRectangle> results;
+  for (const auto& [rectangle, no_overlap_2d_constraints] : active_rectangles) {
+    ActiveRectangle& result = results.emplace_back();
+    result.x_interval = rectangle.first;
+    result.y_interval = rectangle.second;
+    result.no_overlap_2d_constraints = {no_overlap_2d_constraints.begin(),
+                                        no_overlap_2d_constraints.end()};
+  }
+  return results;
 }
 
 std::vector<std::vector<int>>
@@ -1718,8 +1726,8 @@ Neighborhood DecompositionGraphNeighborhoodGenerator::Generate(
     VLOG(2) << "#relaxed " << relaxed_variables.size() << " #zero_score "
             << num_zero_score << " max_width " << max_width
             << " (size,min_width)_after_100 (" << size_at_min_width_after_100
-            << "," << min_width_after_100 << ") " << " final_width "
-            << pq.Size();
+            << "," << min_width_after_100 << ") "
+            << " final_width " << pq.Size();
   }
 
   return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
@@ -2256,14 +2264,125 @@ Neighborhood SchedulingResourceWindowsNeighborhoodGenerator::Generate(
 Neighborhood RandomRectanglesPackingNeighborhoodGenerator::Generate(
     const CpSolverResponse& initial_solution, SolveData& data,
     absl::BitGenRef random) {
-  std::vector<std::pair<int, int>> rectangles_to_freeze =
+  std::vector<ActiveRectangle> rectangles_to_freeze =
       helper_.GetActiveRectangles(initial_solution);
   GetRandomSubset(1.0 - data.difficulty, &rectangles_to_freeze, random);
 
   absl::flat_hash_set<int> variables_to_freeze;
-  for (const auto& [x, y] : rectangles_to_freeze) {
-    InsertVariablesFromConstraint(helper_.ModelProto(), x, variables_to_freeze);
-    InsertVariablesFromConstraint(helper_.ModelProto(), y, variables_to_freeze);
+  for (const ActiveRectangle& rectangle : rectangles_to_freeze) {
+    InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.x_interval,
+                                  variables_to_freeze);
+    InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.y_interval,
+                                  variables_to_freeze);
+  }
+
+  return helper_.FixGivenVariables(initial_solution, variables_to_freeze);
+}
+
+Neighborhood RectanglesPackingRelaxTwoNeighborhoodsGenerator::Generate(
+    const CpSolverResponse& initial_solution, SolveData& data,
+    absl::BitGenRef random) {
+  // First pick a pair of rectangles.
+  std::vector<ActiveRectangle> all_active_rectangles =
+      helper_.GetActiveRectangles(initial_solution);
+  if (all_active_rectangles.size() <= 2) return helper_.FullNeighborhood();
+
+  const int first_idx =
+      absl::Uniform<int>(random, 0, all_active_rectangles.size());
+  int second_idx =
+      absl::Uniform<int>(random, 0, all_active_rectangles.size() - 1);
+  if (second_idx >= first_idx) {
+    second_idx++;
+  }
+
+  const ActiveRectangle& chosen_rectangle_1 = all_active_rectangles[first_idx];
+  const ActiveRectangle& chosen_rectangle_2 = all_active_rectangles[second_idx];
+
+  const auto get_rectangle = [&initial_solution, helper = &helper_](
+                                 const ActiveRectangle& rectangle) {
+    const int x_interval_idx = rectangle.x_interval;
+    const int y_interval_idx = rectangle.y_interval;
+    const ConstraintProto& x_interval_ct =
+        helper->ModelProto().constraints(x_interval_idx);
+    const ConstraintProto& y_interval_ct =
+        helper->ModelProto().constraints(y_interval_idx);
+    return Rectangle{.x_min = GetLinearExpressionValue(
+                         x_interval_ct.interval().start(), initial_solution),
+                     .x_max = GetLinearExpressionValue(
+                         x_interval_ct.interval().end(), initial_solution),
+                     .y_min = GetLinearExpressionValue(
+                         y_interval_ct.interval().start(), initial_solution),
+                     .y_max = GetLinearExpressionValue(
+                         y_interval_ct.interval().end(), initial_solution)};
+  };
+
+  // TODO(user): This computes the distance between the center of the
+  // rectangles. We could use the real distance between the closest points, but
+  // not sure it is worth the extra complexity.
+  const auto compute_rectangle_distance = [](const Rectangle& rect1,
+                                             const Rectangle& rect2) {
+    return (static_cast<double>(rect1.x_min.value()) + rect1.x_max.value() -
+            rect2.x_min.value() - rect2.x_max.value()) *
+           (static_cast<double>(rect1.y_min.value()) + rect1.y_max.value() -
+            rect2.y_min.value() - rect2.y_max.value());
+  };
+  const Rectangle rect1 = get_rectangle(chosen_rectangle_1);
+  const Rectangle rect2 = get_rectangle(chosen_rectangle_2);
+
+  // Now compute a neighborhood around each rectangle. Note that we only
+  // consider two rectangles as potential neighbors if they are part of the same
+  // no_overlap_2d constraint.
+  absl::flat_hash_set<int> variables_to_freeze;
+  std::vector<std::pair<int, double>> distances1;
+  std::vector<std::pair<int, double>> distances2;
+  distances1.reserve(all_active_rectangles.size());
+  distances2.reserve(all_active_rectangles.size());
+  for (int i = 0; i < all_active_rectangles.size(); ++i) {
+    const ActiveRectangle& rectangle = all_active_rectangles[i];
+    InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.x_interval,
+                                  variables_to_freeze);
+    InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.y_interval,
+                                  variables_to_freeze);
+
+    const Rectangle rect = get_rectangle(rectangle);
+    const bool same_no_overlap_as_rect1 =
+        absl::c_any_of(chosen_rectangle_1.no_overlap_2d_constraints,
+                       [&rectangle](const int c) {
+                         return rectangle.no_overlap_2d_constraints.contains(c);
+                       });
+    const bool same_no_overlap_as_rect2 =
+        absl::c_any_of(chosen_rectangle_2.no_overlap_2d_constraints,
+                       [&rectangle](const int c) {
+                         return rectangle.no_overlap_2d_constraints.contains(c);
+                       });
+    if (same_no_overlap_as_rect1) {
+      distances1.push_back({i, compute_rectangle_distance(rect1, rect)});
+    }
+    if (same_no_overlap_as_rect2) {
+      distances2.push_back({i, compute_rectangle_distance(rect2, rect)});
+    }
+  }
+  const int num_to_sample_each =
+      data.difficulty * all_active_rectangles.size() / 2;
+  std::sort(distances1.begin(), distances1.end(),
+            [](const auto& a, const auto& b) { return a.second < b.second; });
+  std::sort(distances2.begin(), distances2.end(),
+            [](const auto& a, const auto& b) { return a.second < b.second; });
+  absl::flat_hash_set<int> variables_to_relax;
+  for (auto& samples : {distances1, distances2}) {
+    const int num_potential_samples = samples.size();
+    for (int i = 0; i < std::min(num_potential_samples, num_to_sample_each);
+         ++i) {
+      const int rectangle_idx = samples[i].first;
+      const ActiveRectangle& rectangle = all_active_rectangles[rectangle_idx];
+      InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.x_interval,
+                                    variables_to_relax);
+      InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.y_interval,
+                                    variables_to_relax);
+    }
+  }
+  for (const int v : variables_to_relax) {
+    variables_to_freeze.erase(v);
   }
 
   return helper_.FixGivenVariables(initial_solution, variables_to_freeze);
@@ -2272,13 +2391,13 @@ Neighborhood RandomRectanglesPackingNeighborhoodGenerator::Generate(
 Neighborhood RandomPrecedencesPackingNeighborhoodGenerator::Generate(
     const CpSolverResponse& initial_solution, SolveData& data,
     absl::BitGenRef random) {
-  std::vector<std::pair<int, int>> rectangles_to_relax =
+  std::vector<ActiveRectangle> rectangles_to_relax =
       helper_.GetActiveRectangles(initial_solution);
   GetRandomSubset(data.difficulty, &rectangles_to_relax, random);
   std::vector<int> intervals_to_relax;
-  for (const auto& [x, y] : rectangles_to_relax) {
-    intervals_to_relax.push_back(x);
-    intervals_to_relax.push_back(y);
+  for (const ActiveRectangle& rect : rectangles_to_relax) {
+    intervals_to_relax.push_back(rect.x_interval);
+    intervals_to_relax.push_back(rect.y_interval);
   }
   gtl::STLSortAndRemoveDuplicates(&intervals_to_relax);
 
@@ -2289,13 +2408,14 @@ Neighborhood RandomPrecedencesPackingNeighborhoodGenerator::Generate(
 Neighborhood SlicePackingNeighborhoodGenerator::Generate(
     const CpSolverResponse& initial_solution, SolveData& data,
     absl::BitGenRef random) {
-  const std::vector<std::pair<int, int>> active_rectangles =
+  const std::vector<ActiveRectangle> active_rectangles =
       helper_.GetActiveRectangles(initial_solution);
   const bool use_first_dimension = absl::Bernoulli(random, 0.5);
   std::vector<int> projected_intervals;
   projected_intervals.reserve(active_rectangles.size());
-  for (const auto& [x, y] : active_rectangles) {
-    projected_intervals.push_back(use_first_dimension ? x : y);
+  for (const ActiveRectangle& rect : active_rectangles) {
+    projected_intervals.push_back(use_first_dimension ? rect.x_interval
+                                                      : rect.y_interval);
   }
 
   const TimePartition partition = PartitionIndicesAroundRandomTimeWindow(
@@ -2310,10 +2430,10 @@ Neighborhood SlicePackingNeighborhoodGenerator::Generate(
   for (int index = 0; index < active_rectangles.size(); ++index) {
     if (indices_to_fix[index]) {
       InsertVariablesFromConstraint(helper_.ModelProto(),
-                                    active_rectangles[index].first,
+                                    active_rectangles[index].x_interval,
                                     variables_to_freeze);
       InsertVariablesFromConstraint(helper_.ModelProto(),
-                                    active_rectangles[index].second,
+                                    active_rectangles[index].y_interval,
                                     variables_to_freeze);
     }
   }
diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h
index dd92ff04d3..00856cbe1c 100644
--- a/ortools/sat/cp_model_lns.h
+++ b/ortools/sat/cp_model_lns.h
@@ -219,7 +219,14 @@ class NeighborhoodGeneratorHelper : public SubSolver {
   // lns_focus_on_performed_intervals. If true, this method returns the list of
   // performed rectangles in the solution. If false, it returns all rectangles
   // of the model.
-  std::vector<std::pair<int, int>> GetActiveRectangles(
+  struct ActiveRectangle {
+    int x_interval;
+    int y_interval;
+    // The set of no_overlap_2d constraints that both x_interval and y_interval
+    // are participating in.
+    absl::flat_hash_set<int> no_overlap_2d_constraints;
+  };
+  std::vector<ActiveRectangle> GetActiveRectangles(
       const CpSolverResponse& initial_solution) const;
 
   // Returns the set of unique intervals list appearing in a no_overlap,
@@ -356,6 +363,8 @@ class NeighborhoodGenerator {
       : name_(name), helper_(*helper), difficulty_(0.5) {}
   virtual ~NeighborhoodGenerator() = default;
 
+  using ActiveRectangle = NeighborhoodGeneratorHelper::ActiveRectangle;
+
   // Adds solve data about one "solved" neighborhood.
   struct SolveData {
     // The status of the sub-solve.
@@ -706,6 +715,21 @@ class RandomRectanglesPackingNeighborhoodGenerator
                         SolveData& data, absl::BitGenRef random) final;
 };
 
+// Only make sense for problems with no_overlap_2d constraints. This selects two
+// random rectangles and relax them alongside the closest rectangles to each one
+// of them. The idea is that this will find a better solution when there is a
+// cost function that would be improved by swapping the two rectangles.
+class RectanglesPackingRelaxTwoNeighborhoodsGenerator
+    : public NeighborhoodGenerator {
+ public:
+  explicit RectanglesPackingRelaxTwoNeighborhoodsGenerator(
+      NeighborhoodGeneratorHelper const* helper, absl::string_view name)
+      : NeighborhoodGenerator(name, helper) {}
+
+  Neighborhood Generate(const CpSolverResponse& initial_solution,
+                        SolveData& data, absl::BitGenRef random) final;
+};
+
 // Only make sense for problems with no_overlap_2d constraints. This select a
 // random set of rectangles (i.e. a pair of intervals) of the problem according
 // to the difficulty. Then add all implied precedences from the current
diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc
index 4b6720e106..1a1eeb358d 100644
--- a/ortools/sat/cp_model_presolve.cc
+++ b/ortools/sat/cp_model_presolve.cc
@@ -1704,6 +1704,60 @@ bool CpModelPresolver::PresolveIntProd(ConstraintProto* ct) {
     }
   }
 
+  if (ct->int_prod().exprs().size() == 2) {
+    const auto is_boolean_affine =
+        [context = context_](const LinearExpressionProto& expr) {
+          return expr.vars().size() == 1 && context->MinOf(expr.vars(0)) == 0 &&
+                 context->MaxOf(expr.vars(0)) == 1;
+        };
+    const LinearExpressionProto* boolean_linear = nullptr;
+    const LinearExpressionProto* other_linear = nullptr;
+    if (is_boolean_affine(ct->int_prod().exprs(0))) {
+      boolean_linear = &ct->int_prod().exprs(0);
+      other_linear = &ct->int_prod().exprs(1);
+    } else if (is_boolean_affine(ct->int_prod().exprs(1))) {
+      boolean_linear = &ct->int_prod().exprs(1);
+      other_linear = &ct->int_prod().exprs(0);
+    }
+    if (boolean_linear) {
+      // We have:
+      // (u + b * v) * other_expr = B, where `b` is a boolean variable.
+      //
+      // We can rewrite this as:
+      //   u * other_expr = B, if b = false;
+      //   (u + v) * other_expr = B, if b = true
+      ConstraintProto* constraint_for_false =
+          context_->working_model->add_constraints();
+      ConstraintProto* constraint_for_true =
+          context_->working_model->add_constraints();
+      constraint_for_true->add_enforcement_literal(boolean_linear->vars(0));
+      constraint_for_false->add_enforcement_literal(
+          NegatedRef(boolean_linear->vars(0)));
+      LinearConstraintProto* linear_for_false =
+          constraint_for_false->mutable_linear();
+      LinearConstraintProto* linear_for_true =
+          constraint_for_true->mutable_linear();
+
+      linear_for_false->add_domain(0);
+      linear_for_false->add_domain(0);
+      AddLinearExpressionToLinearConstraint(
+          *other_linear, boolean_linear->offset(), linear_for_false);
+      AddLinearExpressionToLinearConstraint(ct->int_prod().target(), -1,
+                                            linear_for_false);
+
+      linear_for_true->add_domain(0);
+      linear_for_true->add_domain(0);
+      AddLinearExpressionToLinearConstraint(
+          *other_linear, boolean_linear->offset() + boolean_linear->coeffs(0),
+          linear_for_true);
+      AddLinearExpressionToLinearConstraint(ct->int_prod().target(), -1,
+                                            linear_for_true);
+      context_->UpdateRuleStats("int_prod: boolean affine term");
+      context_->UpdateNewConstraintsVariableUsage();
+      return RemoveConstraint(ct);
+    }
+  }
+
   // For now, we only presolve the case where all variables are Booleans.
   const LinearExpressionProto target_expr = ct->int_prod().target();
   int target;
@@ -12472,6 +12526,27 @@ void CopyEverythingExceptVariablesAndConstraintsFieldsIntoContext(
   }
   if (in_model.has_solution_hint()) {
     *context->working_model->mutable_solution_hint() = in_model.solution_hint();
+
+    // We make sure the hint is within the variables domain.
+    //
+    // This allows to avoid overflow because we know evaluating constraints on
+    // the variables domains should be safe thanks to the initial validation.
+    const int num_terms = in_model.solution_hint().vars().size();
+    for (int i = 0; i < num_terms; ++i) {
+      const int var = in_model.solution_hint().vars(i);
+      const int64_t value = in_model.solution_hint().values(i);
+      const auto& domain = in_model.variables(var).domain();
+      if (domain.empty()) continue;  // UNSAT.
+      const int64_t min = domain[0];
+      const int64_t max = domain[domain.size() - 1];
+      if (value < min) {
+        context->UpdateRuleStats("hint: moved var hint within its domain.");
+        context->working_model->mutable_solution_hint()->set_values(i, min);
+      } else if (value > max) {
+        context->working_model->mutable_solution_hint()->set_values(i, max);
+        context->UpdateRuleStats("hint: moved var hint within its domain.");
+      }
+    }
   }
 }
 
diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc
index 728524ed36..dd60f67fa3 100644
--- a/ortools/sat/cp_model_solver.cc
+++ b/ortools/sat/cp_model_solver.cc
@@ -1663,6 +1663,12 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) {
                 helper, name_filter.LastName()),
             lns_params, helper, shared));
       }
+      if (name_filter.Keep("packing_swap_lns")) {
+        reentrant_interleaved_subsolvers.push_back(std::make_unique<LnsSolver>(
+            std::make_unique<RectanglesPackingRelaxTwoNeighborhoodsGenerator>(
+                helper, name_filter.LastName()),
+            lns_params, helper, shared));
+      }
       if (name_filter.Keep("packing_precedences_lns")) {
         reentrant_interleaved_subsolvers.push_back(std::make_unique<LnsSolver>(
             std::make_unique<RandomPrecedencesPackingNeighborhoodGenerator>(
diff --git a/ortools/sat/feasibility_jump.cc b/ortools/sat/feasibility_jump.cc
index 8219baf393..1083b3e6e0 100644
--- a/ortools/sat/feasibility_jump.cc
+++ b/ortools/sat/feasibility_jump.cc
@@ -86,7 +86,8 @@ bool JumpTable::JumpIsUpToDate(int var) const {
   if (abs(score - scores_[var]) / std::max(abs(score), 1.0) > 1e-2) {
     score_ok = false;
     LOG(ERROR) << "Incorrect score for var " << var << ": " << scores_[var]
-               << " (should be " << score << ") " << " delta = " << delta;
+               << " (should be " << score << ") "
+               << " delta = " << delta;
   }
   return delta == deltas_[var] && score_ok;
 }
diff --git a/ortools/sat/feasibility_pump.cc b/ortools/sat/feasibility_pump.cc
index 5ff2c7dd85..9a0cb4a761 100644
--- a/ortools/sat/feasibility_pump.cc
+++ b/ortools/sat/feasibility_pump.cc
@@ -360,14 +360,14 @@ void FeasibilityPump::L1DistanceMinimize() {
       const ColIndex norm_lhs_slack_variable =
           lp_data_.GetSlackVariable(norm_lhs_constraints_[col]);
       const double lhs_scaling_factor =
-          scaler_.VariableScalingFactor(norm_lhs_slack_variable);
+          scaler_.VariableScalingFactorWithSlack(norm_lhs_slack_variable);
       lp_data_.SetVariableBounds(
           norm_lhs_slack_variable, -glop::kInfinity,
           lhs_scaling_factor * integer_solution_[col.value()]);
       const ColIndex norm_rhs_slack_variable =
           lp_data_.GetSlackVariable(norm_rhs_constraints_[col]);
       const double rhs_scaling_factor =
-          scaler_.VariableScalingFactor(norm_rhs_slack_variable);
+          scaler_.VariableScalingFactorWithSlack(norm_rhs_slack_variable);
       lp_data_.SetVariableBounds(
           norm_rhs_slack_variable, -glop::kInfinity,
           -rhs_scaling_factor * integer_solution_[col.value()]);
diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc
index 1f5288284e..ec2d67e595 100644
--- a/ortools/sat/linear_constraint_manager.cc
+++ b/ortools/sat/linear_constraint_manager.cc
@@ -92,6 +92,15 @@ bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints(
   int new_size = 0;
   for (int i = 0; i < num_rows; ++i) {
     const ConstraintIndex constraint_index = lp_constraints_[i];
+    if (constraint_infos_[constraint_index].constraint.num_terms == 0) {
+      // Remove empty constraint.
+      //
+      // TODO(user): If the constraint is infeasible we could detect unsat
+      // right away, but hopefully this is a case where the propagation part
+      // of the solver can detect that too.
+      constraint_infos_[constraint_index].is_in_lp = false;
+      continue;
+    }
 
     // Constraints that are not tight in the current solution have a basic
     // status. We remove the ones that have been inactive in the last recent
@@ -245,7 +254,8 @@ bool LinearConstraintManager::AddCut(LinearConstraint ct, std::string type_name,
 
   // Only add cut with sufficient efficacy.
   if (violation / l2_norm < 1e-4) {
-    VLOG(3) << "BAD Cut '" << type_name << "'" << " size=" << ct.num_terms
+    VLOG(3) << "BAD Cut '" << type_name << "'"
+            << " size=" << ct.num_terms
             << " max_magnitude=" << ComputeInfinityNorm(ct)
             << " norm=" << l2_norm << " violation=" << violation
             << " eff=" << violation / l2_norm << " " << extra_info;
@@ -764,6 +774,8 @@ bool LinearConstraintManager::ChangeLp(glop::BasisState* solution_state,
 void LinearConstraintManager::AddAllConstraintsToLp() {
   for (ConstraintIndex i(0); i < constraint_infos_.size(); ++i) {
     if (constraint_infos_[i].is_in_lp) continue;
+    if (constraint_infos_[i].constraint.num_terms == 0) continue;
+
     constraint_infos_[i].is_in_lp = true;
     lp_constraints_.push_back(i);
   }
diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc
index 63237ee2f7..d308b1ae9d 100644
--- a/ortools/sat/linear_programming_constraint.cc
+++ b/ortools/sat/linear_programming_constraint.cc
@@ -31,6 +31,7 @@
 #include "absl/log/check.h"
 #include "absl/numeric/int128.h"
 #include "absl/strings/str_cat.h"
+#include "absl/strings/str_format.h"
 #include "absl/strings/string_view.h"
 #include "absl/types/span.h"
 #include "ortools/algorithms/binary_search.h"
@@ -45,6 +46,7 @@
 #include "ortools/lp_data/lp_data_utils.h"
 #include "ortools/lp_data/lp_types.h"
 #include "ortools/lp_data/scattered_vector.h"
+#include "ortools/lp_data/sparse.h"
 #include "ortools/lp_data/sparse_column.h"
 #include "ortools/sat/cp_model_mapping.h"
 #include "ortools/sat/cuts.h"
@@ -373,8 +375,6 @@ void LinearProgrammingConstraint::SetObjectiveCoefficient(IntegerVariable ivar,
 // for TSP for instance where the number of edges is large, but only a small
 // fraction will be used in the optimal solution.
 bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
-  simplex_.NotifyThatMatrixIsChangedForNextSolve();
-
   // Fill integer_lp_.
   integer_lp_.clear();
   integer_lp_cols_.clear();
@@ -384,6 +384,10 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
   const auto& all_constraints = constraint_manager_.AllConstraints();
   for (const auto index : constraint_manager_.LpConstraints()) {
     const LinearConstraint& ct = all_constraints[index].constraint;
+    if (ct.lb > ct.ub) {
+      VLOG(1) << "Trivial infeasible bound in an LP constraint";
+      return false;
+    }
 
     integer_lp_.push_back(LinearConstraintInternal());
     LinearConstraintInternal& new_ct = integer_lp_.back();
@@ -391,16 +395,15 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
     new_ct.ub = ct.ub;
     new_ct.lb_is_trivial = all_constraints[index].lb_is_trivial;
     new_ct.ub_is_trivial = all_constraints[index].ub_is_trivial;
-    const int size = ct.num_terms;
-    if (ct.lb > ct.ub) {
-      VLOG(1) << "Trivial infeasible bound in an LP constraint";
-      return false;
-    }
 
     IntegerValue infinity_norm = 0;
     infinity_norm = std::max(infinity_norm, IntTypeAbs(ct.lb));
     infinity_norm = std::max(infinity_norm, IntTypeAbs(ct.ub));
     new_ct.start_in_buffer = integer_lp_cols_.size();
+
+    // TODO(user): Make sure we don't have empty constraint!
+    // this currently can happen in some corner cases.
+    const int size = ct.num_terms;
     new_ct.num_terms = size;
     for (int i = 0; i < size; ++i) {
       // We only use positive variable inside this class.
@@ -418,12 +421,6 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
         integer_lp_cols_.data() + new_ct.start_in_buffer + new_ct.num_terms));
   }
 
-  // Copy the integer_lp_ into lp_data_.
-  lp_data_.Clear();
-  for (int i = 0; i < integer_variables_.size(); ++i) {
-    CHECK_EQ(glop::ColIndex(i), lp_data_.CreateNewVariable());
-  }
-
   // We remove fixed variables from the objective. This should help the LP
   // scaling, but also our integer reason computation.
   int new_size = 0;
@@ -438,49 +435,26 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
     objective_infinity_norm_ =
         std::max(objective_infinity_norm_, IntTypeAbs(entry.second));
     integer_objective_[new_size++] = entry;
-    lp_data_.SetObjectiveCoefficient(entry.first, ToDouble(entry.second));
   }
+  integer_objective_.resize(new_size);
   objective_infinity_norm_ =
       std::max(objective_infinity_norm_, IntTypeAbs(integer_objective_offset_));
-  integer_objective_.resize(new_size);
-  lp_data_.SetObjectiveOffset(ToDouble(integer_objective_offset_));
 
-  for (const LinearConstraintInternal& ct : integer_lp_) {
-    const ConstraintIndex row = lp_data_.CreateNewConstraint();
+  // Scale everything.
+  // TODO(user): As we have an idea of the LP optimal after the first solves,
+  // maybe we can adapt the scaling accordingly.
+  ComputeIntegerLpScalingFactors();
 
-    // TODO(user): Using trivial bound might be good for things like
-    // sum bool <= 1 since setting the slack in [0, 1] can lead to bound flip in
-    // the simplex. However if the bound is large, maybe it make more sense to
-    // use +/- infinity.
-    const double infinity = std::numeric_limits<double>::infinity();
-    lp_data_.SetConstraintBounds(
-        row, ct.lb_is_trivial ? -infinity : ToDouble(ct.lb),
-        ct.ub_is_trivial ? +infinity : ToDouble(ct.ub));
-    for (int i = 0; i < ct.num_terms; ++i) {
-      const int index = ct.start_in_buffer + i;
-      lp_data_.SetCoefficient(row, integer_lp_cols_[index],
-                              ToDouble(integer_lp_coeffs_[index]));
-    }
-  }
-  lp_data_.NotifyThatColumnsAreClean();
+  // Tricky: we use level zero bounds here for the second scaling step below.
+  FillLpData();
 
-  // We scale the LP using the level zero bounds that we later override
-  // with the current ones.
-  //
-  // TODO(user): As part of the scaling, we may also want to shift the initial
-  // variable bounds so that each variable contain the value zero in their
-  // domain. Maybe just once and for all at the beginning.
-  const int num_vars = integer_variables_.size();
-  for (int i = 0; i < num_vars; i++) {
-    const IntegerVariable cp_var = integer_variables_[i];
-    const double lb = ToDouble(integer_trail_->LevelZeroLowerBound(cp_var));
-    const double ub = ToDouble(integer_trail_->LevelZeroUpperBound(cp_var));
-    lp_data_.SetVariableBounds(glop::ColIndex(i), lb, ub);
-  }
+  // Fills the helper.
+  scaler_.ConfigureFromFactors(row_factors_, col_factors_);
+  scaler_.AverageCostScaling(&obj_with_slack_);
+  scaler_.ContainOneBoundScaling(simplex_.MutableLowerBounds(),
+                                 simplex_.MutableUpperBounds());
 
-  // TODO(user): As we have an idea of the LP optimal after the first solves,
-  // maybe we can adapt the scaling accordingly.
-  scaler_.Scale(simplex_params_, &lp_data_);
+  // Since we used level zero bounds above, fix them.
   UpdateBoundsOfLpVariables();
 
   // Set the information for the step to polish the LP basis. All our variables
@@ -488,6 +462,7 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
   // binary variables.
   if (parameters_.polish_lp_solution()) {
     simplex_.ClearIntegralityScales();
+    const int num_vars = integer_variables_.size();
     for (int i = 0; i < num_vars; ++i) {
       const IntegerVariable cp_var = integer_variables_[i];
       const IntegerValue lb = integer_trail_->LevelZeroLowerBound(cp_var);
@@ -499,13 +474,184 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
     }
   }
 
-  lp_data_.NotifyThatColumnsAreClean();
-  VLOG(3) << "LP relaxation: " << lp_data_.GetDimensionString() << ". "
+  VLOG(3) << "LP relaxation: " << integer_lp_.size() << " x "
+          << integer_variables_.size() << ". "
           << constraint_manager_.AllConstraints().size()
           << " Managed constraints.";
   return true;
 }
 
+// TODO(user): This is a duplicate of glop scaling code, but it allows to
+// work directly on our representation...
+void LinearProgrammingConstraint::ComputeIntegerLpScalingFactors() {
+  const int num_rows = integer_lp_.size();
+  const int num_cols = integer_variables_.size();
+
+  // Assign vectors.
+  const double infinity = std::numeric_limits<double>::infinity();
+  row_factors_.assign(num_rows, 1.0);
+  col_factors_.assign(num_cols, 1.0);
+
+  // Cache pointers to avoid refetching them.
+  IntegerValue* coeffs = integer_lp_coeffs_.data();
+  glop::ColIndex* cols = integer_lp_cols_.data();
+  double* row_factors = row_factors_.data();
+  double* col_factors = col_factors_.data();
+
+  col_min_.assign(num_cols, infinity);
+  col_max_.assign(num_cols, 0.0);
+  double* col_min = col_min_.data();
+  double* col_max = col_max_.data();
+
+  for (int i = 0; i < 4; ++i) {
+    // Scale row geometrically.
+    for (int row = 0; row < num_rows; ++row) {
+      double min_scaled = +infinity;
+      double max_scaled = 0.0;
+      const LinearConstraintInternal& ct = integer_lp_[RowIndex(row)];
+      for (int i = 0; i < ct.num_terms; ++i) {
+        const int index = ct.start_in_buffer + i;
+        const int col = cols[index].value();
+        const double coeff = static_cast<double>(coeffs[index].value());
+        const double scaled_magnitude = col_factors[col] * std::abs(coeff);
+        min_scaled = std::min(min_scaled, scaled_magnitude);
+        max_scaled = std::max(max_scaled, scaled_magnitude);
+      }
+
+      if (ct.num_terms == 0) continue;
+      const Fractional factor(std::sqrt(max_scaled * min_scaled));
+      row_factors[row] = 1.0 / factor;
+    }
+
+    // Scale columns geometrically.
+    for (int row = 0; row < num_rows; ++row) {
+      const double row_factor = row_factors[row];
+      const LinearConstraintInternal& ct = integer_lp_[RowIndex(row)];
+      for (int i = 0; i < ct.num_terms; ++i) {
+        const int index = ct.start_in_buffer + i;
+        const int col = cols[index].value();
+        const double coeff = static_cast<double>(coeffs[index].value());
+        const double scaled_magnitude = row_factor * std::abs(coeff);
+        col_min[col] = std::min(col_min[col], scaled_magnitude);
+        col_max[col] = std::max(col_max[col], scaled_magnitude);
+      }
+    }
+    for (int col = 0; col < num_cols; ++col) {
+      if (col_min[col] == infinity) continue;  // Empty.
+      col_factors[col] = 1.0 / std::sqrt(col_min[col] * col_max[col]);
+
+      // Reset, in case we have many fixed variable, faster than assign again.
+      col_min[col] = infinity;
+      col_max[col] = 0;
+    }
+  }
+
+  // Now we equilibrate (i.e. just divide by the max) the row
+  for (int row = 0; row < num_rows; ++row) {
+    double max_scaled = 0.0;
+    const LinearConstraintInternal& ct = integer_lp_[RowIndex(row)];
+    for (int i = 0; i < ct.num_terms; ++i) {
+      const int index = ct.start_in_buffer + i;
+      const int col = cols[index].value();
+      const double coeff = static_cast<double>(coeffs[index].value());
+      const double scaled_magnitude = col_factors[col] * std::abs(coeff);
+      max_scaled = std::max(max_scaled, scaled_magnitude);
+    }
+    if (ct.num_terms == 0) continue;
+    row_factors[row] = 1.0 / max_scaled;
+  }
+
+  // And finally the columns.
+  for (int row = 0; row < num_rows; ++row) {
+    const double row_factor = row_factors[row];
+    const LinearConstraintInternal& ct = integer_lp_[RowIndex(row)];
+    for (int i = 0; i < ct.num_terms; ++i) {
+      const int index = ct.start_in_buffer + i;
+      const int col = cols[index].value();
+      const double coeff = static_cast<double>(coeffs[index].value());
+      const double scaled_magnitude = row_factor * std::abs(coeff);
+      col_max[col] = std::max(col_max[col], scaled_magnitude);
+    }
+  }
+  for (int col = 0; col < num_cols; ++col) {
+    if (col_max[col] == 0) continue;  // Empty.
+    col_factors[col] = 1.0 / col_max[col];
+  }
+}
+
+void LinearProgrammingConstraint::FillLpData() {
+  const int num_rows = integer_lp_.size();
+  const int num_cols = integer_variables_.size();
+  IntegerValue* coeffs = integer_lp_coeffs_.data();
+  glop::ColIndex* cols = integer_lp_cols_.data();
+  double* row_factors = row_factors_.data();
+  double* col_factors = col_factors_.data();
+
+  // Now fill the tranposed matrix
+  glop::CompactSparseMatrix* data = simplex_.MutableTransposedMatrixWithSlack();
+  data->Reset(glop::RowIndex(num_cols + num_rows));
+  for (int row = 0; row < num_rows; ++row) {
+    const LinearConstraintInternal& ct = integer_lp_[RowIndex(row)];
+    const double row_factor = row_factors[row];
+    for (int i = 0; i < ct.num_terms; ++i) {
+      const int index = ct.start_in_buffer + i;
+      const int col = cols[index].value();
+      const double coeff = static_cast<double>(coeffs[index].value());
+      const double scaled_coeff = row_factor * col_factors[col] * coeff;
+      data->AddEntryToCurrentColumn(RowIndex(col), scaled_coeff);
+    }
+
+    // Add slack.
+    data->AddEntryToCurrentColumn(RowIndex(num_cols + row), 1.0);
+
+    // Close column.
+    data->CloseCurrentColumn();
+  }
+
+  // Fill and scale the objective.
+  const glop::ColIndex num_cols_with_slacks(num_rows + num_cols);
+  obj_with_slack_.assign(num_cols_with_slacks, 0.0);
+  for (const auto [col, value] : integer_objective_) {
+    obj_with_slack_[col] = ToDouble(value) * col_factors[col.value()];
+  }
+
+  // Fill and scales the bound.
+  simplex_.MutableLowerBounds()->resize(num_cols_with_slacks);
+  simplex_.MutableUpperBounds()->resize(num_cols_with_slacks);
+  Fractional* lb_with_slack = simplex_.MutableLowerBounds()->data();
+  Fractional* ub_with_slack = simplex_.MutableUpperBounds()->data();
+  const double infinity = std::numeric_limits<double>::infinity();
+  for (int row = 0; row < integer_lp_.size(); ++row) {
+    const LinearConstraintInternal& ct = integer_lp_[glop::RowIndex(row)];
+
+    // TODO(user): Using trivial bound might be good for things like
+    // sum bool <= 1 since setting the slack in [0, 1] can lead to bound flip in
+    // the simplex. However if the bound is large, maybe it make more sense to
+    // use +/- infinity.
+    const double factor = row_factors[row];
+    lb_with_slack[num_cols + row] =
+        ct.ub_is_trivial ? -infinity : ToDouble(-ct.ub) * factor;
+    ub_with_slack[num_cols + row] =
+        ct.lb_is_trivial ? +infinity : ToDouble(-ct.lb) * factor;
+  }
+
+  // We scale the LP using the level zero bounds that we later override
+  // with the current ones.
+  //
+  // TODO(user): As part of the scaling, we may also want to shift the initial
+  // variable bounds so that each variable contain the value zero in their
+  // domain. Maybe just once and for all at the beginning.
+  const int num_vars = integer_variables_.size();
+  for (int i = 0; i < num_vars; i++) {
+    const IntegerVariable cp_var = integer_variables_[i];
+    const double factor = col_factors[i];
+    lb_with_slack[i] =
+        ToDouble(integer_trail_->LevelZeroLowerBound(cp_var)) * factor;
+    ub_with_slack[i] =
+        ToDouble(integer_trail_->LevelZeroUpperBound(cp_var)) * factor;
+  }
+}
+
 void LinearProgrammingConstraint::FillReducedCostReasonIn(
     const glop::DenseRow& reduced_costs,
     std::vector<IntegerLiteral>* integer_reason) {
@@ -682,12 +828,17 @@ double LinearProgrammingConstraint::GetSolutionReducedCost(
 
 void LinearProgrammingConstraint::UpdateBoundsOfLpVariables() {
   const int num_vars = integer_variables_.size();
+  Fractional* lb_with_slack = simplex_.MutableLowerBounds()->data();
+  Fractional* ub_with_slack = simplex_.MutableUpperBounds()->data();
   for (int i = 0; i < num_vars; i++) {
     const IntegerVariable cp_var = integer_variables_[i];
-    const double lb = ToDouble(integer_trail_->LowerBound(cp_var));
-    const double ub = ToDouble(integer_trail_->UpperBound(cp_var));
+    const double lb =
+        static_cast<double>(integer_trail_->LowerBound(cp_var).value());
+    const double ub =
+        static_cast<double>(integer_trail_->UpperBound(cp_var).value());
     const double factor = scaler_.VariableScalingFactor(glop::ColIndex(i));
-    lp_data_.SetVariableBounds(glop::ColIndex(i), lb * factor, ub * factor);
+    lb_with_slack[i] = lb * factor;
+    ub_with_slack[i] = ub * factor;
   }
 }
 
@@ -697,7 +848,12 @@ bool LinearProgrammingConstraint::SolveLp() {
     lp_at_level_zero_is_final_ = false;
   }
 
-  const auto status = simplex_.Solve(lp_data_, time_limit_);
+  const double unscaling_factor = 1.0 / scaler_.ObjectiveScalingFactor();
+  const double offset_before_unscaling =
+      ToDouble(integer_objective_offset_) * scaler_.ObjectiveScalingFactor();
+  const auto status = simplex_.MinimizeFromTransposedMatrixWithSlack(
+      obj_with_slack_, unscaling_factor, offset_before_unscaling, time_limit_);
+
   state_ = simplex_.GetState();
   total_num_simplex_iterations_ += simplex_.GetNumberOfIterations();
   if (!status.ok()) {
@@ -711,19 +867,14 @@ bool LinearProgrammingConstraint::SolveLp() {
             << average_degeneracy_.CurrentAverage();
   }
 
-  // By default we assume the matrix is unchanged.
-  // This will be reset by CreateLpFromConstraintManager().
-  simplex_.NotifyThatMatrixIsUnchangedForNextSolve();
-
   const int status_as_int = static_cast<int>(simplex_.GetProblemStatus());
   if (status_as_int >= num_solves_by_status_.size()) {
     num_solves_by_status_.resize(status_as_int + 1);
   }
   num_solves_++;
   num_solves_by_status_[status_as_int]++;
-  VLOG(2) << lp_data_.GetDimensionString()
-          << " lvl:" << trail_->CurrentDecisionLevel() << " "
-          << simplex_.GetProblemStatus()
+  VLOG(2) << DimensionString() << " lvl:" << trail_->CurrentDecisionLevel()
+          << " " << simplex_.GetProblemStatus()
           << " iter:" << simplex_.GetNumberOfIterations()
           << " obj:" << simplex_.GetObjectiveValue() << " scaled:"
           << objective_definition_->ScaleObjective(
@@ -1285,11 +1436,17 @@ void LinearProgrammingConstraint::AddCGCuts() {
   const bool old_gomory = true;
 
   // Note that the index is permuted and do not correspond to a row.
-  const RowIndex num_rows = lp_data_.num_constraints();
+  const RowIndex num_rows(integer_lp_.size());
   for (RowIndex index(0); index < num_rows; ++index) {
     if (time_limit_->LimitReached()) break;
 
     const ColIndex basis_col = simplex_.GetBasis(index);
+
+    // If this variable is a slack, we ignore it. This is because the
+    // corresponding row is not tight under the given lp values.
+    if (old_gomory && basis_col >= integer_variables_.size()) continue;
+
+    // TODO(user): If the variable is a slack, the unscaling is wrong!
     const Fractional lp_value = GetVariableValueAtCpScale(basis_col);
 
     // Only consider fractional basis element. We ignore element that are close
@@ -1300,10 +1457,6 @@ void LinearProgrammingConstraint::AddCGCuts() {
     // also be just under it.
     if (std::abs(lp_value - std::round(lp_value)) < 0.01) continue;
 
-    // If this variable is a slack, we ignore it. This is because the
-    // corresponding row is not tight under the given lp values.
-    if (old_gomory && basis_col >= integer_variables_.size()) continue;
-
     // TODO(user): Avoid code duplication between the sparse/dense path.
     tmp_lp_multipliers_.clear();
     const glop::ScatteredRow& lambda = simplex_.GetUnitRowLeftInverse(index);
@@ -1464,9 +1617,12 @@ void LinearProgrammingConstraint::AddMirCuts() {
 
   // We compute all the rows that are tight, these will be used as the base row
   // for the MIR_n procedure below.
-  const int num_rows = lp_data_.num_constraints().value();
+  const int num_cols = integer_variables_.size();
+  const int num_rows = integer_lp_.size();
   std::vector<std::pair<RowIndex, IntegerValue>> base_rows;
   util_intops::StrongVector<RowIndex, double> row_weights(num_rows, 0.0);
+  Fractional* lb_with_slack = simplex_.MutableLowerBounds()->data();
+  Fractional* ub_with_slack = simplex_.MutableUpperBounds()->data();
   for (RowIndex row(0); row < num_rows; ++row) {
     // We only consider tight rows.
     // We use both the status and activity to have as much options as possible.
@@ -1476,12 +1632,14 @@ void LinearProgrammingConstraint::AddMirCuts() {
     // cannot be good.
     const auto status = simplex_.GetConstraintStatus(row);
     const double activity = simplex_.GetConstraintActivity(row);
-    if (activity > lp_data_.constraint_upper_bounds()[row] - 1e-4 ||
+    const double ct_lb = -ub_with_slack[num_cols + row.value()];
+    const double ct_ub = -lb_with_slack[num_cols + row.value()];
+    if (activity > ct_ub - 1e-4 ||
         status == glop::ConstraintStatus::AT_UPPER_BOUND ||
         status == glop::ConstraintStatus::FIXED_VALUE) {
       base_rows.push_back({row, IntegerValue(1)});
     }
-    if (activity < lp_data_.constraint_lower_bounds()[row] + 1e-4 ||
+    if (activity < ct_lb + 1e-4 ||
         status == glop::ConstraintStatus::AT_LOWER_BOUND ||
         status == glop::ConstraintStatus::FIXED_VALUE) {
       base_rows.push_back({row, IntegerValue(-1)});
@@ -1514,6 +1672,7 @@ void LinearProgrammingConstraint::AddMirCuts() {
   std::vector<double> weights;
   util_intops::StrongVector<RowIndex, bool> used_rows;
   std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers;
+  const auto matrix = simplex_.MatrixWithSlack().view();
   for (const std::pair<RowIndex, IntegerValue>& entry : base_rows) {
     if (time_limit_->LimitReached()) break;
     if (dtime_num_entries > 1e7) break;
@@ -1567,8 +1726,7 @@ void LinearProgrammingConstraint::AddMirCuts() {
         if (dense_cut[col] == 0) continue;
 
         max_magnitude = std::max(max_magnitude, IntTypeAbs(dense_cut[col]));
-        const int col_degree =
-            lp_data_.GetSparseColumn(col).num_entries().value();
+        const int col_degree = matrix.ColumnNumEntries(col).value();
         if (col_degree <= 1) continue;
         if (simplex_.GetVariableStatus(col) != glop::VariableStatus::BASIC) {
           continue;
@@ -1592,8 +1750,9 @@ void LinearProgrammingConstraint::AddMirCuts() {
       // What rows can we add to eliminate var_to_eliminate?
       std::vector<RowIndex> possible_rows;
       weights.clear();
-      for (const auto entry : lp_data_.GetSparseColumn(var_to_eliminate)) {
-        const RowIndex row = entry.row();
+      for (const auto entry_index : matrix.Column(var_to_eliminate)) {
+        const RowIndex row = matrix.EntryRow(entry_index);
+        const glop::Fractional coeff = matrix.EntryCoefficient(entry_index);
 
         // We disallow all the rows that contain a variable that we already
         // eliminated (or are about to). This mean that we choose rows that
@@ -1608,14 +1767,14 @@ void LinearProgrammingConstraint::AddMirCuts() {
         // still be chosen after the tight-one in most situation.
         bool add_row = false;
         if (!integer_lp_[row].ub_is_trivial) {
-          if (entry.coefficient() > 0.0) {
+          if (coeff > 0.0) {
             if (dense_cut[var_to_eliminate] < 0) add_row = true;
           } else {
             if (dense_cut[var_to_eliminate] > 0) add_row = true;
           }
         }
         if (!integer_lp_[row].lb_is_trivial) {
-          if (entry.coefficient() > 0.0) {
+          if (coeff > 0.0) {
             if (dense_cut[var_to_eliminate] > 0) add_row = true;
           } else {
             if (dense_cut[var_to_eliminate] < 0) add_row = true;
@@ -2409,11 +2568,8 @@ void LinearProgrammingConstraint::ReducedCostStrengtheningDeductions(
     double cp_objective_delta) {
   deductions_.clear();
 
-  // TRICKY: while simplex_.GetObjectiveValue() use the objective scaling factor
-  // stored in the lp_data_, all the other functions like GetReducedCost() or
-  // GetVariableValue() do not.
   const double lp_objective_delta =
-      cp_objective_delta / lp_data_.objective_scaling_factor();
+      cp_objective_delta / scaler_.ObjectiveScalingFactor();
   const int num_vars = integer_variables_.size();
   for (int i = 0; i < num_vars; i++) {
     const IntegerVariable cp_var = integer_variables_[i];
@@ -2596,5 +2752,10 @@ absl::Span<const IntegerValue> LinearProgrammingConstraint::IntegerLpRowCoeffs(
   return {integer_lp_coeffs_.data() + start, num_terms};
 }
 
+std::string LinearProgrammingConstraint::DimensionString() const {
+  return absl::StrFormat("%d rows, %d columns, %d entries", integer_lp_.size(),
+                         integer_variables_.size(), integer_lp_coeffs_.size());
+}
+
 }  // namespace sat
 }  // namespace operations_research
diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h
index ae98877708..bdf34fc686 100644
--- a/ortools/sat/linear_programming_constraint.h
+++ b/ortools/sat/linear_programming_constraint.h
@@ -190,7 +190,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   const std::vector<IntegerVariable>& integer_variables() const {
     return integer_variables_;
   }
-  std::string DimensionString() const { return lp_data_.GetDimensionString(); }
+  std::string DimensionString() const;
 
   // Returns a IntegerLiteral guided by the underlying LP constraints.
   //
@@ -408,6 +408,15 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   absl::Span<const glop::ColIndex> IntegerLpRowCols(glop::RowIndex row) const;
   absl::Span<const IntegerValue> IntegerLpRowCoeffs(glop::RowIndex row) const;
 
+  void ComputeIntegerLpScalingFactors();
+  void FillLpData();
+
+  // For ComputeIntegerLpScalingFactors().
+  std::vector<double> row_factors_;
+  std::vector<double> col_factors_;
+  std::vector<double> col_max_;
+  std::vector<double> col_min_;
+
   // This epsilon is related to the precision of the value/reduced_cost returned
   // by the LP once they have been scaled back into the CP domain. So for large
   // domain or cost coefficient, we may have some issues.
@@ -456,8 +465,9 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   // Underlying LP solver API.
   glop::GlopParameters simplex_params_;
   glop::BasisState state_;
-  glop::LinearProgram lp_data_;
+  glop::DenseRow obj_with_slack_;
   glop::RevisedSimplex simplex_;
+
   int64_t next_simplex_iter_ = 500;
 
   // For the scaling.
diff --git a/ortools/sat/lp_utils.cc b/ortools/sat/lp_utils.cc
index 847d64cf08..a95b8b70b1 100644
--- a/ortools/sat/lp_utils.cc
+++ b/ortools/sat/lp_utils.cc
@@ -1507,8 +1507,8 @@ bool ConvertBinaryMPModelProtoToBooleanProblem(const MPModelProto& mp_model,
     // Abort if the variable is not binary.
     if (!is_binary) {
       LOG(WARNING) << "The variable #" << var_id << " with name "
-                   << mp_var.name() << " is not binary. " << "lb: " << lb
-                   << " ub: " << ub;
+                   << mp_var.name() << " is not binary. "
+                   << "lb: " << lb << " ub: " << ub;
       return false;
     }
   }
diff --git a/ortools/sat/precedences.cc b/ortools/sat/precedences.cc
index 34b1559a0d..18b2a363d0 100644
--- a/ortools/sat/precedences.cc
+++ b/ortools/sat/precedences.cc
@@ -647,8 +647,9 @@ void PrecedencesPropagator::AddArc(
     // A self-arc is either plain SAT or plain UNSAT or it forces something on
     // the given offset_var or presence_literal_index. In any case it could be
     // presolved in something more efficient.
-    VLOG(1) << "Self arc! This could be presolved. " << "var:" << tail
-            << " offset:" << offset << " offset_var:" << offset_var
+    VLOG(1) << "Self arc! This could be presolved. "
+            << "var:" << tail << " offset:" << offset
+            << " offset_var:" << offset_var
             << " conditioned_by:" << presence_literals;
   }
 
diff --git a/ortools/sat/probing.cc b/ortools/sat/probing.cc
index 941c78c830..37d9b7d9d3 100644
--- a/ortools/sat/probing.cc
+++ b/ortools/sat/probing.cc
@@ -890,8 +890,9 @@ bool FailedLiteralProbingRound(ProbingOptions options, Model* model) {
   const bool limit_reached = time_limit->LimitReached() ||
                              time_limit->GetElapsedDeterministicTime() > limit;
   LOG_IF(INFO, options.log_info)
-      << "Probing. " << " num_probed: " << num_probed << " num_fixed: +"
-      << num_newly_fixed << " (" << num_fixed << "/" << num_variables << ")"
+      << "Probing. "
+      << " num_probed: " << num_probed << " num_fixed: +" << num_newly_fixed
+      << " (" << num_fixed << "/" << num_variables << ")"
       << " explicit_fix:" << num_explicit_fix
       << " num_conflicts:" << num_conflicts
       << " new_binary_clauses: " << num_new_binary
diff --git a/ortools/sat/sat_base.h b/ortools/sat/sat_base.h
index c342af7105..ced510384b 100644
--- a/ortools/sat/sat_base.h
+++ b/ortools/sat/sat_base.h
@@ -625,8 +625,9 @@ inline bool SatPropagator::PropagatePreconditionsAreSatisfied(
   if (propagation_trail_index_ < trail.Index() &&
       trail.Info(trail[propagation_trail_index_].Variable()).level !=
           trail.CurrentDecisionLevel()) {
-    LOG(INFO) << "Issue in '" << name_ << "':" << " propagation_trail_index_="
-              << propagation_trail_index_ << " trail_.Index()=" << trail.Index()
+    LOG(INFO) << "Issue in '" << name_ << "':"
+              << " propagation_trail_index_=" << propagation_trail_index_
+              << " trail_.Index()=" << trail.Index()
               << " level_at_propagation_index="
               << trail.Info(trail[propagation_trail_index_].Variable()).level
               << " current_decision_level=" << trail.CurrentDecisionLevel();
diff --git a/ortools/sat/sat_inprocessing.cc b/ortools/sat/sat_inprocessing.cc
index 1d6dd85cd2..24a7de25b8 100644
--- a/ortools/sat/sat_inprocessing.cc
+++ b/ortools/sat/sat_inprocessing.cc
@@ -704,8 +704,8 @@ bool StampingSimplifier::ComputeStampsForNextRound(bool log_info) {
 
   // TODO(user): compute some dtime, it is always zero currently.
   time_limit_->AdvanceDeterministicTime(dtime_);
-  LOG_IF(INFO, log_info) << "Prestamping." << " num_fixed: " << num_fixed_
-                         << " dtime: " << dtime_
+  LOG_IF(INFO, log_info) << "Prestamping."
+                         << " num_fixed: " << num_fixed_ << " dtime: " << dtime_
                          << " wtime: " << wall_timer.Get();
   return true;
 }
@@ -1259,7 +1259,8 @@ bool BoundedVariableElimination::DoOneRound(bool log_info) {
   dtime_ += 1e-8 * num_inspected_literals_;
   time_limit_->AdvanceDeterministicTime(dtime_);
   log_info |= VLOG_IS_ON(1);
-  LOG_IF(INFO, log_info) << "BVE." << " num_fixed: "
+  LOG_IF(INFO, log_info) << "BVE."
+                         << " num_fixed: "
                          << trail_->Index() - saved_trail_index
                          << " num_simplified_literals: " << num_simplifications_
                          << " num_blocked_clauses_: " << num_blocked_clauses_
diff --git a/ortools/sat/var_domination.cc b/ortools/sat/var_domination.cc
index f4368db2d5..ae589ac3e7 100644
--- a/ortools/sat/var_domination.cc
+++ b/ortools/sat/var_domination.cc
@@ -1304,8 +1304,8 @@ void ScanModelForDominanceDetection(PresolveContext& context,
     }
   }
   if (num_unconstrained_refs == 0 && num_dominated_refs == 0) return;
-  VLOG(1) << "Dominance:" << " num_unconstrained_refs="
-          << num_unconstrained_refs
+  VLOG(1) << "Dominance:"
+          << " num_unconstrained_refs=" << num_unconstrained_refs
           << " num_dominated_refs=" << num_dominated_refs
           << " num_dominance_relations=" << num_dominance_relations;
 }
diff --git a/ortools/sat/work_assignment_test.cc b/ortools/sat/work_assignment_test.cc
new file mode 100644
index 0000000000..ba4a6d9533
--- /dev/null
+++ b/ortools/sat/work_assignment_test.cc
@@ -0,0 +1,545 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/sat/work_assignment.h"
+
+#include <vector>
+
+#include "absl/strings/string_view.h"
+#include "gtest/gtest.h"
+#include "ortools/base/gmock.h"
+#include "ortools/base/parse_text_proto.h"
+#include "ortools/sat/cp_model.h"
+#include "ortools/sat/cp_model.pb.h"
+#include "ortools/sat/cp_model_checker.h"
+#include "ortools/sat/cp_model_loader.h"
+#include "ortools/sat/cp_model_solver.h"
+#include "ortools/sat/integer.h"
+#include "ortools/sat/model.h"
+#include "ortools/sat/sat_parameters.pb.h"
+#include "ortools/sat/synchronization.h"
+
+namespace operations_research {
+namespace sat {
+namespace {
+
+TEST(ProtoTrailTest, PushLevel) {
+  ProtoTrail p;
+  p.PushLevel({0, 0}, 0, 1);
+
+  EXPECT_EQ(p.MaxLevel(), 1);
+  EXPECT_EQ(p.Decision(1), ProtoLiteral(0, 0));
+  EXPECT_EQ(p.ObjectiveLb(1), 0);
+}
+
+TEST(ProtoTrailTest, AddImplications) {
+  ProtoTrail p;
+  p.PushLevel({0, 0}, 0, 1);
+  p.PushLevel({1, 0}, 1, 2);
+  p.PushLevel({2, 0}, 2, 3);
+  p.PushLevel({3, 0}, 2, 4);
+
+  p.AddImplication(2, {5, 0});
+  p.AddImplication(3, {6, 0});
+
+  EXPECT_THAT(p.Implications(2), testing::ElementsAre(ProtoLiteral(5, 0)));
+  EXPECT_THAT(p.Implications(3), testing::ElementsAre(ProtoLiteral(6, 0)));
+  p.SetLevelImplied(3);
+  EXPECT_THAT(p.Implications(2),
+              testing::UnorderedElementsAre(
+                  ProtoLiteral(5, 0), ProtoLiteral(2, 0), ProtoLiteral(6, 0)));
+}
+
+TEST(ProtoTrailTest, SetLevel1Implied) {
+  ProtoTrail p;
+  p.PushLevel({0, 0}, 0, 1);
+  p.PushLevel({1, 0}, 1, 2);
+  p.PushLevel({2, 0}, 2, 3);
+
+  p.SetLevelImplied(1);
+
+  EXPECT_THAT(p.NodeIds(0), testing::ElementsAre(1));
+  EXPECT_THAT(p.NodeIds(1), testing::ElementsAre(2));
+  EXPECT_THAT(p.NodeIds(2), testing::ElementsAre(3));
+  EXPECT_EQ(p.MaxLevel(), 2);
+  EXPECT_EQ(p.Decision(1), ProtoLiteral(1, 0));
+  EXPECT_EQ(p.Decision(2), ProtoLiteral(2, 0));
+  EXPECT_EQ(p.ObjectiveLb(1), 1);
+  EXPECT_EQ(p.ObjectiveLb(2), 2);
+}
+
+TEST(ProtoTrailTest, SetMidLevelImplied) {
+  ProtoTrail p;
+  p.PushLevel({0, 0}, 0, 1);
+  p.PushLevel({1, 0}, 1, 2);
+  p.PushLevel({2, 0}, 2, 3);
+
+  p.SetLevelImplied(2);
+
+  EXPECT_THAT(p.NodeIds(0), testing::IsEmpty());
+  EXPECT_THAT(p.NodeIds(1), testing::ElementsAre(1, 2));
+  EXPECT_THAT(p.NodeIds(2), testing::ElementsAre(3));
+  EXPECT_EQ(p.MaxLevel(), 2);
+  EXPECT_EQ(p.Decision(1), ProtoLiteral(0, 0));
+  EXPECT_EQ(p.Decision(2), ProtoLiteral(2, 0));
+  EXPECT_EQ(p.ObjectiveLb(1), 1);
+  EXPECT_EQ(p.ObjectiveLb(2), 2);
+}
+
+TEST(ProtoTrailTest, SetFinalLevelImplied) {
+  ProtoTrail p;
+  p.PushLevel({0, 0}, 0, 1);
+  p.PushLevel({1, 0}, 1, 2);
+  p.PushLevel({2, 0}, 2, 3);
+
+  p.SetLevelImplied(3);
+
+  EXPECT_THAT(p.NodeIds(0), testing::IsEmpty());
+  EXPECT_THAT(p.NodeIds(1), testing::ElementsAre(1));
+  EXPECT_THAT(p.NodeIds(2), testing::ElementsAre(2, 3));
+  EXPECT_EQ(p.MaxLevel(), 2);
+  EXPECT_EQ(p.Decision(1), ProtoLiteral(0, 0));
+  EXPECT_EQ(p.Decision(2), ProtoLiteral(1, 0));
+  EXPECT_EQ(p.ObjectiveLb(1), 0);
+  EXPECT_EQ(p.ObjectiveLb(2), 2);
+}
+
+TEST(ProtoTrailTest, SetMultiLevelImplied) {
+  ProtoTrail p;
+  p.PushLevel({0, 0}, 0, 1);
+  p.PushLevel({1, 0}, 1, 2);
+  p.PushLevel({2, 0}, 2, 3);
+
+  p.SetLevelImplied(3);
+  p.SetLevelImplied(1);
+
+  EXPECT_EQ(p.MaxLevel(), 1);
+  EXPECT_EQ(p.Decision(1), ProtoLiteral(1, 0));
+  EXPECT_EQ(p.ObjectiveLb(1), 2);
+}
+
+TEST(ProtoTrailTest, Clear) {
+  ProtoTrail p;
+  p.PushLevel({0, 0}, 0, 1);
+  p.PushLevel({1, 0}, 1, 2);
+  p.PushLevel({2, 0}, 2, 3);
+
+  p.Clear();
+
+  EXPECT_EQ(p.MaxLevel(), 0);
+}
+
+class SharedTreeSolveTest : public testing::TestWithParam<absl::string_view> {
+ public:
+  SatParameters GetParams() {
+    SatParameters params;
+    params.set_num_workers(4);
+    params.set_shared_tree_num_workers(4);
+    params.set_cp_model_presolve(false);
+    params.MergeFrom(
+        google::protobuf::contrib::parse_proto::ParseTextProtoOrDie(
+            GetParam()));
+    return params;
+  }
+};
+INSTANTIATE_TEST_SUITE_P(
+    SharedTreeParams, SharedTreeSolveTest,
+    testing::Values("shared_tree_worker_enable_trail_sharing:false",
+                    "shared_tree_worker_enable_trail_sharing:true"));
+
+TEST_P(SharedTreeSolveTest, SmokeTest) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var + 5 * bool_var);
+  Model model;
+  SatParameters params = GetParams();
+  model.Add(NewSatParameters(params));
+
+  CpSolverResponse response = SolveCpModel(model_builder.Build(), &model);
+
+  EXPECT_EQ(model.GetOrCreate<SharedTreeManager>()->NumWorkers(),
+            params.shared_tree_num_workers());
+  ASSERT_EQ(response.status(), OPTIMAL)
+      << "Validation: " << ValidateCpModel(model_builder.Build());
+  EXPECT_EQ(response.objective_value(), 5 + 3);
+  EXPECT_EQ(SolutionBooleanValue(response, bool_var), true);
+  EXPECT_EQ(SolutionIntegerValue(response, int_var), 3);
+}
+
+TEST_P(SharedTreeSolveTest, FeasiblePidgeonHoleSmokeTest) {
+  CpModelBuilder model_builder;
+  const int pidgeons = 10;
+  const int holes = 10;
+  std::vector<LinearExpr> count_per_hole(holes);
+  IntVar max_pidgeon_hole_product =
+      model_builder.NewIntVar({0, pidgeons * holes});
+  for (int i = 0; i < pidgeons; ++i) {
+    LinearExpr count_per_pidgeon;
+    for (int j = 0; j < holes; ++j) {
+      auto var = model_builder.NewBoolVar();
+      count_per_hole[j] += LinearExpr(var);
+      count_per_pidgeon += LinearExpr(var);
+      model_builder
+          .AddGreaterOrEqual(max_pidgeon_hole_product, (i + 1) * (j + 1))
+          .OnlyEnforceIf(var);
+    }
+    model_builder.AddEquality(count_per_pidgeon, 1);
+  }
+  for (const auto& count : count_per_hole) {
+    model_builder.AddLessOrEqual(count, 1);
+  }
+  Model model;
+  SatParameters params = GetParams();
+  model.Add(NewSatParameters(params));
+
+  CpSolverResponse response = SolveCpModel(model_builder.Build(), &model);
+
+  EXPECT_EQ(model.GetOrCreate<SharedTreeManager>()->NumWorkers(), 4);
+  EXPECT_EQ(response.status(), OPTIMAL);
+}
+
+TEST_P(SharedTreeSolveTest, InfeasiblePidgeonHoleSmokeTest) {
+  CpModelBuilder model_builder;
+  const int pidgeons = 10;
+  const int holes = 9;
+  std::vector<LinearExpr> count_per_hole(holes);
+  IntVar max_pidgeon_hole_product =
+      model_builder.NewIntVar({0, pidgeons * holes});
+  for (int i = 0; i < pidgeons; ++i) {
+    LinearExpr count_per_pidgeon;
+    for (int j = 0; j < holes; ++j) {
+      auto var = model_builder.NewBoolVar();
+      count_per_hole[j] += LinearExpr(var);
+      count_per_pidgeon += LinearExpr(var);
+      model_builder
+          .AddGreaterOrEqual(max_pidgeon_hole_product, (i + 1) * (j + 1))
+          .OnlyEnforceIf(var);
+    }
+    model_builder.AddEquality(count_per_pidgeon, 1);
+  }
+  for (const auto& count : count_per_hole) {
+    model_builder.AddLessOrEqual(count, 1);
+  }
+  Model model;
+  SatParameters params = GetParams();
+  model.Add(NewSatParameters(params));
+
+  CpSolverResponse response = SolveCpModel(model_builder.Build(), &model);
+
+  EXPECT_EQ(model.GetOrCreate<SharedTreeManager>()->NumWorkers(), 4);
+  EXPECT_EQ(response.status(), INFEASIBLE);
+}
+
+TEST(SharedTreeManagerTest, SplitTest) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var);
+  Model model;
+  SatParameters params;
+  params.set_num_workers(4);
+  params.set_shared_tree_num_workers(4);
+  params.set_cp_model_presolve(false);
+  model.Add(NewSatParameters(params));
+  LoadVariables(model_builder.Build(), false, &model);
+  auto* shared_tree_manager = model.GetOrCreate<SharedTreeManager>();
+  ProtoTrail shared_trail;
+
+  shared_tree_manager->ProposeSplit(shared_trail, {-1, 0});
+
+  EXPECT_EQ(shared_trail.MaxLevel(), 1);
+}
+
+TEST(SharedTreeManagerTest, RestartTest) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var);
+  Model model;
+  SatParameters params;
+  params.set_num_workers(4);
+  params.set_shared_tree_num_workers(4);
+  params.set_cp_model_presolve(false);
+  model.Add(NewSatParameters(params));
+  LoadVariables(model_builder.Build(), false, &model);
+  auto* shared_tree_manager = model.GetOrCreate<SharedTreeManager>();
+  ProtoTrail shared_trail;
+
+  shared_tree_manager->ProposeSplit(shared_trail, {-1, 0});
+  shared_tree_manager->Restart();
+  shared_tree_manager->SyncTree(shared_trail);
+
+  EXPECT_EQ(shared_trail.MaxLevel(), 0);
+}
+
+TEST(SharedTreeManagerTest, RestartTestWithLevelZeroImplications) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var);
+  Model model;
+  SatParameters params;
+  params.set_num_workers(4);
+  params.set_shared_tree_num_workers(4);
+  params.set_cp_model_presolve(false);
+  model.Add(NewSatParameters(params));
+  LoadVariables(model_builder.Build(), false, &model);
+  auto* shared_tree_manager = model.GetOrCreate<SharedTreeManager>();
+  ProtoTrail shared_trail;
+
+  shared_tree_manager->ProposeSplit(shared_trail, {-1, 0});
+  shared_tree_manager->CloseTree(shared_trail, 1);
+  shared_tree_manager->SyncTree(shared_trail);
+  shared_tree_manager->ReplaceTree(shared_trail);
+  shared_tree_manager->Restart();
+  shared_tree_manager->SyncTree(shared_trail);
+
+  EXPECT_EQ(shared_trail.NodeIds(0).size(), 0);
+  EXPECT_EQ(shared_trail.MaxLevel(), 0);
+}
+
+TEST(SharedTreeManagerTest, SharedBranchingTest) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var);
+  Model model;
+  SatParameters params;
+  params.set_num_workers(2);
+  params.set_shared_tree_num_workers(2);
+  params.set_cp_model_presolve(false);
+  model.Add(NewSatParameters(params));
+  LoadVariables(model_builder.Build(), false, &model);
+  auto* shared_tree_manager = model.GetOrCreate<SharedTreeManager>();
+  ProtoTrail trail1, trail2;
+
+  shared_tree_manager->ProposeSplit(trail1, {-1, 0});
+  shared_tree_manager->ReplaceTree(trail2);
+
+  EXPECT_EQ(trail1.MaxLevel(), 1);
+  EXPECT_EQ(trail2.MaxLevel(), 1);
+  EXPECT_EQ(trail1.Decision(1), trail2.Decision(1).Negated());
+}
+
+TEST(SharedTreeManagerTest, ObjectiveLbSplitTest) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var);
+  Model model;
+  SatParameters params;
+  params.set_num_workers(4);
+  params.set_shared_tree_num_workers(4);
+  params.set_cp_model_presolve(false);
+  params.set_shared_tree_split_strategy(
+      SatParameters::SPLIT_STRATEGY_OBJECTIVE_LB);
+  model.Add(NewSatParameters(params));
+  LoadVariables(model_builder.Build(), false, &model);
+  auto* response_manager = model.GetOrCreate<SharedResponseManager>();
+  response_manager->InitializeObjective(model_builder.Build());
+  auto* shared_tree_manager = model.GetOrCreate<SharedTreeManager>();
+  ProtoTrail trail1, trail2;
+
+  shared_tree_manager->ProposeSplit(trail1, {-1, 0});
+  ASSERT_EQ(trail1.MaxLevel(), 1);
+  trail1.SetObjectiveLb(1, 2);
+  shared_tree_manager->SyncTree(trail1);
+  shared_tree_manager->ReplaceTree(trail2);
+  ASSERT_EQ(trail2.MaxLevel(), 1);
+  trail2.SetObjectiveLb(1, 1);
+  shared_tree_manager->SyncTree(trail2);
+  // Reject this split because it is not at the global lower bound.
+  shared_tree_manager->ProposeSplit(trail1, {int_var.index(), 3});
+
+  EXPECT_EQ(response_manager->GetInnerObjectiveLowerBound(), 1);
+  EXPECT_EQ(shared_tree_manager->NumNodes(), 3);
+}
+
+TEST(SharedTreeManagerTest, DiscrepancySplitTestOneLeafPerWorker) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var);
+  Model model;
+  SatParameters params;
+  params.set_num_workers(4);
+  params.set_shared_tree_num_workers(4);
+  params.set_shared_tree_open_leaves_per_worker(1);
+  params.set_cp_model_presolve(false);
+  params.set_shared_tree_split_strategy(
+      SatParameters::SPLIT_STRATEGY_DISCREPANCY);
+  model.Add(NewSatParameters(params));
+  LoadVariables(model_builder.Build(), false, &model);
+  auto* response_manager = model.GetOrCreate<SharedResponseManager>();
+  response_manager->InitializeObjective(model_builder.Build());
+  auto* shared_tree_manager = model.GetOrCreate<SharedTreeManager>();
+  ProtoTrail trail1, trail2;
+
+  shared_tree_manager->ProposeSplit(trail1, {-1, 0});
+  shared_tree_manager->SyncTree(trail1);
+  shared_tree_manager->ReplaceTree(trail2);
+  shared_tree_manager->ProposeSplit(trail2, {int_var.index(), 3});
+  shared_tree_manager->ProposeSplit(trail1, {int_var.index(), 3});
+  // Reject this split: 2 depth + 1 discrepancy is not minimal.
+  shared_tree_manager->ProposeSplit(trail2, {int_var.index(), 5});
+  // Reject this split: 2 depth  + 0 discrepancy is not minimal.
+  shared_tree_manager->ProposeSplit(trail1, {int_var.index(), 5});
+
+  EXPECT_EQ(trail1.MaxLevel(), 2);
+  EXPECT_EQ(trail2.MaxLevel(), 2);
+  EXPECT_EQ(shared_tree_manager->NumNodes(), 7);
+}
+
+TEST(SharedTreeManagerTest, DiscrepancySplitTest) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var);
+  Model model;
+  SatParameters params;
+  params.set_num_workers(2);
+  params.set_shared_tree_num_workers(2);
+  params.set_shared_tree_open_leaves_per_worker(2);
+  params.set_cp_model_presolve(false);
+  params.set_shared_tree_split_strategy(
+      SatParameters::SPLIT_STRATEGY_DISCREPANCY);
+  model.Add(NewSatParameters(params));
+  LoadVariables(model_builder.Build(), false, &model);
+  auto* response_manager = model.GetOrCreate<SharedResponseManager>();
+  response_manager->InitializeObjective(model_builder.Build());
+  auto* shared_tree_manager = model.GetOrCreate<SharedTreeManager>();
+  ProtoTrail trail1, trail2;
+
+  shared_tree_manager->ProposeSplit(trail1, {-1, 0});
+  shared_tree_manager->SyncTree(trail1);
+  shared_tree_manager->ReplaceTree(trail2);
+  shared_tree_manager->ProposeSplit(trail2, {int_var.index(), 3});
+  shared_tree_manager->ProposeSplit(trail1, {int_var.index(), 3});
+  // Reject this split: 2 depth + 1 discrepancy is not minimal.
+  shared_tree_manager->ProposeSplit(trail2, {int_var.index(), 5});
+  // Reject this split: 2 depth  + 0 discrepancy is not minimal.
+  shared_tree_manager->ProposeSplit(trail1, {int_var.index(), 5});
+
+  EXPECT_EQ(trail1.MaxLevel(), 2);
+  EXPECT_EQ(trail2.MaxLevel(), 2);
+  EXPECT_EQ(shared_tree_manager->NumNodes(), 7);
+}
+
+TEST(SharedTreeManagerTest, BalancedSplitTestOneLeafPerWorker) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var);
+  Model model;
+  SatParameters params;
+  params.set_num_workers(5);
+  params.set_shared_tree_num_workers(5);
+  params.set_shared_tree_open_leaves_per_worker(1);
+  params.set_cp_model_presolve(false);
+  params.set_shared_tree_split_strategy(
+      SatParameters::SPLIT_STRATEGY_BALANCED_TREE);
+  model.Add(NewSatParameters(params));
+  LoadVariables(model_builder.Build(), false, &model);
+  auto* response_manager = model.GetOrCreate<SharedResponseManager>();
+  response_manager->InitializeObjective(model_builder.Build());
+  auto* shared_tree_manager = model.GetOrCreate<SharedTreeManager>();
+  ProtoTrail trail1, trail2;
+
+  shared_tree_manager->ProposeSplit(trail1, {-1, 0});
+  shared_tree_manager->SyncTree(trail1);
+  shared_tree_manager->ReplaceTree(trail2);
+  shared_tree_manager->SyncTree(trail2);
+  shared_tree_manager->ProposeSplit(trail1, {int_var.index(), 3});
+  // Reject this split because it creates an unbalanced tree
+  shared_tree_manager->ProposeSplit(trail1, {int_var.index(), 5});
+  shared_tree_manager->ProposeSplit(trail2, {int_var.index(), 3});
+
+  EXPECT_EQ(shared_tree_manager->NumNodes(), 7);
+  EXPECT_EQ(trail1.MaxLevel(), 2);
+  EXPECT_EQ(trail2.MaxLevel(), 2);
+}
+
+TEST(SharedTreeManagerTest, BalancedSplitTest) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var);
+  Model model;
+  SatParameters params;
+  params.set_num_workers(3);
+  params.set_shared_tree_num_workers(3);
+  params.set_shared_tree_open_leaves_per_worker(2);
+  params.set_cp_model_presolve(false);
+  params.set_shared_tree_split_strategy(
+      SatParameters::SPLIT_STRATEGY_BALANCED_TREE);
+  model.Add(NewSatParameters(params));
+  LoadVariables(model_builder.Build(), false, &model);
+  auto* response_manager = model.GetOrCreate<SharedResponseManager>();
+  response_manager->InitializeObjective(model_builder.Build());
+  auto* shared_tree_manager = model.GetOrCreate<SharedTreeManager>();
+  ProtoTrail trail1, trail2;
+
+  shared_tree_manager->ProposeSplit(trail1, {-1, 0});
+  shared_tree_manager->SyncTree(trail1);
+  shared_tree_manager->ReplaceTree(trail2);
+  shared_tree_manager->SyncTree(trail2);
+  shared_tree_manager->ProposeSplit(trail1, {int_var.index(), 3});
+  // Reject this split because it creates an unbalanced tree
+  shared_tree_manager->ProposeSplit(trail1, {int_var.index(), 5});
+  shared_tree_manager->ProposeSplit(trail2, {int_var.index(), 3});
+
+  EXPECT_EQ(shared_tree_manager->NumNodes(), 7);
+  EXPECT_EQ(trail1.MaxLevel(), 2);
+  EXPECT_EQ(trail2.MaxLevel(), 2);
+}
+
+TEST(SharedTreeManagerTest, CloseTreeTest) {
+  CpModelBuilder model_builder;
+  auto bool_var = model_builder.NewBoolVar();
+  auto int_var = model_builder.NewIntVar({0, 7});
+  model_builder.AddLessOrEqual(int_var, 3).OnlyEnforceIf(bool_var);
+  model_builder.Maximize(int_var);
+  Model model;
+  SatParameters params;
+  params.set_num_workers(4);
+  params.set_shared_tree_num_workers(4);
+  params.set_cp_model_presolve(false);
+  model.Add(NewSatParameters(params));
+  LoadVariables(model_builder.Build(), false, &model);
+  auto* shared_tree_manager = model.GetOrCreate<SharedTreeManager>();
+  ProtoTrail trail1, trail2, trail3;
+  shared_tree_manager->ProposeSplit(trail1, {-1, 0});
+  shared_tree_manager->ReplaceTree(trail2);
+  shared_tree_manager->ProposeSplit(trail1, {1, 0});
+  shared_tree_manager->CloseTree(trail1, 1);
+  shared_tree_manager->ReplaceTree(trail1);
+
+  EXPECT_EQ(trail1.MaxLevel(), 0);
+  EXPECT_EQ(trail2.MaxLevel(), 1);
+  EXPECT_EQ(trail2.Decision(1), ProtoLiteral(0, 1));
+}
+// TODO(user): Test objective propagation.
+}  // namespace
+}  // namespace sat
+}  // namespace operations_research