From ed9b82dc919991daf57a7983eef56d9a2637c252 Mon Sep 17 00:00:00 2001
From: Laurent Perron <lperron@google.com>
Date: Thu, 15 Feb 2024 08:46:25 +0100
Subject: [PATCH] pdlp uses SolverLogger

---
 ortools/pdlp/BUILD.bazel                      |   2 +
 ortools/pdlp/iteration_stats.cc               |  91 ++--
 ortools/pdlp/iteration_stats.h                |   7 +-
 ortools/pdlp/iteration_stats_test.cc          |   1 +
 ortools/pdlp/primal_dual_hybrid_gradient.cc   | 515 +++++++++++-------
 ortools/pdlp/primal_dual_hybrid_gradient.h    |   9 +
 .../pdlp/primal_dual_hybrid_gradient_test.cc  |  66 ++-
 ortools/pdlp/sharded_optimization_utils.cc    |  24 +-
 ortools/pdlp/sharded_optimization_utils.h     |   5 +-
 .../pdlp/sharded_optimization_utils_test.cc   |   8 +
 ortools/pdlp/sharded_quadratic_program.cc     |  46 +-
 ortools/pdlp/sharded_quadratic_program.h      |   6 +-
 ortools/pdlp/solve_log.proto                  |  23 +-
 13 files changed, 524 insertions(+), 279 deletions(-)

