revert to same result

Author: Yanyu000

check/TestPdlp.cpp

@@ -336,7 +336,7 @@ TEST_CASE("pdlp-restart-add-row", "[pdlp]") {
 
 TEST_CASE("hi-pdlp", "[pdlp]") {
   std::string model =
-      "adlittle";  //"adlittle";//"afiro";// shell// stair //25fv47 //fit2p
+      "25fv47";  //"adlittle";//"afiro";// shell// stair //25fv47 //fit2p
   std::string model_file =
       std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
   Highs h;
@@ -359,7 +359,7 @@ TEST_CASE("hi-pdlp", "[pdlp]") {
   h.setOptionValue("pdlp_scaling_mode", pdlp_scaling);
   h.setOptionValue("pdlp_step_size_strategy", 1);
   h.setOptionValue("pdlp_restart_strategy", 2);
-  h.setOptionValue("pdlp_iteration_limit", 3520);
+  // h.setOptionValue("pdlp_iteration_limit", 3520);
   //  h.setOptionValue("log_dev_level", kHighsLogDevLevelVerbose);
   auto start_hipdlp = std::chrono::high_resolution_clock::now();
   HighsStatus run_status = h.run();

highs/pdlp/cupdlp/cupdlp_defs.h

@@ -41,7 +41,7 @@ extern "C" {
 #define CUPDLP_DUMP_LINESEARCH_STATS (1)
 #define CUPDLP_INEXACT_EPS (1e-4)
 
-#define PDLP_DEBUG_LOG (0)
+#define PDLP_DEBUG_LOG (1)
 
 typedef struct CUPDLP_CUDA_DENSE_VEC CUPDLPvec;
 typedef struct CUPDLP_DENSE_MATRIX CUPDLPdense;

highs/pdlp/hipdlp/pdhg.cc

@@ -164,6 +164,7 @@ void PDLPSolver::preprocessLp() {
 
   int num_new_cols = 0;
   int nEqs = 0;
+  sense_origin_ = (original_lp_->sense_ == ObjSense::kMinimize) ? 1 : -1;
   constraint_types_.resize(nRows_orig);
 
   // 1. First pass: Classify constraints and count slack variables needed
@@ -433,21 +434,27 @@ PostSolveRetcode PDLPSolver::postprocess(HighsSolution& solution) {
     }
   }
 
-  std::vector<double> ax_current_(unscaled_processed_lp_.num_row_);
-  linalg::Ax(unscaled_processed_lp_, x_current_, ax_current_);
+  std::vector<double> ax_original(original_lp_->num_row_, 0.0);
+
+  // Get original matrix in column-wise format
+  HighsSparseMatrix orig_matrix = original_lp_->a_matrix_;
+  orig_matrix.ensureColwise();
+
+  // Compute Ax using only the original columns (not slack variables)
+  for (int col = 0; col < original_num_col_; ++col) {
+      double x_val = x_current_[col];  // Use unscaled x values
+      
+      for (int el = orig_matrix.start_[col]; 
+          el < orig_matrix.start_[col + 1]; ++el) {
+          int row = orig_matrix.index_[el];
+          double a_val = orig_matrix.value_[el];
+          ax_original[row] += a_val * x_val;
+      }
+  }
 
-  int slack_variable_idx = original_num_col_;
+  // Now ax_original contains the correct row activity values
   for (int orig_row = 0; orig_row < original_lp_->num_row_; ++orig_row) {
-    int reordered_row = constraint_new_idx_[orig_row];
-
-    if (constraint_types_[orig_row] == BOUND ||
-        constraint_types_[orig_row] == FREE) {
-      solution.row_value[orig_row] = x_current_[slack_variable_idx++];
-    } else if (constraint_types_[orig_row] == LEQ) {
-      solution.row_value[orig_row] = -ax_current_[reordered_row];
-    } else {  // EQ, GEQ
-      solution.row_value[orig_row] = ax_current_[reordered_row];
-    }
+      solution.row_value[orig_row] = ax_original[orig_row];
   }
 
   // 6. Recover Dual Column Values (Reduced Costs)
@@ -461,6 +468,7 @@ PostSolveRetcode PDLPSolver::postprocess(HighsSolution& solution) {
       return PostSolveRetcode::DIMENSION_MISMATCH;
     }
     solution.col_dual[i] = dSlackPos_[i] - dSlackNeg_[i];
+    solution.col_dual[i] *= sense_origin_;
   }
 
   solution.value_valid = true;  // to do
@@ -635,8 +643,12 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
 
         restart_scheme_.last_restart_iter_ = iter;
         // Recompute Ax and ATy for the restarted iterates
+        hipdlpTimerStart(kHipdlpClockMatrixMultiply);
         linalg::Ax(lp, x_current_, Ax_cache_);
+        hipdlpTimerStop(kHipdlpClockMatrixMultiply);
+        hipdlpTimerStart(kHipdlpClockMatrixTransposeMultiply);
         linalg::ATy(lp, y_current_, ATy_cache_);
+        hipdlpTimerStop(kHipdlpClockMatrixTransposeMultiply);
 
         restart_scheme_.SetLastRestartIter(iter);
       }
@@ -679,9 +691,12 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
 
     // Compute ATy for the new iterate
     Ax_cache_ = Ax_next_;
+    /*
     hipdlpTimerStart(kHipdlpClockMatrixTransposeMultiply);
     linalg::ATy(lp, y_next_, ATy_cache_);
     hipdlpTimerStop(kHipdlpClockMatrixTransposeMultiply);
+    */
+    ATy_cache_ = ATy_next_;
 
     hipdlpTimerStop(kHipdlpClockIterateUpdate);
 
@@ -866,16 +881,23 @@ void PDLPSolver::computeDualSlacks(const std::vector<double>& ATy_vector) {
   if (dSlackNeg_.size() != lp_.num_col_) dSlackNeg_.resize(lp_.num_col_);
 
   for (HighsInt i = 0; i < lp_.num_col_; ++i) {
-    double reduced_cost = lp_.col_cost_[i] - ATy_vector[i];
-
-    // Slack for lower bound l_i <= x_i (dual variable is non-negative)
-    dSlackPos_[i] =
-        (lp_.col_lower_[i] > -kHighsInf) ? std::max(0.0, reduced_cost) : 0.0;
+    double dual_residual = lp_.col_cost_[i] - ATy_vector[i];
 
-    // Slack for upper bound x_i <= u_i (dual variable is non-positive, so slack
-    // is non-negative)
-    dSlackNeg_[i] =
-        (lp_.col_upper_[i] < kHighsInf) ? std::max(0.0, -reduced_cost) : 0.0;
+    // Compute positive slack (for lower bounds)
+        // CUPDLP: max(dual_residual, 0) * hasLower
+        if (lp_.col_lower_[i] > -kHighsInf) {
+            dSlackPos_[i] = std::max(0.0, dual_residual);
+        } else {
+            dSlackPos_[i] = 0.0;
+        }
+        
+        // Compute negative slack (for upper bounds)
+        // CUPDLP: -min(dual_residual, 0) * hasUpper
+        if (lp_.col_upper_[i] < kHighsInf) {
+            dSlackNeg_[i] = std::max(0.0, -dual_residual);
+        } else {
+            dSlackNeg_[i] = 0.0;
+        }
   }
 }
 
@@ -1239,7 +1261,6 @@ void PDLPSolver::setup(const HighsOptions& options, HighsTimer& timer) {
   restart_scheme_.passParams(&params_);
   // Copy what's needed to use HiGHS logging
   params_.log_options_ = options.log_options;
-
   // log the options
   logger_.print_params(params_);
 }
@@ -1249,9 +1270,6 @@ void PDLPSolver::scaleProblem() {
   scaling_.passParams(&params_);
   scaling_.Initialize(lp_);
   scaling_.scaleProblem();
-  //print the norm of lp_.col_cost_
-  unscaled_c_norm_ = linalg::vector_norm(lp_.col_cost_);
-  std::cout << "Unscaled cost norm = " << unscaled_c_norm_ << std::endl;
 }
 
 void PDLPSolver::unscaleSolution(std::vector<double>& x,
@@ -1264,8 +1282,10 @@ void PDLPSolver::unscaleSolution(std::vector<double>& x,
   const std::vector<double>& col_scale = scaling_.GetColScaling();
   if (!dSlackPos_.empty() && col_scale.size() == dSlackPos_.size()) {
     for (size_t i = 0; i < dSlackPos_.size(); ++i) {
-      dSlackPos_[i] /= col_scale[i];
-      dSlackNeg_[i] /= col_scale[i];
+      std::cout << "dSlackPos_ before unscale[" << i << "] = " << dSlackPos_[i]
+                << std::endl;
+      dSlackPos_[i] *= col_scale[i];
+      dSlackNeg_[i] *= col_scale[i];
     }
   }
 }
@@ -1403,6 +1423,10 @@ void PDLPSolver::updateIteratesFixed() {
   hipdlpTimerStart(kHipdlpClockProjectY);
   y_next_ = updateY(y_current_, Ax_cache_, Ax_next_, stepsize_.dual_step);
   hipdlpTimerStop(kHipdlpClockProjectY);
+
+  hipdlpTimerStart(kHipdlpClockMatrixTransposeMultiply);
+  linalg::ATy(lp_, y_next_, ATy_next_);
+  hipdlpTimerStop(kHipdlpClockMatrixTransposeMultiply);
 }
 
 void PDLPSolver::updateIteratesAdaptive() {

highs/pdlp/hipdlp/pdhg.hpp

@@ -121,6 +121,7 @@ class PDLPSolver {
   std::vector<bool> is_equality_row_;
   std::vector<int> constraint_new_idx_;
   std::vector<ConstraintType> constraint_types_;
+  int sense_origin_ = 1;
 
   // --- Solver State ---
   int final_iter_count_ = 0;