wip

Author: Yanyu000

check/TestPdlp.cpp

@@ -344,7 +344,7 @@ TEST_CASE("hi-pdlp", "[pdlp]") {
   h.setOptionValue("pdlp_restart_strategy", 0);
   h.setOptionValue("pdlp_step_size_strategy", 0);
 
-  h.setOptionValue("pdlp_iteration_limit", 1000000);
+  h.setOptionValue("pdlp_iteration_limit", 15);
   //  h.setOptionValue("log_dev_level", kHighsLogDevLevelVerbose);
   HighsStatus run_status = h.run();
   //  REQUIRE(run_status == HighsStatus::kOk);

highs/pdlp/cupdlp/cupdlp_step.c

@@ -18,6 +18,11 @@ void PDHG_primalGradientStep(CUPDLPwork *work, CUPDLPvec *xUpdate,
                              cupdlp_float dPrimalStepSize) {
   CUPDLPproblem *problem = work->problem;
 
+// print norm of A'y
+double aty_norm = 0.0;
+cupdlp_twoNorm(work, problem->nCols, ATy->data, &aty_norm);
+cupdlp_printf("Norm of A'y: %g\n", aty_norm);
+
 #if !defined(CUPDLP_CPU) && USE_KERNELS
   cupdlp_pgrad_cuda(xUpdate->data, x->data, problem->cost,
                     ATy->data, problem->lower, problem->upper, dPrimalStepSize,
@@ -192,6 +197,11 @@ void PDHG_Update_Iterate_Constant_Step_Size(CUPDLPwork *pdhg) {
   CUPDLPvec *atyUpdate = iterates->aty[(iter + 1) % 2];
 
   Ax(pdhg, ax, x);
+  //pint norm of Ax
+  double ax_norm = 0.0;
+  cupdlp_twoNorm(pdhg, problem->nRows, ax->data,
+                  &ax_norm);
+  cupdlp_printf("Norm of Ax: %g\n", ax_norm);
   ATy(pdhg, aty, y);
 
   // x^{k+1} = proj_{X}(x^k - dPrimalStep * (c - A'y^k))

highs/pdlp/hipdlp/pdhg.cc

@@ -350,8 +350,7 @@ void PDLPSolver::Solve(HighsLp& original_lp, const PrimalDualParams& params,
   // We need Ax for the y-update and for computing residuals
   std::vector<double> Ax_current(lp.num_row_, 0.0);
   std::vector<double> ATy_current(lp.num_col_, 0.0);
-  linalg::Ax(lp, x_current_, Ax_current);
-  linalg::ATy(lp, y_current_, ATy_current);
+
   std::vector<double> Ax_new(lp.num_row_, 0.0);
   std::vector<double> ATy_new(lp.num_col_, 0.0);
 
@@ -364,6 +363,12 @@ void PDLPSolver::Solve(HighsLp& original_lp, const PrimalDualParams& params,
   // A single loop handles max iterations, convergence, and restarts.
   for (int iter = 0; iter < params.max_iterations; ++iter) {
     pdlpLog(pdlp_log_file, iter);
+    linalg::Ax(lp, x_current_, Ax_current);
+    //print norm of Ax
+    double ax_norm = linalg::vector_norm(Ax_current);
+    std::cout << "Norm of Ax: " << ax_norm << std::endl;
+
+    linalg::ATy(lp, y_current_, ATy_current);
     if (solver_timer.read() > params.time_limit) {
       logger_.info("Time limit reached.");
       final_iter_count_ = iter;

highs/pdlp/hipdlp/step.cc

@@ -89,7 +89,9 @@ void UpdateX(std::vector<double>& x_new, const std::vector<double>& x_current,
              double eta, double omega) {
   std::vector<double> ATy_cache(lp.num_col_);
   linalg::ATy(lp, y_current, ATy_cache);
-
+  //print the norm of ATy_cache
+  double aty_norm = linalg::vector_norm(ATy_cache);
+  std::cout << "Norm of A'y: " << aty_norm << std::endl;
   for (HighsInt i = 0; i < lp.num_col_; i++) {
     double gradient = lp.col_cost_[i] - ATy_cache[i];
     x_new[i] = linalg::project_box(x_current[i] - (eta / omega) * gradient,