diff --git a/ortools/pdlp/BUILD.bazel b/ortools/pdlp/BUILD.bazel
index 9c74843487..3034d6b01a 100644
--- a/ortools/pdlp/BUILD.bazel
+++ b/ortools/pdlp/BUILD.bazel
@@ -117,6 +117,7 @@ cc_library(
         "//ortools/lp_data",
         "//ortools/lp_data:base",
         "//ortools/lp_data:proto_utils",
+        "//ortools/util:logging",
         "@com_google_absl//absl/algorithm:container",
         "@com_google_absl//absl/status",
         "@com_google_absl//absl/status:statusor",
@@ -252,6 +253,7 @@ cc_library(
         ":sharder",
         "//ortools/base",
         "//ortools/base:threadpool",
+        "//ortools/util:logging",
         "@com_google_absl//absl/memory",
         "@com_google_absl//absl/strings",
         "@eigen//:eigen3",
diff --git a/ortools/pdlp/iteration_stats.cc b/ortools/pdlp/iteration_stats.cc
index b8698ac819..ee94eb0e93 100644
--- a/ortools/pdlp/iteration_stats.cc
+++ b/ortools/pdlp/iteration_stats.cc
@@ -60,8 +60,7 @@ struct ResidualNorms {
 // any corrections (so the returned `.objective_correction == 0` and
 // `.objective_full_correction == 0`). `sharded_qp` is assumed to be the scaled
 // problem. If `use_homogeneous_constraint_bounds` is set to true the residuals
-// are computed with upper and lower bounds zeroed out (note that we only zero
-// out the bounds that are finite in the original problem).
+// are computed with all finite bounds mapped to zero.
 // NOTE: `componentwise_residual_offset` only affects the value of
 // `l_inf_componentwise_residual` in the returned `ResidualNorms`.
 ResidualNorms PrimalResidualNorms(
@@ -221,13 +220,13 @@ ResidualNorms DualResidualNorms(const PrimalDualHybridGradientParams& params,
           // unscaled_primal_gradient = primal_gradient / scale, so the scales
           // cancel out.
           if (primal_gradient_shard[i] == 0.0) continue;
-          const double bound_for_rc = primal_gradient_shard[i] > 0.0
-                                          ? lower_bound_shard[i]
-                                          : upper_bound_shard[i];
+          const double upper_bound = upper_bound_shard[i];
+          const double lower_bound = lower_bound_shard[i];
+          const double bound_for_rc =
+              primal_gradient_shard[i] > 0.0 ? lower_bound : upper_bound;
           dual_full_correction += bound_for_rc * primal_gradient_shard[i];
           VariableBounds effective_bounds = EffectiveVariableBounds(
-              params, primal_solution_shard[i], lower_bound_shard[i],
-              upper_bound_shard[i]);
+              params, primal_solution_shard[i], lower_bound, upper_bound);
           // The dual correction (using the appropriate bound) is applied even
           // if the gradient is handled as a residual, so that the dual
           // objective is convex.
@@ -453,13 +452,46 @@ ConvergenceInformation ComputeConvergenceInformation(
   return result;
 }
 
+namespace {
+
+double PrimalRayMaxSignViolation(const ShardedQuadraticProgram& sharded_qp,
+                                 const VectorXd& col_scaling_vec,
+                                 const VectorXd& scaled_primal_ray) {
+  VectorXd primal_ray_local_max_sign_violation(
+      sharded_qp.PrimalSharder().NumShards());
+  sharded_qp.PrimalSharder().ParallelForEachShard(
+      [&](const Sharder::Shard& shard) {
+        const auto lower_bound_shard =
+            shard(sharded_qp.Qp().variable_lower_bounds);
+        const auto upper_bound_shard =
+            shard(sharded_qp.Qp().variable_upper_bounds);
+        const auto ray_shard = shard(scaled_primal_ray);
+        const auto scale_shard = shard(col_scaling_vec);
+        double local_max = 0.0;
+        for (int64_t i = 0; i < ray_shard.size(); ++i) {
+          if (std::isfinite(lower_bound_shard[i])) {
+            local_max = std::max(local_max, -ray_shard[i] * scale_shard[i]);
+          }
+          if (std::isfinite(upper_bound_shard[i])) {
+            local_max = std::max(local_max, ray_shard[i] * scale_shard[i]);
+          }
+        }
+        primal_ray_local_max_sign_violation[shard.Index()] = local_max;
+      });
+  return primal_ray_local_max_sign_violation.lpNorm<Eigen::Infinity>();
+}
+
+}  // namespace
+
 InfeasibilityInformation ComputeInfeasibilityInformation(
     const PrimalDualHybridGradientParams& params,
     const ShardedQuadraticProgram& scaled_sharded_qp,
     const Eigen::VectorXd& col_scaling_vec,
     const Eigen::VectorXd& row_scaling_vec,
     const Eigen::VectorXd& scaled_primal_ray,
-    const Eigen::VectorXd& scaled_dual_ray, PointType candidate_type) {
+    const Eigen::VectorXd& scaled_dual_ray,
+    const Eigen::VectorXd& primal_solution_for_residual_tests,
+    PointType candidate_type) {
   const QuadraticProgram& qp = scaled_sharded_qp.Qp();
   CHECK_EQ(col_scaling_vec.size(), scaled_sharded_qp.PrimalSize());
   CHECK_EQ(row_scaling_vec.size(), scaled_sharded_qp.DualSize());
@@ -479,8 +511,9 @@ InfeasibilityInformation ComputeInfeasibilityInformation(
   // We don't use `dual_residuals.l_inf_componentwise_residual`, so don't need
   // to set `componentwise_residual_offset` to a meaningful value.
   ResidualNorms dual_residuals = DualResidualNorms(
-      params, scaled_sharded_qp, col_scaling_vec, scaled_primal_ray,
-      scaled_primal_gradient, /*componentwise_residual_offset=*/0.0);
+      params, scaled_sharded_qp, col_scaling_vec,
+      primal_solution_for_residual_tests, scaled_primal_gradient,
+      /*componentwise_residual_offset=*/0.0);
 
   double dual_ray_objective =
       DualObjectiveBoundsTerm(scaled_sharded_qp, scaled_dual_ray) +
@@ -501,32 +534,12 @@ InfeasibilityInformation ComputeInfeasibilityInformation(
       PrimalResidualNorms(scaled_sharded_qp, row_scaling_vec, scaled_primal_ray,
                           /*componentwise_residual_offset=*/0.0,
                           /*use_homogeneous_constraint_bounds=*/true);
-  // `primal_residuals` contains the violations of the linear constraints. The
-  // signs of the components are also constrained by the presence or absence
-  // of variable bounds.
-  VectorXd primal_ray_local_sign_max_violation(
-      scaled_sharded_qp.PrimalSharder().NumShards());
-  scaled_sharded_qp.PrimalSharder().ParallelForEachShard(
-      [&](const Sharder::Shard& shard) {
-        const auto lower_bound_shard =
-            shard(scaled_sharded_qp.Qp().variable_lower_bounds);
-        const auto upper_bound_shard =
-            shard(scaled_sharded_qp.Qp().variable_upper_bounds);
-        const auto ray_shard = shard(scaled_primal_ray);
-        const auto scale_shard = shard(col_scaling_vec);
-        double local_max = 0.0;
-        for (int64_t i = 0; i < ray_shard.size(); ++i) {
-          if (std::isfinite(lower_bound_shard[i])) {
-            local_max = std::max(local_max, -ray_shard[i] * scale_shard[i]);
-          }
-          if (std::isfinite(upper_bound_shard[i])) {
-            local_max = std::max(local_max, ray_shard[i] * scale_shard[i]);
-          }
-        }
-        primal_ray_local_sign_max_violation[shard.Index()] = local_max;
-      });
-  const double primal_ray_sign_max_violation =
-      primal_ray_local_sign_max_violation.lpNorm<Eigen::Infinity>();
+
+  // The primal ray should have been projected onto the feasibility bounds, so
+  // that it has the correct signs.
+  DCHECK_EQ(PrimalRayMaxSignViolation(scaled_sharded_qp, col_scaling_vec,
+                                      scaled_primal_ray),
+            0.0);
 
   if (l_inf_primal > 0.0) {
     VectorXd scaled_objective_product =
@@ -534,10 +547,8 @@ InfeasibilityInformation ComputeInfeasibilityInformation(
     result.set_primal_ray_quadratic_norm(
         LInfNorm(scaled_objective_product, scaled_sharded_qp.PrimalSharder()) /
         l_inf_primal);
-    result.set_max_primal_ray_infeasibility(
-        std::max(primal_residuals.l_inf_residual,
-                 primal_ray_sign_max_violation) /
-        l_inf_primal);
+    result.set_max_primal_ray_infeasibility(primal_residuals.l_inf_residual /
+                                            l_inf_primal);
     result.set_primal_ray_linear_objective(
         Dot(scaled_primal_ray, qp.objective_vector,
             scaled_sharded_qp.PrimalSharder()) /
diff --git a/ortools/pdlp/iteration_stats.h b/ortools/pdlp/iteration_stats.h
index cc9672d64c..978ae829a6 100644
--- a/ortools/pdlp/iteration_stats.h
+++ b/ortools/pdlp/iteration_stats.h
@@ -65,13 +65,18 @@ ConvergenceInformation ComputeConvergenceInformation(
 // row_scaling_vec)`. `scaled_primal_ray` and `scaled_dual_ray` are
 // potential certificates for the scaled problem. The stats are computed with
 // respect to the implicit original problem.
+// `primal_solution_for_residual_tests` is used instead of `scaled_primal_ray`
+// when deciding whether or not to treat a primal gradient as a dual residual
+// or not.
 InfeasibilityInformation ComputeInfeasibilityInformation(
     const PrimalDualHybridGradientParams& params,
     const ShardedQuadraticProgram& scaled_sharded_qp,
     const Eigen::VectorXd& col_scaling_vec,
     const Eigen::VectorXd& row_scaling_vec,
     const Eigen::VectorXd& scaled_primal_ray,
-    const Eigen::VectorXd& scaled_dual_ray, PointType candidate_type);
+    const Eigen::VectorXd& scaled_dual_ray,
+    const Eigen::VectorXd& primal_solution_for_residual_tests,
+    PointType candidate_type);
 
 // Computes the reduced costs vector, objective_matrix * `primal_solution` +
 // objective_vector - constraint_matrix * `dual_solution`, when
diff --git a/ortools/pdlp/iteration_stats_test.cc b/ortools/pdlp/iteration_stats_test.cc
index 6554a17aeb..4dff93ae64 100644
--- a/ortools/pdlp/iteration_stats_test.cc
+++ b/ortools/pdlp/iteration_stats_test.cc
@@ -14,6 +14,7 @@
 #include "ortools/pdlp/iteration_stats.h"
 
 #include <cmath>
+#include <limits>
 #include <optional>
 #include <utility>
 
diff --git a/ortools/pdlp/primal_dual_hybrid_gradient.cc b/ortools/pdlp/primal_dual_hybrid_gradient.cc
index bb9bc84817..05af624274 100644
--- a/ortools/pdlp/primal_dual_hybrid_gradient.cc
+++ b/ortools/pdlp/primal_dual_hybrid_gradient.cc
@@ -56,12 +56,14 @@
 #include "ortools/pdlp/solvers_proto_validation.h"
 #include "ortools/pdlp/termination.h"
 #include "ortools/pdlp/trust_region.h"
+#include "ortools/util/logging.h"
 
 namespace operations_research::pdlp {
 
 namespace {
 
 using ::Eigen::VectorXd;
+using ::operations_research::SolverLogger;
 
 using IterationStatsCallback =
     std::function<void(const IterationCallbackInfo&)>;
@@ -69,7 +71,7 @@ using IterationStatsCallback =
 // Computes a `num_threads` that is capped by the problem size and `num_shards`,
 // if specified, to avoid creating unusable threads.
 int NumThreads(const int num_threads, const int num_shards,
-               const QuadraticProgram& qp) {
+               const QuadraticProgram& qp, SolverLogger& logger) {
   int capped_num_threads = num_threads;
   if (num_shards > 0) {
     capped_num_threads = std::min(capped_num_threads, num_shards);
@@ -80,9 +82,9 @@ int NumThreads(const int num_threads, const int num_shards,
       static_cast<int>(std::min(int64_t{capped_num_threads}, problem_limit));
   capped_num_threads = std::max(capped_num_threads, 1);
   if (capped_num_threads != num_threads) {
-    LOG(WARNING) << "Reducing num_threads from " << num_threads << " to "
-                 << capped_num_threads
-                 << " because additional threads would be useless.";
+    SOLVER_LOG(&logger, "WARNING: Reducing num_threads from ", num_threads,
+               " to ", capped_num_threads,
+               " because additional threads would be useless.");
   }
   return capped_num_threads;
 }
@@ -184,10 +186,6 @@ std::string ConvergenceInformationShortString(
   LOG(FATAL) << "Invalid residual norm " << residual_norm << ".";
 }
 
-void LogInfoWithoutPrefix(absl::string_view message) {
-  LOG(INFO).NoPrefix() << message;
-}
-
 // Logs a description of `iter_stats`, based on the
 // `iter_stats.convergence_information` entry with
 // `.candidate_type()==preferred_candidate` if one exists, otherwise based on
@@ -198,7 +196,7 @@ void LogIterationStats(int verbosity_level, bool use_feasibility_polishing,
                        const IterationStats& iter_stats,
                        const TerminationCriteria& termination_criteria,
                        const QuadraticProgramBoundNorms& bound_norms,
-                       PointType preferred_candidate) {
+                       PointType preferred_candidate, SolverLogger& logger) {
   std::string iteration_string =
       verbosity_level >= 3
           ? absl::StrFormat("%6d %8.1f %6.1f", iter_stats.iteration_number(),
@@ -261,18 +259,18 @@ void LogIterationStats(int verbosity_level, bool use_feasibility_polishing,
             : ConvergenceInformationShortString(
                   *convergence_information, relative_information,
                   termination_criteria.optimality_norm());
-    LogInfoWithoutPrefix(absl::StrCat(phase_string, iterate_string,
-                                      iteration_string, " | ",
-                                      convergence_string));
+    SOLVER_LOG(&logger, phase_string, iterate_string, iteration_string, " | ",
+               convergence_string);
   } else {
     // No convergence information, just log the basic work stats.
-    LogInfoWithoutPrefix(absl::StrCat(
-        phase_string, verbosity_level >= 4 ? "? " : "", iteration_string));
+    SOLVER_LOG(&logger, phase_string, verbosity_level >= 4 ? "? " : "",
+               iteration_string);
   }
 }
 
 void LogIterationStatsHeader(int verbosity_level,
-                             bool use_feasibility_polishing) {
+                             bool use_feasibility_polishing,
+                             SolverLogger& logger) {
   std::string work_labels =
       verbosity_level >= 3
           ? absl::StrFormat("%6s %8s %6s", "iter#", "kkt_pass", "time")
@@ -286,9 +284,9 @@ void LogIterationStatsHeader(int verbosity_level,
                 "dual_var_l2")
           : absl::StrFormat("%10s %10s %10s | %10s %10s", "rel_p_res",
                             "rel_d_res", "rel_gap", "prim_obj", "dual_obj");
-  LogInfoWithoutPrefix(absl::StrCat(use_feasibility_polishing ? "f " : "",
-                                    verbosity_level >= 4 ? "I " : "",
-                                    work_labels, " | ", convergence_labels));
+  SOLVER_LOG(&logger, use_feasibility_polishing ? "f " : "",
+             verbosity_level >= 4 ? "I " : "", work_labels, " | ",
+             convergence_labels);
 }
 
 enum class InnerStepOutcome {
@@ -319,12 +317,14 @@ class PreprocessSolver {
  public:
   // Assumes that `qp` and `params` are valid.
   // Note that the `qp` is intentionally passed by value.
+  // `logger` is captured, and must outlive the class.
   // NOTE: Many `PreprocessSolver` methods accept a `params` argument. This is
   // passed as an argument instead of stored as a member variable to support
   // using different `params` in different contexts with the same
   // `PreprocessSolver` object.
   explicit PreprocessSolver(QuadraticProgram qp,
-                            const PrimalDualHybridGradientParams& params);
+                            const PrimalDualHybridGradientParams& params,
+                            SolverLogger* logger);
 
   // Not copyable or movable (because `glop::MainLpPreprocessor` isn't).
   PreprocessSolver(const PreprocessSolver&) = delete;
@@ -383,7 +383,8 @@ class PreprocessSolver {
   // NOTE: `result` is passed by value. To avoid unnecessary copying, move this
   // object (i.e. use `std::move()`).
   SolverResult ConstructOriginalSolverResult(
-      const PrimalDualHybridGradientParams& params, SolverResult result) const;
+      const PrimalDualHybridGradientParams& params, SolverResult result,
+      SolverLogger& logger) const;
 
   const ShardedQuadraticProgram& ShardedWorkingQp() const {
     return sharded_qp_;
@@ -416,6 +417,8 @@ class PreprocessSolver {
     return original_bound_norms_;
   }
 
+  SolverLogger& Logger() { return logger_; }
+
  private:
   struct PresolveInfo {
     explicit PresolveInfo(ShardedQuadraticProgram original_qp,
@@ -453,7 +456,7 @@ class PreprocessSolver {
                                 VectorXd& starting_primal_solution,
                                 VectorXd& starting_dual_solution);
 
-  void LogQuadraticProgramStats(const QuadraticProgramStats& stats);
+  void LogQuadraticProgramStats(const QuadraticProgramStats& stats) const;
 
   double InitialPrimalWeight(const PrimalDualHybridGradientParams& params,
                              double l2_norm_primal_linear_objective,
@@ -498,6 +501,7 @@ class PreprocessSolver {
   // A counter used to trigger writing iteration headers.
   int log_counter_ = 0;
   absl::Time time_of_last_log_ = absl::InfinitePast();
+  SolverLogger& logger_;
   IterationStatsCallback iteration_stats_callback_;
 };
 
@@ -709,18 +713,22 @@ class Solver {
 };
 
 PreprocessSolver::PreprocessSolver(QuadraticProgram qp,
-                                   const PrimalDualHybridGradientParams& params)
-    : num_threads_(NumThreads(params.num_threads(), params.num_shards(), qp)),
+                                   const PrimalDualHybridGradientParams& params,
+                                   SolverLogger* logger)
+    : num_threads_(
+          NumThreads(params.num_threads(), params.num_shards(), qp, *logger)),
       num_shards_(NumShards(num_threads_, params.num_shards())),
-      sharded_qp_(std::move(qp), num_threads_, num_shards_) {}
+      sharded_qp_(std::move(qp), num_threads_, num_shards_),
+      logger_(*logger) {}
 
 SolverResult ErrorSolverResult(const TerminationReason reason,
-                               const std::string& message) {
+                               const std::string& message,
+                               SolverLogger& logger) {
   SolveLog error_log;
   error_log.set_termination_reason(reason);
   error_log.set_termination_string(message);
-  LOG(WARNING) << "The solver did not run because of invalid input: "
-               << message;
+  SOLVER_LOG(&logger,
+             "The solver did not run because of invalid input: ", message);
   return SolverResult{.solve_log = error_log};
 }
 
@@ -730,28 +738,29 @@ SolverResult ErrorSolverResult(const TerminationReason reason,
 // but that cause problems for other code such as glop's presolve.
 std::optional<SolverResult> CheckProblemStats(
     const QuadraticProgramStats& problem_stats, const double objective_offset,
-    bool check_excessively_small_values) {
+    bool check_excessively_small_values, SolverLogger& logger) {
   const double kExcessiveInputValue = 1e50;
   const double kExcessivelySmallInputValue = 1e-50;
   const double kMaxDynamicRange = 1e20;
   if (std::isnan(problem_stats.constraint_matrix_l2_norm())) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
-                             "Constraint matrix has a NAN.");
+                             "Constraint matrix has a NAN.", logger);
   }
   if (problem_stats.constraint_matrix_abs_max() > kExcessiveInputValue) {
     return ErrorSolverResult(
         TERMINATION_REASON_INVALID_PROBLEM,
         absl::StrCat("Constraint matrix has a non-zero with absolute value ",
                      problem_stats.constraint_matrix_abs_max(),
-                     " which exceeds limit of ", kExcessiveInputValue, "."));
+                     " which exceeds limit of ", kExcessiveInputValue, "."),
+        logger);
   }
   if (problem_stats.constraint_matrix_abs_max() >
       kMaxDynamicRange * problem_stats.constraint_matrix_abs_min()) {
-    LOG(WARNING) << "Constraint matrix has largest absolute value "
-                 << problem_stats.constraint_matrix_abs_max()
-                 << " and smallest non-zero absolute value "
-                 << problem_stats.constraint_matrix_abs_min()
-                 << " performance may suffer.";
+    SOLVER_LOG(
+        &logger, "WARNING: Constraint matrix has largest absolute value ",
+        problem_stats.constraint_matrix_abs_max(),
+        " and smallest non-zero absolute value ",
+        problem_stats.constraint_matrix_abs_min(), " performance may suffer.");
   }
   if (problem_stats.constraint_matrix_col_min_l_inf_norm() > 0 &&
       problem_stats.constraint_matrix_col_min_l_inf_norm() <
@@ -761,7 +770,8 @@ std::optional<SolverResult> CheckProblemStats(
         absl::StrCat("Constraint matrix has a column with Linf norm ",
                      problem_stats.constraint_matrix_col_min_l_inf_norm(),
                      " which is less than limit of ",
-                     kExcessivelySmallInputValue, "."));
+                     kExcessivelySmallInputValue, "."),
+        logger);
   }
   if (problem_stats.constraint_matrix_row_min_l_inf_norm() > 0 &&
       problem_stats.constraint_matrix_row_min_l_inf_norm() <
@@ -771,11 +781,12 @@ std::optional<SolverResult> CheckProblemStats(
         absl::StrCat("Constraint matrix has a row with Linf norm ",
                      problem_stats.constraint_matrix_row_min_l_inf_norm(),
                      " which is less than limit of ",
-                     kExcessivelySmallInputValue, "."));
+                     kExcessivelySmallInputValue, "."),
+        logger);
   }
   if (std::isnan(problem_stats.combined_bounds_l2_norm())) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
-                             "Constraint bounds vector has a NAN.");
+                             "Constraint bounds vector has a NAN.", logger);
   }
   if (problem_stats.combined_bounds_max() > kExcessiveInputValue) {
     return ErrorSolverResult(
@@ -783,7 +794,8 @@ std::optional<SolverResult> CheckProblemStats(
         absl::StrCat("Combined constraint bounds vector has a non-zero with "
                      "absolute value ",
                      problem_stats.combined_bounds_max(),
-                     " which exceeds limit of ", kExcessiveInputValue, "."));
+                     " which exceeds limit of ", kExcessiveInputValue, "."),
+        logger);
   }
   if (check_excessively_small_values &&
       problem_stats.combined_bounds_min() > 0 &&
@@ -794,19 +806,22 @@ std::optional<SolverResult> CheckProblemStats(
                      "absolute value ",
                      problem_stats.combined_bounds_min(),
                      " which is less than the limit of ",
-                     kExcessivelySmallInputValue, "."));
+                     kExcessivelySmallInputValue, "."),
+        logger);
   }
   if (problem_stats.combined_bounds_max() >
       kMaxDynamicRange * problem_stats.combined_bounds_min()) {
-    LOG(WARNING)
-        << "Combined constraint bounds vector has largest absolute value "
-        << problem_stats.combined_bounds_max()
-        << " and smallest non-zero absolute value "
-        << problem_stats.combined_bounds_min() << "; performance may suffer.";
+    SOLVER_LOG(&logger,
+               "WARNING: Combined constraint bounds vector has largest "
+               "absolute value ",
+               problem_stats.combined_bounds_max(),
+               " and smallest non-zero absolute value ",
+               problem_stats.combined_bounds_min(),
+               "; performance may suffer.");
   }
   if (std::isnan(problem_stats.combined_variable_bounds_l2_norm())) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
-                             "Variable bounds vector has a NAN.");
+                             "Variable bounds vector has a NAN.", logger);
   }
   if (problem_stats.combined_variable_bounds_max() > kExcessiveInputValue) {
     return ErrorSolverResult(
@@ -814,7 +829,8 @@ std::optional<SolverResult> CheckProblemStats(
         absl::StrCat("Combined variable bounds vector has a non-zero with "
                      "absolute value ",
                      problem_stats.combined_variable_bounds_max(),
-                     " which exceeds limit of ", kExcessiveInputValue, "."));
+                     " which exceeds limit of ", kExcessiveInputValue, "."),
+        logger);
   }
   if (check_excessively_small_values &&
       problem_stats.combined_variable_bounds_min() > 0 &&
@@ -826,46 +842,51 @@ std::optional<SolverResult> CheckProblemStats(
                      "absolute value ",
                      problem_stats.combined_variable_bounds_min(),
                      " which is less than the limit of ",
-                     kExcessivelySmallInputValue, "."));
+                     kExcessivelySmallInputValue, "."),
+        logger);
   }
   if (problem_stats.combined_variable_bounds_max() >
       kMaxDynamicRange * problem_stats.combined_variable_bounds_min()) {
-    LOG(WARNING)
-        << "Combined variable bounds vector has largest absolute value "
-        << problem_stats.combined_variable_bounds_max()
-        << " and smallest non-zero absolute value "
-        << problem_stats.combined_variable_bounds_min()
-        << "; performance may suffer.";
+    SOLVER_LOG(
+        &logger,
+        "WARNING: Combined variable bounds vector has largest absolute value ",
+        problem_stats.combined_variable_bounds_max(),
+        " and smallest non-zero absolute value ",
+        problem_stats.combined_variable_bounds_min(),
+        "; performance may suffer.");
   }
   if (problem_stats.variable_bound_gaps_max() >
       kMaxDynamicRange * problem_stats.variable_bound_gaps_min()) {
-    LOG(WARNING) << "Variable bound gap vector has largest absolute value "
-                 << problem_stats.variable_bound_gaps_max()
-                 << " and smallest non-zero absolute value "
-                 << problem_stats.variable_bound_gaps_min()
-                 << "; performance may suffer.";
+    SOLVER_LOG(&logger,
+               "WARNING: Variable bound gap vector has largest absolute value ",
+               problem_stats.variable_bound_gaps_max(),
+               " and smallest non-zero absolute value ",
+               problem_stats.variable_bound_gaps_min(),
+               "; performance may suffer.");
   }
   if (std::isnan(objective_offset)) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
