Very sensitive to power method result; logging of HiPDLP primal and dual step lengths probably incorrect

Author: jajhall

check/TestPdlp.cpp

@@ -339,16 +339,19 @@ TEST_CASE("hi-pdlp", "[pdlp]") {
   REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
   h.setOptionValue("solver", kHiPdlpString);
   h.setOptionValue("kkt_tolerance", kkt_tolerance);
+
   h.setOptionValue("pdlp_scaling_mode", 0);
   h.setOptionValue("pdlp_restart_strategy", 0);
   h.setOptionValue("pdlp_step_size_strategy", 0);
-  //  h.setOptionValue("pdlp_iteration_limit", 10000);
+
+  h.setOptionValue("pdlp_iteration_limit", 100000);
   //  h.setOptionValue("log_dev_level", kHighsLogDevLevelVerbose);
   HighsStatus run_status = h.run();
   //  REQUIRE(run_status == HighsStatus::kOk);
   //  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
-  const bool cupdlp_test = false;
+  const bool cupdlp_test = true;
   if (cupdlp_test) {
+    h.clearSolver();
     h.setOptionValue("solver", kCuPdlpString);
     h.setOptionValue("pdlp_features_off", kPdlpAllFeaturesOff);
     run_status = h.run();

highs/pdlp/cupdlp/cupdlp_step.c

@@ -95,8 +95,8 @@ cupdlp_retcode PDHG_Power_Method(CUPDLPwork *work, cupdlp_float *lambda) {
   //
   // y is stored in aty->data
   //
-  int log_iters = work->settings->nLogLevel > 0;
-  log_iters = 1;
+  int log_iters = work->settings->nLogLevel > 1;
+  //  log_iters = 1;
 
   if (log_iters)
     cupdlp_printf("It       lambda   dl_lambda    residual\n");
@@ -318,6 +318,7 @@ cupdlp_retcode PDHG_Init_Step_Sizes(CUPDLPwork *pdhg) {
 
   if (stepsize->eLineSearchMethod == PDHG_FIXED_LINESEARCH) {
     CUPDLP_CALL(PDHG_Power_Method(pdhg, &stepsize->dPrimalStep));
+    cupdlp_float power_method_lambda = stepsize->dPrimalStep;
     // PDLP Intial primal weight = norm(cost) / norm(rhs) = sqrt(beta)
     // cupdlp_float a = twoNormSquared(problem->cost, problem->nCols);
     // cupdlp_float b = twoNormSquared(problem->rhs, problem->nRows);
@@ -336,6 +337,8 @@ cupdlp_retcode PDHG_Init_Step_Sizes(CUPDLPwork *pdhg) {
     stepsize->dDualStep = stepsize->dPrimalStep;
     stepsize->dPrimalStep /= sqrt(stepsize->dBeta);
     stepsize->dDualStep *= sqrt(stepsize->dBeta);
+    cupdlp_printf("Initial step sizes from power method lambda = %g: primal = %g; dual = %g\n",
+		  power_method_lambda, stepsize->dPrimalStep, stepsize->dDualStep);	  
   } else {
     stepsize->dTheta = 1.0;
 

highs/pdlp/hipdlp/pdhg.cc

@@ -317,7 +317,8 @@ void PDLPSolver::Solve(HighsLp & original_lp, const PrimalDualParams& params, st
     working_params.omega = omega;
     working_params.eta = fixed_eta;
     std::cout << "Using power method step sizes: eta = " << fixed_eta << ", omega = " << omega << std::endl;
-
+    printf("Initial step sizes from power method lambda = %g: primal = %g; dual = %g\n",
+	   op_norm_sq, fixed_eta/omega, fixed_eta*omega);
     // --- 1. Initialization ---
     Initialize(lp, x, y); // Sets initial x, y and results_
     restart_scheme_.Initialize(params, results_);
@@ -805,7 +806,7 @@ HighsStatus PDLPSolver::PowerMethod(HighsLp &lp, double& op_norm_sq) {
     const HighsInt kATAPowerMethod = 1;
     const HighsInt kCuPdlpAATPowerMethod = 2;
 
-    const HighsInt power_method = kCuPdlpAATPowerMethod;
+    const HighsInt power_method = kYanyuPowerMethod;//kCuPdlpAATPowerMethod;//
     // Allocate memory for matrix-vector products
     std::vector<double> y_vec;
     std::vector<double> z_vec;
@@ -820,7 +821,7 @@ HighsStatus PDLPSolver::PowerMethod(HighsLp &lp, double& op_norm_sq) {
     double op_norm_sq_old = 0.0;
     LogLevel log_level = logger_.getLogLevel();
     int log_iters = log_level == LogLevel::kVerbose || log_level == LogLevel::kDebug;
-    log_iters = 1;
+    //    log_iters = 1;
 
     if (log_iters) printf("It       lambda   dl_lambda\n");