add projection before

Author: Yanyu000

highs/pdlp/hipdlp/linalg.cc

@@ -21,6 +21,19 @@ double project_non_negative(double x) {
   return std::max(0.0, x);
 }
 
+void project_bounds(const HighsLp& lp, std::vector<double>& x) {
+  for (HighsInt i = 0; i < lp.num_col_; ++i) {
+    // Project to upper bound
+    if (x[i] > lp.col_upper_[i]) {
+      x[i] = lp.col_upper_[i];
+    }
+    // Project to lower bound
+    if (x[i] < lp.col_lower_[i]) {
+      x[i] = lp.col_lower_[i];
+    }
+  }
+}
+
 void Ax(const HighsLp& lp, const std::vector<double>& x,
         std::vector<double>& result) {
   std::fill(result.begin(), result.end(), 0.0);

highs/pdlp/hipdlp/linalg.hpp

@@ -18,6 +18,7 @@
 namespace linalg {
 double project_box(double x, double l, double u);
 double project_non_negative(double y);
+void project_bounds(const HighsLp& lp, std::vector<double>& x);
 
 // Function to compute A*x for a given HighsLp and vector x
 void Ax(const HighsLp& lp, const std::vector<double>& x,

highs/pdlp/hipdlp/pdhg.cc

@@ -377,17 +377,20 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
   // Use member variables for iterates for cleaner state management
   x_current_ = x;
   y_current_ = y;
+  linalg::project_bounds(lp_, x_current_);
+  x_sum_ = x_current_;
+  x_avg_ = x_current_;
+  linalg::Ax(lp, x_current_, Ax_cache_);
+  std::vector<double> Ax_avg = Ax_cache_;
+  std::vector<double> ATy_avg(lp.num_col_, 0.0);
+  
   num_rejected_steps_ = 0;
   bool first_malitsky_iteration = true;
   ratio_last_two_step_sizes_ = 1.0;
   bool using_malitsky_averaging =
       (params_.step_size_strategy == StepSizeStrategy::MALITSKY_POCK);
   bool primal_average_initialized = false; 
 
-  // Initialize vectors for matrix-vector products
-  std::vector<double> Ax_avg(lp.num_row_, 0.0);
-  std::vector<double> ATy_avg(lp.num_col_, 0.0);
-
   // Store iterates at last restart for primal weight update
   x_at_last_restart_ = x_current_;
   y_at_last_restart_ = y_current_;
@@ -400,7 +403,7 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
   debug_pdlp_data_.ax_average_norm = 0.0;
   debug_pdlp_data_.aty_average_norm = 0.0;
   debug_pdlp_data_.x_average_norm = 0.0;
-  debug_pdlp_data_.ax_norm = 0.0;
+  debug_pdlp_data_.ax_norm = linalg::vector_norm(Ax_cache_);
 
   for (int iter = 0; iter < params_.max_iterations; ++iter) {
     debugPdlpIterLog(debug_pdlp_log_file_, iter, &debug_pdlp_data_,