-                             "Objective offset is NAN.");
+                             "Objective offset is NAN.", logger);
   }
   if (std::abs(objective_offset) > kExcessiveInputValue) {
     return ErrorSolverResult(
         TERMINATION_REASON_INVALID_PROBLEM,
         absl::StrCat("Objective offset ", objective_offset,
                      " has absolute value which exceeds limit of ",
-                     kExcessiveInputValue, "."));
+                     kExcessiveInputValue, "."),
+        logger);
   }
   if (std::isnan(problem_stats.objective_vector_l2_norm())) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
-                             "Objective vector has a NAN.");
+                             "Objective vector has a NAN.", logger);
   }
   if (problem_stats.objective_vector_abs_max() > kExcessiveInputValue) {
     return ErrorSolverResult(
         TERMINATION_REASON_INVALID_PROBLEM,
         absl::StrCat("Objective vector has a non-zero with absolute value ",
                      problem_stats.objective_vector_abs_max(),
-                     " which exceeds limit of ", kExcessiveInputValue, "."));
+                     " which exceeds limit of ", kExcessiveInputValue, "."),
+        logger);
   }
   if (check_excessively_small_values &&
       problem_stats.objective_vector_abs_min() > 0 &&
@@ -875,41 +896,43 @@ std::optional<SolverResult> CheckProblemStats(
         absl::StrCat("Objective vector has a non-zero with absolute value ",
                      problem_stats.objective_vector_abs_min(),
                      " which is less than the limit of ",
-                     kExcessivelySmallInputValue, "."));
+                     kExcessivelySmallInputValue, "."),
+        logger);
   }
   if (problem_stats.objective_vector_abs_max() >
       kMaxDynamicRange * problem_stats.objective_vector_abs_min()) {
-    LOG(WARNING) << "Objective vector has largest absolute value "
-                 << problem_stats.objective_vector_abs_max()
-                 << " and smallest non-zero absolute value "
-                 << problem_stats.objective_vector_abs_min()
-                 << "; performance may suffer.";
+    SOLVER_LOG(&logger, "WARNING: Objective vector has largest absolute value ",
+               problem_stats.objective_vector_abs_max(),
+               " and smallest non-zero absolute value ",
+               problem_stats.objective_vector_abs_min(),
+               "; performance may suffer.");
   }
   if (std::isnan(problem_stats.objective_matrix_l2_norm())) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
-                             "Objective matrix has a NAN.");
+                             "Objective matrix has a NAN.", logger);
   }
   if (problem_stats.objective_matrix_abs_max() > kExcessiveInputValue) {
     return ErrorSolverResult(
         TERMINATION_REASON_INVALID_PROBLEM,
         absl::StrCat("Objective matrix has a non-zero with absolute value ",
                      problem_stats.objective_matrix_abs_max(),
-                     " which exceeds limit of ", kExcessiveInputValue, "."));
+                     " which exceeds limit of ", kExcessiveInputValue, "."),
+        logger);
   }
   if (problem_stats.objective_matrix_abs_max() >
       kMaxDynamicRange * problem_stats.objective_matrix_abs_min()) {
-    LOG(WARNING) << "Objective matrix has largest absolute value "
-                 << problem_stats.objective_matrix_abs_max()
-                 << " and smallest non-zero absolute value "
-                 << problem_stats.objective_matrix_abs_min()
-                 << "; performance may suffer.";
+    SOLVER_LOG(&logger, "WARNING: Objective matrix has largest absolute value ",
+               problem_stats.objective_matrix_abs_max(),
+               " and smallest non-zero absolute value ",
+               problem_stats.objective_matrix_abs_min(),
+               "; performance may suffer.");
   }
   return std::nullopt;
 }
 
 std::optional<SolverResult> CheckInitialSolution(
     const ShardedQuadraticProgram& sharded_qp,
-    const PrimalAndDualSolution& initial_solution) {
+    const PrimalAndDualSolution& initial_solution, SolverLogger& logger) {
   const double kExcessiveInputValue = 1e50;
   if (initial_solution.primal_solution.size() != sharded_qp.PrimalSize()) {
     return ErrorSolverResult(
@@ -917,12 +940,13 @@ std::optional<SolverResult> CheckInitialSolution(
         absl::StrCat("Initial primal solution has size ",
                      initial_solution.primal_solution.size(),
                      " which differs from problem primal size ",
-                     sharded_qp.PrimalSize()));
+                     sharded_qp.PrimalSize()),
+        logger);
   }
   if (std::isnan(
           Norm(initial_solution.primal_solution, sharded_qp.PrimalSharder()))) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_INITIAL_SOLUTION,
