Added pdlpLog to cupdlp_utils.c and calling it from both cupdlp_solve.c and pdhg.cpp

Author: jajhall

check/TestPdlp.cpp

@@ -338,7 +338,7 @@ TEST_CASE("hi-pdlp", "[pdlp]") {
   //  h.setOptionValue("output_flag", dev_run);
   REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
   h.setOptionValue("solver", kHiPdlpString);
-  h.setOptionValue("kkt_tolerance", 1e-6);  // kkt_tolerance);
+  h.setOptionValue("kkt_tolerance", kkt_tolerance);
 
   h.setOptionValue("pdlp_scaling_mode", 0);
   h.setOptionValue("pdlp_restart_strategy", 0);

highs/pdlp/cupdlp/cupdlp_solver.c

@@ -897,6 +897,9 @@ cupdlp_retcode PDHG_Solve(const cupdlp_int* has_variables, CUPDLPwork *pdhg) {
   timers->nIter = 0;
   timers->dSolvingBeg = getTimeStamp();
 
+  FILE* pdlp_log_file = fopen("cuPDLP.log", "w");
+  assert(pdlp_log_file);
+
   // PDHG_Init_Data does nothing!
   PDHG_Init_Data(pdhg);
 
@@ -915,6 +918,7 @@ cupdlp_retcode PDHG_Solve(const cupdlp_int* has_variables, CUPDLPwork *pdhg) {
   const int iter_log_between_header = 50;
   int iter_log_since_header = iter_log_between_header;
   for (timers->nIter = 0; timers->nIter < settings->nIterLim; ++timers->nIter) {
+    pdlpLog(pdlp_log_file, timers->nIter);
     PDHG_Compute_SolvingTime(pdhg);
 #if CUPDLP_DUMP_ITERATES_STATS && CUPDLP_DEBUG
     PDHG_Dump_Stats(pdhg);
@@ -1150,6 +1154,7 @@ cupdlp_retcode PDHG_Solve(const cupdlp_int* has_variables, CUPDLPwork *pdhg) {
 #endif
 
 exit_cleanup:
+  fclose(pdlp_log_file);
   return retcode;
 }
 

highs/pdlp/cupdlp/cupdlp_utils.c

@@ -1781,3 +1781,7 @@ void writeSol(const char *fout, cupdlp_int nCols, cupdlp_int nRows,
   // Close the file
   fclose(fptr);
 }
+
+void pdlpLog(FILE* file, const int iter_num) {
+  fprintf(file, "%6d\n", iter_num);
+}

highs/pdlp/cupdlp/cupdlp_utils.h

@@ -193,6 +193,8 @@ void writeSol(const char *fout, cupdlp_int nCols, cupdlp_int nRows,
               cupdlp_float *col_value, cupdlp_float *col_dual,
               cupdlp_float *row_value, cupdlp_float *row_dual);
 
+void pdlpLog(FILE* file, const int iter_num);
+
 #ifdef __cplusplus
 }
 #endif

highs/pdlp/hipdlp/pdhg.cc

@@ -1,15 +1,15 @@
 #include "pdhg.hpp"
 
-#include "linalg.hpp"
-#include "restart.hpp"
-#include "step.hpp"
-// #include "Highs.h"
-// #include "highs_interface.hpp"
 #include <cmath>
 #include <iostream>
 #include <random>
 #include <tuple>
 
+#include "linalg.hpp"
+#include "pdlp/cupdlp/cupdlp.h"  // For pdlpLogging
+#include "restart.hpp"
+#include "step.hpp"
+
 #define PDHG_CHECK_INTERVAL 40
 
 using namespace std;
@@ -266,7 +266,7 @@ void PDLPSolver::Solve(HighsLp& original_lp, const PrimalDualParams& params,
   Timer solver_timer;
   if (original_lp.num_row_ == 0 || original_lp.a_matrix_.start_.size() < 2) {
     std::cout << "Empty LP problem." << std::endl;
-    return;
+    return solveReturn();
   }
 
   HighsSolution solution;
@@ -323,8 +323,7 @@ void PDLPSolver::Solve(HighsLp& original_lp, const PrimalDualParams& params,
 
   working_params.omega = std::sqrt(step_size.dual_step / step_size.primal_step);
   working_params.eta = std::sqrt(step_size.primal_step * step_size.dual_step);
-  current_eta_ =
-      working_params.eta;  // Initial step size for adaptive strategy
+  current_eta_ = working_params.eta;  // Initial step size for adaptive strategy
   std::cout << "Using power method step sizes: eta = " << working_params.eta
             << ", omega = " << working_params.omega << std::endl;
 
@@ -358,9 +357,13 @@ void PDLPSolver::Solve(HighsLp& original_lp, const PrimalDualParams& params,
 
   logger_.print_iteration_header();
 
+  pdlp_log_file = fopen("HiPDLP.log", "w");
+  assert(pdlp_log_file);
+
   // --- 2. Main PDHG Loop ---
   // A single loop handles max iterations, convergence, and restarts.
   for (int iter = 0; iter < params.max_iterations; ++iter) {
+    pdlpLog(pdlp_log_file, iter);
     if (solver_timer.read() > params.time_limit) {
       logger_.info("Time limit reached.");
       final_iter_count_ = iter;
@@ -381,7 +384,7 @@ void PDLPSolver::Solve(HighsLp& original_lp, const PrimalDualParams& params,
       }
       results_.primal_obj = primal_objective;
 
-      return;  // Exit the function
+      return solveReturn();  // Exit the function
     }
 
     // --- 3. Core PDHG Update Step ---
@@ -412,7 +415,7 @@ void PDLPSolver::Solve(HighsLp& original_lp, const PrimalDualParams& params,
           x = x_avg_;
           y = y_avg_;
           // scaling_.UnscaleSolution(x, y);
-          return;
+          return solveReturn();
         }
     }
 
@@ -463,7 +466,7 @@ void PDLPSolver::Solve(HighsLp& original_lp, const PrimalDualParams& params,
         results_ = average_results;
         results_.term_code = TerminationStatus::OPTIMAL;
         Postsolve(original_lp, lp, x, y, solution);
-        return;
+        return solveReturn();
       }
 
       bool current_converged = CheckConvergence(
@@ -491,7 +494,7 @@ void PDLPSolver::Solve(HighsLp& original_lp, const PrimalDualParams& params,
         results_ = current_results;
         results_.term_code = TerminationStatus::OPTIMAL;
         Postsolve(original_lp, lp, x, y, solution);
-        return;
+        return solveReturn();
       }
 
       // --- 6. Restart Check ---
@@ -551,7 +554,11 @@ void PDLPSolver::Solve(HighsLp& original_lp, const PrimalDualParams& params,
   results_.primal_obj = primal_objective;
   results_.term_code = TerminationStatus::TIMEOUT;
   Postsolve(original_lp, lp, x_avg_, y_avg_, solution);
-  return;
+  return solveReturn();
+}
+
+void PDLPSolver::solveReturn() {
+  if (pdlp_log_file) fclose(pdlp_log_file);
 }
 
 void PDLPSolver::Initialize(const HighsLp& lp, std::vector<double>& x,
@@ -792,10 +799,10 @@ bool PDLPSolver::CheckConvergence(const HighsLp& lp,
       ComputeDualityGap(lp, x, y, lambda);
   results.duality_gap = duality_gap;
 
-  results.relative_obj_gap = std::abs(duality_gap) / (1.0 + std::abs(cTx) +
-                                                       std::abs(qTy) +
-                                                       std::abs(lTlambda_plus) +
-                                                       std::abs(uTlambda_minus));
+  results.relative_obj_gap =
+      std::abs(duality_gap) /
+      (1.0 + std::abs(cTx) + std::abs(qTy) + std::abs(lTlambda_plus) +
+       std::abs(uTlambda_minus));
 
   bool primal_feasible = primal_feasibility <= epsilon * (1 + q_norm);
   bool dual_feasible = dual_feasibility <= epsilon * (1 + c_norm);

highs/pdlp/hipdlp/pdhg.hpp

@@ -25,7 +25,6 @@
 #include "scaling.hpp"
 #include "solver_results.hpp"
 #include "step.hpp"
-// #include "highs_interface.hpp"
 
 // --- Define Macros ---
 // Enable or disable GPU usage
@@ -85,6 +84,8 @@ class PDLPSolver {
   ScalingParams scaling_params_;
   Scaling scaling_;
 
+  FILE* pdlp_log_file = nullptr;
+
   // HighsStatus TransformGxLeqToGeq(HighsLp& lp);
   void PreprocessLp(const HighsLp& original_lp, HighsLp& processed_lp);
   void Postsolve(const HighsLp& original_lp, HighsLp& processed_lp,
@@ -93,6 +94,7 @@ class PDLPSolver {
                  HighsSolution& solution);
 
   // Helper functions
+  void solveReturn();
   void Initialize(const HighsLp& lp, std::vector<double>& x,
                   std::vector<double>& y);
   bool RunPresolve(const HighsLp& original_lp, Highs& highs,