-                             "Initial primal solution has a NAN.");
+                             "Initial primal solution has a NAN.", logger);
   }
   if (const double norm = LInfNorm(initial_solution.primal_solution,
                                    sharded_qp.PrimalSharder());
@@ -931,7 +955,8 @@ std::optional<SolverResult> CheckInitialSolution(
         TERMINATION_REASON_INVALID_INITIAL_SOLUTION,
         absl::StrCat(
             "Initial primal solution has an entry with absolute value ", norm,
-            " which exceeds limit of ", kExcessiveInputValue));
+            " which exceeds limit of ", kExcessiveInputValue),
+        logger);
   }
   if (initial_solution.dual_solution.size() != sharded_qp.DualSize()) {
     return ErrorSolverResult(
@@ -939,12 +964,13 @@ std::optional<SolverResult> CheckInitialSolution(
         absl::StrCat("Initial dual solution has size ",
                      initial_solution.dual_solution.size(),
                      " which differs from problem dual size ",
-                     sharded_qp.DualSize()));
+                     sharded_qp.DualSize()),
+        logger);
   }
   if (std::isnan(
           Norm(initial_solution.dual_solution, sharded_qp.DualSharder()))) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_INITIAL_SOLUTION,
-                             "Initial dual solution has a NAN.");
+                             "Initial dual solution has a NAN.", logger);
   }
   if (const double norm =
           LInfNorm(initial_solution.dual_solution, sharded_qp.DualSharder());
@@ -952,7 +978,8 @@ std::optional<SolverResult> CheckInitialSolution(
     return ErrorSolverResult(
         TERMINATION_REASON_INVALID_INITIAL_SOLUTION,
         absl::StrCat("Initial dual solution has an entry with absolute value ",
-                     norm, " which exceeds limit of ", kExcessiveInputValue));
+                     norm, " which exceeds limit of ", kExcessiveInputValue),
+        logger);
   }
   return std::nullopt;
 }
@@ -977,7 +1004,8 @@ SolverResult PreprocessSolver::PreprocessAndSolve(
         "The input problem has invalid bounds (after replacing large "
         "constraint bounds with infinity): some variable or constraint has "
         "lower_bound > upper_bound, lower_bound == inf, or upper_bound == "
-        "-inf.");
+        "-inf.",
+        logger_);
   }
   if (Qp().objective_matrix.has_value() &&
       !sharded_qp_.PrimalSharder().ParallelTrueForAllShards(
@@ -987,20 +1015,21 @@ SolverResult PreprocessSolver::PreprocessAndSolve(
           })) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
                              "The objective is not convex (i.e., the objective "
-                             "matrix contains negative or NAN entries).");
+                             "matrix contains negative or NAN entries).",
+                             logger_);
   }
   *solve_log.mutable_original_problem_stats() = ComputeStats(sharded_qp_);
   const QuadraticProgramStats& original_problem_stats =
       solve_log.original_problem_stats();
   if (auto maybe_result =
           CheckProblemStats(original_problem_stats, Qp().objective_offset,
-                            params.presolve_options().use_glop());
+                            params.presolve_options().use_glop(), logger_);
       maybe_result.has_value()) {
     return *maybe_result;
   }
   if (initial_solution.has_value()) {
     if (auto maybe_result =
-            CheckInitialSolution(sharded_qp_, *initial_solution);
+            CheckInitialSolution(sharded_qp_, *initial_solution, logger_);
         maybe_result.has_value()) {
       return *maybe_result;
     }
@@ -1010,8 +1039,7 @@ SolverResult PreprocessSolver::PreprocessAndSolve(
       params.presolve_options().use_glop() ? "presolving and " : "",
       "rescaling:");
   if (params.verbosity_level() >= 1) {
-    LogInfoWithoutPrefix(
-        absl::StrCat("Problem stats before ", preprocessing_string));
+    SOLVER_LOG(&logger_, "Problem stats before ", preprocessing_string);
     LogQuadraticProgramStats(solve_log.original_problem_stats());
   }
   iteration_stats_callback_ = std::move(iteration_stats_callback);
@@ -1049,17 +1077,21 @@ SolverResult PreprocessSolver::PreprocessAndSolve(
     } else {
       if (*maybe_terminate == TERMINATION_REASON_OPTIMAL) {
         final_termination_reason = TERMINATION_REASON_NUMERICAL_ERROR;
-        LOG(WARNING) << "Presolve claimed to solve the LP optimally but the "
-                        "solution doesn't satisfy the optimality criteria.";
+        SOLVER_LOG(
+            &logger_,
+            "WARNING: Presolve claimed to solve the LP optimally but the "
+            "solution doesn't satisfy the optimality criteria.");
       } else {
         final_termination_reason = *maybe_terminate;
       }
     }
     return ConstructOriginalSolverResult(
-        params, ConstructSolverResult(
-                    std::move(working_primal), std::move(working_dual),
-                    std::move(iteration_stats), final_termination_reason,
-                    POINT_TYPE_PRESOLVER_SOLUTION, std::move(solve_log)));
+        params,
+        ConstructSolverResult(
+            std::move(working_primal), std::move(working_dual),
+            std::move(iteration_stats), final_termination_reason,
+            POINT_TYPE_PRESOLVER_SOLUTION, std::move(solve_log)),
+        logger_);
   }
 
   VectorXd starting_primal_solution;
@@ -1081,8 +1113,7 @@ SolverResult PreprocessSolver::PreprocessAndSolve(
                            starting_dual_solution);
   *solve_log.mutable_preprocessed_problem_stats() = ComputeStats(sharded_qp_);
   if (params.verbosity_level() >= 1) {
-    LogInfoWithoutPrefix(
-        absl::StrCat("Problem stats after ", preprocessing_string));
+    SOLVER_LOG(&logger_, "Problem stats after ", preprocessing_string);
     LogQuadraticProgramStats(solve_log.preprocessed_problem_stats());
   }
 
@@ -1131,7 +1162,7 @@ SolverResult PreprocessSolver::PreprocessAndSolve(
   solve_log.set_preprocessing_time_sec(timer.Get());
   SolverResult result = solver.Solve(IterationType::kNormal, interrupt_solve,
                                      std::move(solve_log));
-  return ConstructOriginalSolverResult(params, std::move(result));
+  return ConstructOriginalSolverResult(params, std::move(result), logger_);
 }
 
 glop::GlopParameters PreprocessSolver::PreprocessorParameters(
@@ -1151,7 +1182,7 @@ glop::GlopParameters PreprocessSolver::PreprocessorParameters(
 }
 
 TerminationReason GlopStatusToTerminationReason(
-    const glop::ProblemStatus glop_status) {
+    const glop::ProblemStatus glop_status, SolverLogger& logger) {
   switch (glop_status) {
     case glop::ProblemStatus::OPTIMAL:
       return TERMINATION_REASON_OPTIMAL;
@@ -1167,7 +1198,8 @@ TerminationReason GlopStatusToTerminationReason(
     case glop::ProblemStatus::PRIMAL_UNBOUNDED:
       return TERMINATION_REASON_PRIMAL_OR_DUAL_INFEASIBLE;
     default:
-      LOG(WARNING) << "Unexpected preprocessor status " << glop_status;
+      SOLVER_LOG(&logger, "WARNING: Unexpected preprocessor status ",
+                 glop_status);
       return TERMINATION_REASON_OTHER;
   }
 }
@@ -1180,20 +1212,23 @@ std::optional<TerminationReason> PreprocessSolver::ApplyPresolveIfEnabled(
     return std::nullopt;
   }
   if (!IsLinearProgram(Qp())) {
-    LOG(WARNING)
-        << "Skipping presolve, which is only supported for linear programs";
+    SOLVER_LOG(&logger_,
+               "WARNING: Skipping presolve, which is only supported for linear "
+               "programs");
     return std::nullopt;
   }
   absl::StatusOr<MPModelProto> model = QpToMpModelProto(Qp());
   if (!model.ok()) {
-    LOG(WARNING)
-        << "Skipping presolve because of error converting to MPModelProto: "
-        << model.status();
+    SOLVER_LOG(&logger_,
+               "WARNING: Skipping presolve because of error converting to "
+               "MPModelProto: ",
+               model.status());
     return std::nullopt;
   }
   if (initial_solution->has_value()) {
-    LOG(WARNING) << "Ignoring initial solution. Initial solutions "
-                    "are ignored when presolve is on.";
+    SOLVER_LOG(&logger_,
+               "WARNING: Ignoring initial solution. Initial solutions are "
+               "ignored when presolve is on.");
     initial_solution->reset();
   }
   glop::LinearProgram glop_lp;
@@ -1226,7 +1261,8 @@ std::optional<TerminationReason> PreprocessSolver::ApplyPresolveIfEnabled(
   if (presolve_info_->preprocessor.status() != glop::ProblemStatus::INIT) {
     col_scaling_vec_ = OnesVector(sharded_qp_.PrimalSharder());
     row_scaling_vec_ = OnesVector(sharded_qp_.DualSharder());
-    return GlopStatusToTerminationReason(presolve_info_->preprocessor.status());
+    return GlopStatusToTerminationReason(presolve_info_->preprocessor.status(),
+                                         logger_);
   }
   return std::nullopt;
 }
@@ -1248,51 +1284,63 @@ void PreprocessSolver::ComputeAndApplyRescaling(
 }
 
 void PreprocessSolver::LogQuadraticProgramStats(
-    const QuadraticProgramStats& stats) {
-  LogInfoWithoutPrefix(
-      absl::StrFormat("There are %i variables, %i constraints, and %i "
-                      "constraint matrix nonzeros.",
-                      stats.num_variables(), stats.num_constraints(),
-                      stats.constraint_matrix_num_nonzeros()));
+    const QuadraticProgramStats& stats) const {
+  SOLVER_LOG(&logger_,
+             absl::StrFormat("There are %i variables, %i constraints, and %i "
+                             "constraint matrix nonzeros.",
+                             stats.num_variables(), stats.num_constraints(),
+                             stats.constraint_matrix_num_nonzeros()));
   if (Qp().constraint_matrix.nonZeros() > 0) {
-    LogInfoWithoutPrefix(absl::StrFormat(
-        "Absolute values of nonzero constraint matrix elements: largest=%f, "
-        "smallest=%f, avg=%f",
-        stats.constraint_matrix_abs_max(), stats.constraint_matrix_abs_min(),
-        stats.constraint_matrix_abs_avg()));
-    LogInfoWithoutPrefix(
+    SOLVER_LOG(&logger_,
+               absl::StrFormat("Absolute values of nonzero constraint matrix "
+                               "elements: largest=%f, "
+                               "smallest=%f, avg=%f",
+                               stats.constraint_matrix_abs_max(),
+                               stats.constraint_matrix_abs_min(),
+                               stats.constraint_matrix_abs_avg()));
+    SOLVER_LOG(
+        &logger_,
         absl::StrFormat("Constraint matrix, infinity norm: max(row & col)=%f, "
                         "min_col=%f, min_row=%f",
                         stats.constraint_matrix_abs_max(),
                         stats.constraint_matrix_col_min_l_inf_norm(),
                         stats.constraint_matrix_row_min_l_inf_norm()));
-    LogInfoWithoutPrefix(absl::StrFormat(
-        "Constraint bounds statistics (max absolute value per row): "
-        "largest=%f, smallest=%f, avg=%f, l2_norm=%f",
-        stats.combined_bounds_max(), stats.combined_bounds_min(),
-        stats.combined_bounds_avg(), stats.combined_bounds_l2_norm()));
+    SOLVER_LOG(
+        &logger_,
+        absl::StrFormat(
+            "Constraint bounds statistics (max absolute value per row): "
+            "largest=%f, smallest=%f, avg=%f, l2_norm=%f",
+            stats.combined_bounds_max(), stats.combined_bounds_min(),
+            stats.combined_bounds_avg(), stats.combined_bounds_l2_norm()));
   }
   if (!IsLinearProgram(Qp())) {
-    LogInfoWithoutPrefix(absl::StrFormat(
-        "There are %i nonzero diagonal coefficients in the objective matrix.",
-        stats.objective_matrix_num_nonzeros()));
-    LogInfoWithoutPrefix(absl::StrFormat(
-        "Absolute values of nonzero objective matrix elements: largest=%f, "
-        "smallest=%f, avg=%f",
-        stats.objective_matrix_abs_max(), stats.objective_matrix_abs_min(),
-        stats.objective_matrix_abs_avg()));
+    SOLVER_LOG(&logger_,
+               absl::StrFormat("There are %i nonzero diagonal coefficients in "
+                               "the objective matrix.",
+                               stats.objective_matrix_num_nonzeros()));
+    SOLVER_LOG(
+        &logger_,
+        absl::StrFormat(
+            "Absolute values of nonzero objective matrix elements: largest=%f, "
+            "smallest=%f, avg=%f",
+            stats.objective_matrix_abs_max(), stats.objective_matrix_abs_min(),
+            stats.objective_matrix_abs_avg()));
   }
-  LogInfoWithoutPrefix(absl::StrFormat(
-      "Absolute values of objective vector elements: largest=%f, smallest=%f, "
-      "avg=%f, l2_norm=%f",
-      stats.objective_vector_abs_max(), stats.objective_vector_abs_min(),
-      stats.objective_vector_abs_avg(), stats.objective_vector_l2_norm()));
-  LogInfoWithoutPrefix(absl::StrFormat(
-      "Gaps between variable upper and lower bounds: #finite=%i of %i, "
-      "largest=%f, smallest=%f, avg=%f",
-      stats.variable_bound_gaps_num_finite(), stats.num_variables(),
-      stats.variable_bound_gaps_max(), stats.variable_bound_gaps_min(),
-      stats.variable_bound_gaps_avg()));
+  SOLVER_LOG(&logger_, absl::StrFormat("Absolute values of objective vector "
+                                       "elements: largest=%f, smallest=%f, "
+                                       "avg=%f, l2_norm=%f",
+                                       stats.objective_vector_abs_max(),
+                                       stats.objective_vector_abs_min(),
+                                       stats.objective_vector_abs_avg(),
+                                       stats.objective_vector_l2_norm()));
+  SOLVER_LOG(
+      &logger_,
+      absl::StrFormat(
+          "Gaps between variable upper and lower bounds: #finite=%i of %i, "
+          "largest=%f, smallest=%f, avg=%f",
+          stats.variable_bound_gaps_num_finite(), stats.num_variables(),
+          stats.variable_bound_gaps_max(), stats.variable_bound_gaps_min(),
+          stats.variable_bound_gaps_avg()));
 }
 
 double PreprocessSolver::InitialPrimalWeight(
@@ -1508,22 +1556,22 @@ PreprocessSolver::UpdateIterationStatsAndCheckTermination(
            absl::Seconds(params.log_interval_seconds()))) {
     if (log_counter_ == 0) {
       LogIterationStatsHeader(params.verbosity_level(),
-                              params.use_feasibility_polishing());
+                              params.use_feasibility_polishing(), logger_);
     }
     LogIterationStats(params.verbosity_level(),
                       params.use_feasibility_polishing(), iteration_type, stats,
                       params.termination_criteria(), original_bound_norms_,
-                      POINT_TYPE_AVERAGE_ITERATE);
+                      POINT_TYPE_AVERAGE_ITERATE, logger_);
     if (params.verbosity_level() >= 4) {
       // If the convergence information doesn't contain an average iterate, the
       // previous `LogIterationStats()` will report some other iterate (usually
       // the current one) so don't repeat logging it.
       if (GetConvergenceInformation(stats, POINT_TYPE_AVERAGE_ITERATE) !=
           std::nullopt) {
-        LogIterationStats(params.verbosity_level(),
-                          params.use_feasibility_polishing(), iteration_type,
-                          stats, params.termination_criteria(),
-                          original_bound_norms_, POINT_TYPE_CURRENT_ITERATE);
+        LogIterationStats(
+            params.verbosity_level(), params.use_feasibility_polishing(),
+            iteration_type, stats, params.termination_criteria(),
+            original_bound_norms_, POINT_TYPE_CURRENT_ITERATE, logger_);
       }
     }
     time_of_last_log_ = logging_time;
@@ -1577,11 +1625,19 @@ void PreprocessSolver::ComputeConvergenceAndInfeasibilityFromWorkingSolution(
           candidate_type);
     }
     if (infeasibility_information != nullptr) {
+      VectorXd primal_copy =
+          CloneVector(original.primal_solution,
+                      presolve_info_->sharded_original_qp.PrimalSharder());
+      ProjectToPrimalVariableBounds(presolve_info_->sharded_original_qp,
+                                    primal_copy,
+                                    /*use_feasibility_bounds=*/true);
+      // Only iterate differences should be able to violate the dual variable
+      // bounds, and iterate differences aren't used when presolve is enabled.
       *infeasibility_information = ComputeInfeasibilityInformation(
           params, presolve_info_->sharded_original_qp,
           presolve_info_->trivial_col_scaling_vec,
-          presolve_info_->trivial_row_scaling_vec, original.primal_solution,
-          original.dual_solution, candidate_type);
+          presolve_info_->trivial_row_scaling_vec, primal_copy,
+          original.dual_solution, original.primal_solution, candidate_type);
     }
   } else {
     if (convergence_information != nullptr) {
@@ -1591,9 +1647,23 @@ void PreprocessSolver::ComputeConvergenceAndInfeasibilityFromWorkingSolution(
           dual_epsilon_ratio, candidate_type);
     }
     if (infeasibility_information != nullptr) {
-      *infeasibility_information = ComputeInfeasibilityInformation(
-          params, sharded_qp_, col_scaling_vec_, row_scaling_vec_,
-          working_primal, working_dual, candidate_type);
+      VectorXd primal_copy =
+          CloneVector(working_primal, sharded_qp_.PrimalSharder());
+      ProjectToPrimalVariableBounds(sharded_qp_, primal_copy,
+                                    /*use_feasibility_bounds=*/true);
+      if (candidate_type == POINT_TYPE_ITERATE_DIFFERENCE) {
+        // Iterate differences might not satisfy the dual variable bounds.
+        VectorXd dual_copy =
+            CloneVector(working_dual, sharded_qp_.DualSharder());
+        ProjectToDualVariableBounds(sharded_qp_, dual_copy);
+        *infeasibility_information = ComputeInfeasibilityInformation(
+            params, sharded_qp_, col_scaling_vec_, row_scaling_vec_,
+            primal_copy, dual_copy, working_primal, candidate_type);
+      } else {
+        *infeasibility_information = ComputeInfeasibilityInformation(
+            params, sharded_qp_, col_scaling_vec_, row_scaling_vec_,
+            primal_copy, working_dual, working_primal, candidate_type);
+      }
     }
   }
 }
@@ -1601,7 +1671,8 @@ void PreprocessSolver::ComputeConvergenceAndInfeasibilityFromWorkingSolution(
 // `result` is used both as the input and as the temporary that will be
 // returned.
 SolverResult PreprocessSolver::ConstructOriginalSolverResult(
-    const PrimalDualHybridGradientParams& params, SolverResult result) const {
+    const PrimalDualHybridGradientParams& params, SolverResult result,
+    SolverLogger& logger) const {
   const bool use_zero_primal_objective =
       result.solve_log.termination_reason() ==
       TERMINATION_REASON_PRIMAL_INFEASIBLE;
@@ -1611,6 +1682,14 @@ SolverResult PreprocessSolver::ConstructOriginalSolverResult(
         {.primal_solution = std::move(result.primal_solution),
          .dual_solution = std::move(result.dual_solution)});
     result.primal_solution = std::move(original_solution.primal_solution);
+    if (result.solve_log.termination_reason() ==
+        TERMINATION_REASON_DUAL_INFEASIBLE) {
+      ProjectToPrimalVariableBounds(presolve_info_->sharded_original_qp,
+                                    result.primal_solution,
+                                    /*use_feasibility_bounds=*/true);
+    }
+    // If presolve is enabled, no projection of the dual is performed when
+    // checking for infeasibility.
     result.dual_solution = std::move(original_solution.dual_solution);
     // `RecoverOriginalSolution` doesn't recover reduced costs so we need to
     // compute them with respect to the original problem.
@@ -1618,6 +1697,15 @@ SolverResult PreprocessSolver::ConstructOriginalSolverResult(
         params, presolve_info_->sharded_original_qp, result.primal_solution,
         result.dual_solution, use_zero_primal_objective);
   } else {
+    if (result.solve_log.termination_reason() ==
+        TERMINATION_REASON_DUAL_INFEASIBLE) {
+      ProjectToPrimalVariableBounds(sharded_qp_, result.primal_solution,
+                                    /*use_feasibility_bounds=*/true);
+    }
+    if (result.solve_log.termination_reason() ==
+        TERMINATION_REASON_PRIMAL_INFEASIBLE) {
+      ProjectToDualVariableBounds(sharded_qp_, result.dual_solution);
+    }
     result.reduced_costs =
         ReducedCosts(params, sharded_qp_, result.primal_solution,
                      result.dual_solution, use_zero_primal_objective);
@@ -1650,27 +1738,24 @@ SolverResult PreprocessSolver::ConstructOriginalSolverResult(
   }
 
   if (params.verbosity_level() >= 1) {
-    LogInfoWithoutPrefix(absl::StrCat(
-        "Termination reason: ",
-        TerminationReason_Name(result.solve_log.termination_reason())));
-    LogInfoWithoutPrefix(
-        absl::StrCat("Solution point type: ",
-                     PointType_Name(result.solve_log.solution_type())));
-    LogInfoWithoutPrefix("Final solution stats:");
+    SOLVER_LOG(&logger, "Termination reason: ",
+               TerminationReason_Name(result.solve_log.termination_reason()));
+    SOLVER_LOG(&logger, "Solution point type: ",
+               PointType_Name(result.solve_log.solution_type()));
+    SOLVER_LOG(&logger, "Final solution stats:");
     LogIterationStatsHeader(params.verbosity_level(),
-                            params.use_feasibility_polishing());
+                            params.use_feasibility_polishing(), logger);
     LogIterationStats(params.verbosity_level(),
                       params.use_feasibility_polishing(), iteration_type,
                       result.solve_log.solution_stats(),
                       params.termination_criteria(), original_bound_norms_,
-                      result.solve_log.solution_type());
+                      result.solve_log.solution_type(), logger);
     const auto& convergence_info = GetConvergenceInformation(
         result.solve_log.solution_stats(), result.solve_log.solution_type());
     if (convergence_info.has_value()) {
       if (std::isfinite(convergence_info->corrected_dual_objective())) {
-        LogInfoWithoutPrefix(
-            absl::StrCat("Dual objective after infeasibility correction: ",
-                         convergence_info->corrected_dual_objective()));
+        SOLVER_LOG(&logger, "Dual objective after infeasibility correction: ",
+                   convergence_info->corrected_dual_objective());
       }
     }
   }
@@ -2037,9 +2122,8 @@ double Solver::ComputeNewPrimalWeight() const {
       std::exp(smoothing_param * std::log(unsmoothed_new_primal_weight) +
                (1.0 - smoothing_param) * std::log(primal_weight_));
   if (params_.verbosity_level() >= 4) {
-    LogInfoWithoutPrefix(absl::StrCat("New computed primal weight is ",
-                                      new_primal_weight, " at iteration ",
-                                      iterations_completed_));
+    SOLVER_LOG(&preprocess_solver_->Logger(), "New computed primal weight is ",
+               new_primal_weight, " at iteration ", iterations_completed_);
   }
   return new_primal_weight;
 }
@@ -2081,16 +2165,16 @@ void Solver::ApplyRestartChoice(const RestartChoice restart_to_apply) {
       return;
     case RESTART_CHOICE_WEIGHTED_AVERAGE_RESET:
       if (params_.verbosity_level() >= 4) {
-        LogInfoWithoutPrefix(absl::StrCat(
-            "Restarted to current on iteration ", iterations_completed_,
-            " after ", primal_average_.NumTerms(), " iterations"));
+        SOLVER_LOG(&preprocess_solver_->Logger(),
+                   "Restarted to current on iteration ", iterations_completed_,
+                   " after ", primal_average_.NumTerms(), " iterations");
       }
       break;
     case RESTART_CHOICE_RESTART_TO_AVERAGE:
       if (params_.verbosity_level() >= 4) {
-        LogInfoWithoutPrefix(absl::StrCat(
-            "Restarted to average on iteration ", iterations_completed_,
-            " after ", primal_average_.NumTerms(), " iterations"));
+        SOLVER_LOG(&preprocess_solver_->Logger(),
+                   "Restarted to average on iteration ", iterations_completed_,
+                   " after ", primal_average_.NumTerms(), " iterations");
       }
       current_primal_solution_ = primal_average_.ComputeAverage();
       current_dual_solution_ = dual_average_.ComputeAverage();
@@ -2233,14 +2317,16 @@ void Solver::ResetAverageToCurrent() {
 
 void Solver::LogNumericalTermination() const {
   if (params_.verbosity_level() >= 2) {
-    LogInfoWithoutPrefix(absl::StrCat(
-        "Forced numerical termination at iteration ", iterations_completed_));
+    SOLVER_LOG(&preprocess_solver_->Logger(),
+               "Forced numerical termination at iteration ",
+               iterations_completed_);
   }
 }
 
 void Solver::LogInnerIterationLimitHit() const {
-  LOG(WARNING) << "Inner iteration limit reached at iteration "
-               << iterations_completed_;
+  SOLVER_LOG(&preprocess_solver_->Logger(),
+             "WARNING: Inner iteration limit reached at iteration ",
+             iterations_completed_);
 }
 
 InnerStepOutcome Solver::TakeMalitskyPockStep() {
@@ -2490,23 +2576,25 @@ std::optional<SolverResult> Solver::TryFeasibilityPolishing(
 
   if (!ObjectiveGapMet(optimality_criteria, first_convergence_info)) {
     if (params_.verbosity_level() >= 2) {
-      LogInfoWithoutPrefix(
-          "Skipping feasibility polishing because the objective gap "
-          "is too large.");
+      SOLVER_LOG(&preprocess_solver_->Logger(),
+                 "Skipping feasibility polishing because the objective gap "
+                 "is too large.");
     }
     return std::nullopt;
   }
 
   if (params_.verbosity_level() >= 2) {
-    LogInfoWithoutPrefix("Starting primal feasibility polishing");
+    SOLVER_LOG(&preprocess_solver_->Logger(),
+               "Starting primal feasibility polishing");
   }
   SolverResult primal_result = TryPrimalPolishing(
       std::move(average_primal), iteration_limit, interrupt_solve, solve_log);
 
   if (params_.verbosity_level() >= 2) {
-    LogInfoWithoutPrefix(absl::StrCat(
+    SOLVER_LOG(
+        &preprocess_solver_->Logger(),
         "Primal feasibility polishing termination reason: ",
-        TerminationReason_Name(primal_result.solve_log.termination_reason())));
+        TerminationReason_Name(primal_result.solve_log.termination_reason()));
   }
   if (TerminationReasonIsWorkLimit(
           primal_result.solve_log.termination_reason())) {
@@ -2518,21 +2606,24 @@ std::optional<SolverResult> Solver::TryFeasibilityPolishing(
     // we ignore the polishing result indicating infeasibility.
     // `TERMINATION_REASON_NUMERICAL_ERROR` can occur, but would be surprising
     // and interesting. Other termination reasons are probably bugs.
-    LOG(WARNING) << "Primal feasibility polishing terminated with error "
-                 << primal_result.solve_log.termination_reason();
+    SOLVER_LOG(&preprocess_solver_->Logger(),
+               "WARNING: Primal feasibility polishing terminated with error ",
+               primal_result.solve_log.termination_reason());
     return std::nullopt;
   }
 
   if (params_.verbosity_level() >= 2) {
-    LogInfoWithoutPrefix("Starting dual feasibility polishing");
+    SOLVER_LOG(&preprocess_solver_->Logger(),
+               "Starting dual feasibility polishing");
   }
   SolverResult dual_result = TryDualPolishing(
       std::move(average_dual), iteration_limit, interrupt_solve, solve_log);
 
   if (params_.verbosity_level() >= 2) {
-    LogInfoWithoutPrefix(absl::StrCat(
+    SOLVER_LOG(
+        &preprocess_solver_->Logger(),
         "Dual feasibility polishing termination reason: ",
-        TerminationReason_Name(dual_result.solve_log.termination_reason())));
+        TerminationReason_Name(dual_result.solve_log.termination_reason()));
   }
 
   if (TerminationReasonIsWorkLimit(
@@ -2543,8 +2634,9 @@ std::optional<SolverResult> Solver::TryFeasibilityPolishing(
     // Note: The comment in the corresponding location when checking the
     // termination reason for primal feasibility polishing applies here
     // too.
-    LOG(WARNING) << "Dual feasibility polishing terminated with error "
-                 << dual_result.solve_log.termination_reason();
+    SOLVER_LOG(&preprocess_solver_->Logger(),
+               "WARNING: Dual feasibility polishing terminated with error ",
+               dual_result.solve_log.termination_reason());
     return std::nullopt;
   }
 
@@ -2554,15 +2646,18 @@ std::optional<SolverResult> Solver::TryFeasibilityPolishing(
       POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION,
       full_stats.add_convergence_information(), nullptr);
   if (params_.verbosity_level() >= 2) {
-    LogInfoWithoutPrefix("solution stats for polished solution:");
+    SOLVER_LOG(&preprocess_solver_->Logger(),
+               "solution stats for polished solution:");
     LogIterationStatsHeader(params_.verbosity_level(),
-                            /*use_feasibility_polishing=*/true);
+                            /*use_feasibility_polishing=*/true,
+                            preprocess_solver_->Logger());
     LogIterationStats(params_.verbosity_level(),
                       /*use_feasibility_polishing=*/true,
                       IterationType::kFeasibilityPolishingTermination,
                       full_stats, params_.termination_criteria(),
                       preprocess_solver_->OriginalBoundNorms(),
-                      POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION);
+                      POINT_TYPE_FEASIBILITY_POLISHING_SOLUTION,
+                      preprocess_solver_->Logger());
   }
   std::optional<TerminationReasonAndPointType> earned_termination =
       CheckIterateTerminationCriteria(params_.termination_criteria(),
@@ -2813,9 +2908,10 @@ SolverResult Solver::Solve(const IterationType iteration_type,
 SolverResult PrimalDualHybridGradient(
     QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
     const std::atomic<bool>* interrupt_solve,
+    std::function<void(const std::string&)> message_callback,
     IterationStatsCallback iteration_stats_callback) {
   return PrimalDualHybridGradient(std::move(qp), params, std::nullopt,
-                                  interrupt_solve,
+                                  interrupt_solve, std::move(message_callback),
                                   std::move(iteration_stats_callback));
 }
 
@@ -2823,33 +2919,44 @@ SolverResult PrimalDualHybridGradient(
     QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
     std::optional<PrimalAndDualSolution> initial_solution,
     const std::atomic<bool>* interrupt_solve,
+    std::function<void(const std::string&)> message_callback,
     IterationStatsCallback iteration_stats_callback) {
+  SolverLogger logger;
+  logger.EnableLogging(true);
+  if (message_callback) {
+    logger.AddInfoLoggingCallback(std::move(message_callback));
+  } else {
+    logger.SetLogToStdOut(true);
+  }
   const absl::Status params_status =
       ValidatePrimalDualHybridGradientParams(params);
   if (!params_status.ok()) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PARAMETER,
-                             params_status.ToString());
+                             params_status.ToString(), logger);
   }
   if (!qp.constraint_matrix.isCompressed()) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
                              "constraint_matrix must be in compressed format. "
-                             "Call constraint_matrix.makeCompressed()");
+                             "Call constraint_matrix.makeCompressed()",
+                             logger);
   }
   const absl::Status dimensions_status = ValidateQuadraticProgramDimensions(qp);
   if (!dimensions_status.ok()) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
-                             dimensions_status.ToString());
+                             dimensions_status.ToString(), logger);
   }
   if (qp.objective_scaling_factor == 0) {
     return ErrorSolverResult(TERMINATION_REASON_INVALID_PROBLEM,
-                             "The objective scaling factor cannot be zero.");
+                             "The objective scaling factor cannot be zero.",
+                             logger);
   }
   if (params.use_feasibility_polishing() && !IsLinearProgram(qp)) {
     return ErrorSolverResult(
         TERMINATION_REASON_INVALID_PARAMETER,
-        "use_feasibility_polishing is only implemented for linear programs.");
+        "use_feasibility_polishing is only implemented for linear programs.",
+        logger);
   }
-  PreprocessSolver solver(std::move(qp), params);
+  PreprocessSolver solver(std::move(qp), params, &logger);
   return solver.PreprocessAndSolve(params, std::move(initial_solution),
                                    interrupt_solve,
                                    std::move(iteration_stats_callback));
diff --git a/ortools/pdlp/primal_dual_hybrid_gradient.h b/ortools/pdlp/primal_dual_hybrid_gradient.h
index 12199d9365..70d8c3628d 100644
--- a/ortools/pdlp/primal_dual_hybrid_gradient.h
+++ b/ortools/pdlp/primal_dual_hybrid_gradient.h
@@ -17,6 +17,7 @@
 #include <atomic>
 #include <functional>
 #include <optional>
+#include <string>
 
 #include "Eigen/Core"
 #include "ortools/lp_data/lp_data.h"
@@ -128,6 +129,12 @@ struct IterationCallbackInfo {
 // if `interrupt_solve->load()` is true, in which case the solve will terminate
 // with `TERMINATION_REASON_INTERRUPTED_BY_USER`.
 //
+// If `message_callback` is not nullptr, solver logging will be passed to
+// `message_callback`. Otherwise solver logging will be written to stdout.
+// In either case, the amount of logging is controlled by `verbosity_level`.
+// In particular, if `verbosity_level == 0`, there will be no logging in either
+// case.
+//
 // If `iteration_stats_callback` is not nullptr, then at each termination step
 // (when iteration stats are logged), `iteration_stats_callback` will also
 // be called with those iteration stats.
@@ -144,6 +151,7 @@ struct IterationCallbackInfo {
 SolverResult PrimalDualHybridGradient(
     QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
     const std::atomic<bool>* interrupt_solve = nullptr,
+    std::function<void(const std::string&)> message_callback = nullptr,
     std::function<void(const IterationCallbackInfo&)> iteration_stats_callback =
         nullptr);
 
@@ -156,6 +164,7 @@ SolverResult PrimalDualHybridGradient(
     QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
     std::optional<PrimalAndDualSolution> initial_solution,
     const std::atomic<bool>* interrupt_solve = nullptr,
+    std::function<void(const std::string&)> message_callback = nullptr,
     std::function<void(const IterationCallbackInfo&)> iteration_stats_callback =
         nullptr);
 
diff --git a/ortools/pdlp/primal_dual_hybrid_gradient_test.cc b/ortools/pdlp/primal_dual_hybrid_gradient_test.cc
index abe96fc321..952633f583 100644
--- a/ortools/pdlp/primal_dual_hybrid_gradient_test.cc
+++ b/ortools/pdlp/primal_dual_hybrid_gradient_test.cc
@@ -47,6 +47,7 @@
 #include "ortools/pdlp/solvers.pb.h"
 #include "ortools/pdlp/termination.h"
 #include "ortools/pdlp/test_util.h"
+#include "testing/base/public/googletest.h"
 
 namespace operations_research::pdlp {
 namespace {
@@ -443,6 +444,36 @@ TEST_P(PrimalDualHybridGradientLPTest, InfeasibleDual) {
   EXPECT_LE(output.solve_log.iteration_count(), iteration_upperbound);
 }
 
+TEST_P(PrimalDualHybridGradientLPTest, InfeasiblePrimalWithReducedCosts) {
+  // A trivial LP that has reduced costs within the dual ray.
+  //     min x
+  //     Constraint: 2 <= x
+  //     Variable: 0 <= x <= 1
+  QuadraticProgram lp(1, 1);
+  lp.objective_vector = Eigen::VectorXd{{1}};
+  lp.constraint_lower_bounds = Eigen::VectorXd{{2}};
+  lp.constraint_upper_bounds =
+      Eigen::VectorXd{{std::numeric_limits<double>::infinity()}};
+  lp.variable_lower_bounds = Eigen::VectorXd{{0}};
+  lp.variable_upper_bounds = Eigen::VectorXd{{1}};
+  lp.constraint_matrix.coeffRef(0, 0) = 1.0;
+  lp.constraint_matrix.makeCompressed();
+  const int iteration_upperbound = 100;
+  PrimalDualHybridGradientParams params =
+      CreateSolverParamsForFixture(iteration_upperbound,
+                                   /*eps_optimal_absolute=*/1.0e-6);
+  params.set_major_iteration_frequency(5);
+  params.mutable_termination_criteria()->set_eps_primal_infeasible(1.0e-6);
+
+  SolverResult output = PrimalDualHybridGradient(lp, params);
+  EXPECT_EQ(output.solve_log.termination_reason(),
+            TERMINATION_REASON_PRIMAL_INFEASIBLE);
+  const auto& dual = output.dual_solution;
+  EXPECT_GT(dual[0], 0.0);
+  EXPECT_EQ(output.reduced_costs[0], -dual[0]);
+  EXPECT_LT(output.solve_log.iteration_count(), iteration_upperbound);
+}
+
 TEST_P(PrimalDualHybridGradientLPTest, InfeasiblePrimalDual) {
   const int iteration_upperbound = 600;
   PrimalDualHybridGradientParams params =
@@ -1425,6 +1456,7 @@ TEST(PrimalDualHybridGradientTest, CallsCallback) {
   int callback_count = 0;
   SolverResult output = PrimalDualHybridGradient(
       TestLp(), params, /*interrupt_solve=*/nullptr,
+      /*message_callback=*/nullptr,
       [&callback_count](const IterationCallbackInfo& callback_info) {
         ++callback_count;
       });
@@ -1525,7 +1557,8 @@ TEST(PrimalDualHybridGradientTest, RespectsInterruptFromCallback) {
   };
 
   const SolverResult output =
-      PrimalDualHybridGradient(TestLp(), params, &interrupt_solve, callback);
+      PrimalDualHybridGradient(TestLp(), params, &interrupt_solve,
+                               /*message_callback=*/nullptr, callback);
   EXPECT_EQ(output.solve_log.termination_reason(),
             TERMINATION_REASON_INTERRUPTED_BY_USER);
   EXPECT_GE(output.solve_log.iteration_count(), 10);
@@ -1921,6 +1954,7 @@ TEST_F(FeasibilityPolishingPrimalTest, CallsCallbackForAllThreePhases) {
   bool polishing_termination_is_last = false;
   SolverResult output = PrimalDualHybridGradient(
       lp_, params_, /*interrupt_solve=*/nullptr,
+      /*message_callback=*/nullptr,
       [&](const IterationCallbackInfo& callback_info) {
         ++callback_count[callback_info.iteration_type];
         polishing_termination_is_last =
@@ -2061,6 +2095,36 @@ TEST(PresolveTest, PresolveInfeasible) {
   EXPECT_EQ(output.solve_log.iteration_count(), 0);
 }
 
+TEST(PresolveTest, WarnsInitialSolution) {
+  PrimalAndDualSolution initial_solution;
+  initial_solution.primal_solution.setZero(6);
+  initial_solution.dual_solution.setZero(3);
+  PrimalDualHybridGradientParams params;
+  params.mutable_presolve_options()->set_use_glop(true);
+  CaptureTestStdout();
+  SolverResult output = PrimalDualHybridGradient(CorrelationClusteringStarLp(),
+                                                 params, initial_solution);
+  EXPECT_THAT(GetCapturedTestStdout(),
+              AllOf(HasSubstr("initial solution"), HasSubstr("presolve")));
+}
+
+TEST(PresolveTest, WarnsInitialSolutionViaCallback) {
+  PrimalAndDualSolution initial_solution;
+  initial_solution.primal_solution.setZero(6);
+  initial_solution.dual_solution.setZero(3);
+  PrimalDualHybridGradientParams params;
+  params.mutable_presolve_options()->set_use_glop(true);
+  std::vector<std::string> message_output;
+  SolverResult output = PrimalDualHybridGradient(
+      CorrelationClusteringStarLp(), params, initial_solution,
+      /*interrupt_solve=*/nullptr,
+      [&message_output](const std::string& message) {
+        message_output.push_back(message);
+      });
+  EXPECT_THAT(message_output, Contains(AllOf(HasSubstr("initial solution"),
+                                             HasSubstr("presolve"))));
+}
+
 TEST_P(PresolveDualScalingTest, Dualize) {
   auto [dualize, negate_and_scale_objective] = GetParam();
   PrimalDualHybridGradientParams params;
diff --git a/ortools/pdlp/sharded_optimization_utils.cc b/ortools/pdlp/sharded_optimization_utils.cc
index ee9960d40e..7226edf1e1 100644
--- a/ortools/pdlp/sharded_optimization_utils.cc
+++ b/ortools/pdlp/sharded_optimization_utils.cc
@@ -723,13 +723,26 @@ bool HasValidBounds(const ShardedQuadraticProgram& sharded_qp) {
 }
 
 void ProjectToPrimalVariableBounds(const ShardedQuadraticProgram& sharded_qp,
-                                   VectorXd& primal) {
+                                   VectorXd& primal,
+                                   const bool use_feasibility_bounds) {
+  const auto make_finite_values_zero = [](const double x) {
+    return std::isfinite(x) ? 0.0 : x;
+  };
   sharded_qp.PrimalSharder().ParallelForEachShard(
       [&](const Sharder::Shard& shard) {
         const QuadraticProgram& qp = sharded_qp.Qp();
-        shard(primal) = shard(primal)
-                            .cwiseMin(shard(qp.variable_upper_bounds))
-                            .cwiseMax(shard(qp.variable_lower_bounds));
+        if (use_feasibility_bounds) {
+          shard(primal) =
+              shard(primal)
+                  .cwiseMin(shard(qp.variable_upper_bounds)
+                                .unaryExpr(make_finite_values_zero))
+                  .cwiseMax(shard(qp.variable_lower_bounds)
+                                .unaryExpr(make_finite_values_zero));
+        } else {
+          shard(primal) = shard(primal)
+                              .cwiseMin(shard(qp.variable_upper_bounds))
+                              .cwiseMax(shard(qp.variable_lower_bounds));
+        }
       });
 }
 
@@ -742,6 +755,9 @@ void ProjectToDualVariableBounds(const ShardedQuadraticProgram& sharded_qp,
         const auto upper_bound_shard = shard(qp.constraint_upper_bounds);
         auto dual_shard = shard(dual);
 
+        // TODO(user): Investigate whether it is more efficient to
+        // use .cwiseMax() + .cwiseMin() with unaryExpr(s) that map
+        // upper_bound_shard and lower_bound_shard to appropriate values.
         for (int64_t i = 0; i < dual_shard.size(); ++i) {
           if (!std::isfinite(upper_bound_shard[i])) {
             dual_shard[i] = std::max(dual_shard[i], 0.0);
diff --git a/ortools/pdlp/sharded_optimization_utils.h b/ortools/pdlp/sharded_optimization_utils.h
index 6f36d85a95..8ec23878cc 100644
--- a/ortools/pdlp/sharded_optimization_utils.h
+++ b/ortools/pdlp/sharded_optimization_utils.h
@@ -190,8 +190,11 @@ SingularValueAndIterations EstimateMaximumSingularValueOfConstraintMatrix(
 bool HasValidBounds(const ShardedQuadraticProgram& sharded_qp);
 
 // Projects `primal` onto the variable bounds constraints.
+// If `use_feasibility_bounds == true`, all finite variable bounds are replaced
+// by zero.
 void ProjectToPrimalVariableBounds(const ShardedQuadraticProgram& sharded_qp,
-                                   Eigen::VectorXd& primal);
+                                   Eigen::VectorXd& primal,
+                                   bool use_feasibility_bounds = false);
 
 // Projects `dual` onto the dual variable bounds; see
 // https://developers.google.com/optimization/lp/pdlp_math#dual_variable_bounds.
diff --git a/ortools/pdlp/sharded_optimization_utils_test.cc b/ortools/pdlp/sharded_optimization_utils_test.cc
index 9f3f007759..bd9d63231c 100644
--- a/ortools/pdlp/sharded_optimization_utils_test.cc
+++ b/ortools/pdlp/sharded_optimization_utils_test.cc
@@ -631,6 +631,14 @@ TEST(ProjectToPrimalVariableBoundsTest, TestLp) {
   EXPECT_THAT(primal, ElementsAre(-3, -2, 5, 3.5));
 }
 
+TEST(ProjectToPrimalVariableBoundsTest, TestLpWithFeasibility) {
+  ShardedQuadraticProgram qp(TestLp(), /*num_threads=*/2,
+                             /*num_shards=*/2);
+  Eigen::VectorXd primal{{-3, -3, 5, 5}};
+  ProjectToPrimalVariableBounds(qp, primal, /*use_feasibility_bounds=*/true);
+  EXPECT_THAT(primal, ElementsAre(-3, 0, 0, 0));
+}
+
 TEST(ProjectToDualVariableBoundsTest, TestLp) {
   ShardedQuadraticProgram qp(TestLp(), /*num_threads=*/2,
                              /*num_shards=*/2);
diff --git a/ortools/pdlp/sharded_quadratic_program.cc b/ortools/pdlp/sharded_quadratic_program.cc
index 642d5ade9e..d4f72b3106 100644
--- a/ortools/pdlp/sharded_quadratic_program.cc
+++ b/ortools/pdlp/sharded_quadratic_program.cc
@@ -26,6 +26,7 @@
 #include "ortools/base/threadpool.h"
 #include "ortools/pdlp/quadratic_program.h"
 #include "ortools/pdlp/sharder.h"
+#include "ortools/util/logging.h"
 
 namespace operations_research::pdlp {
 
@@ -35,7 +36,8 @@ namespace {
 // a single column.
 void WarnIfMatrixUnbalanced(
     const Eigen::SparseMatrix<double, Eigen::ColMajor, int64_t>& matrix,
-    absl::string_view matrix_name, int64_t density_limit) {
+    absl::string_view matrix_name, int64_t density_limit,
+    operations_research::SolverLogger* logger) {
   if (matrix.cols() == 0) return;
   int64_t worst_column = 0;
   for (int64_t col = 0; col < matrix.cols(); ++col) {
@@ -46,23 +48,35 @@ void WarnIfMatrixUnbalanced(
   if (matrix.col(worst_column).nonZeros() > density_limit) {
     // TODO(user): fix this automatically in presolve instead of asking the
     // user to do it.
-    LOG(WARNING)
-        << "The " << matrix_name << " has "
-        << matrix.col(worst_column).nonZeros() << " non-zeros in its "
-        << worst_column
-        << "th column. For best parallelization all rows and columns should "
-           "have at most order "
-        << density_limit
-        << " non-zeros. Consider rewriting the QP to split the corresponding "
-           "variable or constraint.";
+    if (logger) {
+      SOLVER_LOG(
+          logger, "WARNING: The ", matrix_name, " has ",
+          matrix.col(worst_column).nonZeros(), " non-zeros in its ",
+          worst_column,
+          "th column. For best parallelization all rows and columns should "
+          "have at most order ",
+          density_limit,
+          " non-zeros. Consider rewriting the QP to split the corresponding "
+          "variable or constraint.");
+    } else {
+      LOG(WARNING)
+          << "The " << matrix_name << " has "
+          << matrix.col(worst_column).nonZeros() << " non-zeros in its "
+          << worst_column
+          << "th column. For best parallelization all rows and columns should "
+             "have at most order "
+          << density_limit
+          << " non-zeros. Consider rewriting the QP to split the corresponding "
+             "variable or constraint.";
+    }
   }
 }
 
 }  // namespace
 
-ShardedQuadraticProgram::ShardedQuadraticProgram(QuadraticProgram qp,
-                                                 const int num_threads,
-                                                 const int num_shards)
+ShardedQuadraticProgram::ShardedQuadraticProgram(
+    QuadraticProgram qp, const int num_threads, const int num_shards,
+    operations_research::SolverLogger* logger)
     : qp_(std::move(qp)),
       transposed_constraint_matrix_(qp_.constraint_matrix.transpose()),
       thread_pool_(num_threads == 1
@@ -85,10 +99,10 @@ ShardedQuadraticProgram::ShardedQuadraticProgram(QuadraticProgram qp,
                                        qp_.constraint_lower_bounds.size();
     const int64_t column_density_limit = work_per_iteration / num_threads;
     WarnIfMatrixUnbalanced(qp_.constraint_matrix, "constraint matrix",
-                           column_density_limit);
+                           column_density_limit, logger);
     WarnIfMatrixUnbalanced(transposed_constraint_matrix_,
-                           "transposed constraint matrix",
-                           column_density_limit);
+                           "transposed constraint matrix", column_density_limit,
+                           logger);
   }
 }
 
diff --git a/ortools/pdlp/sharded_quadratic_program.h b/ortools/pdlp/sharded_quadratic_program.h
index 0e0438dd5f..9a80ec5f3a 100644
--- a/ortools/pdlp/sharded_quadratic_program.h
+++ b/ortools/pdlp/sharded_quadratic_program.h
@@ -24,6 +24,7 @@
 #include "ortools/base/threadpool.h"
 #include "ortools/pdlp/quadratic_program.h"
 #include "ortools/pdlp/sharder.h"
+#include "ortools/util/logging.h"
 
 namespace operations_research::pdlp {
 
@@ -37,7 +38,10 @@ class ShardedQuadraticProgram {
  public:
   // Requires `num_shards` >= `num_threads` >= 1.
   // Note that the `qp` is intentionally passed by value.
-  ShardedQuadraticProgram(QuadraticProgram qp, int num_threads, int num_shards);
+  // If `logger` is not nullptr, warns about unbalanced matrices using it;
+  // otherwise warns via Google standard logging.
+  ShardedQuadraticProgram(QuadraticProgram qp, int num_threads, int num_shards,
+                          operations_research::SolverLogger* logger = nullptr);
 
   // Movable but not copyable.
   ShardedQuadraticProgram(const ShardedQuadraticProgram&) = delete;
diff --git a/ortools/pdlp/solve_log.proto b/ortools/pdlp/solve_log.proto
index 37bc8a6174..282048499a 100644
--- a/ortools/pdlp/solve_log.proto
+++ b/ortools/pdlp/solve_log.proto
@@ -204,12 +204,13 @@ message ConvergenceInformation {
 // infeasibility (i.e. has no solution); see also TerminationCriteria.
 message InfeasibilityInformation {
   // Let x_ray be the algorithm's estimate of the primal extreme ray where x_ray
-  // is a vector scaled such that its infinity norm is one. A simple and typical
+  // is a vector that satisfies the sign constraints for a ray, scaled such that
+  // its infinity norm is one (the sign constraints are the variable bound
+  // constraints, with all finite bounds mapped to zero). A simple and typical
   // choice of x_ray is x_ray = x / | x |_∞ where x is the current primal
-  // iterate. For this value compute the maximum absolute error in the primal
-  // linear program with the right hand side and finite variable bounds set to
-  // zero. This error refers to both the linear constraints and sign constraints
-  // on the ray.
+  // iterate projected onto the primal ray sign constraints. For this value
+  // compute the maximum absolute error in the primal linear program with the
+  // right hand side set to zero.
   optional double max_primal_ray_infeasibility = 1;
 
   // The value of the linear part of the primal objective (ignoring additive
@@ -223,12 +224,12 @@ message InfeasibilityInformation {
   optional double primal_ray_quadratic_norm = 3;
 
   // Let (y_ray, r_ray) be the algorithm's estimate of the dual and reduced cost
-  // extreme ray where (y_ray, r_ray) is a vector scaled such that its infinity
-  // norm is one. A simple and typical choice of y_ray is
-  // (y_ray, r_ray) = (y, r) / max(| y |_∞, | r |_∞) where y is the
-  // current dual iterate and r is the current dual reduced costs.
-  // Consider the quadratic program we are solving but with the objective (both
-  // quadratic and linear terms) set to zero. This forms a linear program
+  // extreme ray where (y_ray, r_ray) is a vector (satisfying the dual variable
+  // constraints) scaled such that its infinity norm is one. A simple and
+  // typical choice of y_ray is (y_ray, r_ray) = (y, r) / max(| y |_∞, | r |_∞)
+  // where y is the current dual iterate and r is the current dual reduced
+  // costs. Consider the quadratic program we are solving but with the objective
+  // (both quadratic and linear terms) set to zero. This forms a linear program
   // (label this linear program (1)) with no objective. Take the dual of (1) and
   // compute the maximum absolute value of the constraint error for
   // (y_ray, r_ray) to obtain the value of max_dual_ray_infeasibility.