Now logging after calls to CheckConvergence

Author: jajhall

check/TestPdlp.cpp

@@ -333,7 +333,7 @@ TEST_CASE("pdlp-restart-add-row", "[pdlp]") {
 }
 
 TEST_CASE("hi-pdlp", "[pdlp]") {
-  std::string model = "afiro";  //"adlittle";//"afiro";// 
+  std::string model = "afiro";  //"adlittle";//"afiro";//
   std::string model_file =
       std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
   Highs h;
@@ -342,17 +342,15 @@ TEST_CASE("hi-pdlp", "[pdlp]") {
   h.setOptionValue("solver", kHiPdlpString);
   h.setOptionValue("kkt_tolerance", kkt_tolerance);
 
-  HighsInt pdlp_features_off = 
-  // kPdlpScalingOff + 
-  kPdlpRestartOff + 
-  kPdlpAdaptiveStepSizeOff;
+  HighsInt pdlp_features_off =
+      // kPdlpScalingOff +
+      kPdlpRestartOff + kPdlpAdaptiveStepSizeOff;
   h.setOptionValue("pdlp_features_off", pdlp_features_off);
 
-  HighsInt pdlp_scaling = //0;
-   kPdlpScalingRuiz 
-  //+ kPdlpScalingL2 
-  + kPdlpScalingPC
-  ;
+  HighsInt pdlp_scaling =  // 0;
+      kPdlpScalingRuiz
+      //+ kPdlpScalingL2
+      + kPdlpScalingPC;
   h.setOptionValue("pdlp_scaling_mode", pdlp_scaling);
   h.setOptionValue("pdlp_step_size_strategy", 0);
   h.setOptionValue("pdlp_restart_strategy", 2);

highs/pdlp/CupdlpWrapper.cpp

@@ -207,10 +207,8 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer,
 
   // Print final solution using debugPdlpFinalSolutionLog
   debugPdlpFinalSolutionLog(w->debug_pdlp_log_file_,
-			    highs_solution.col_value.data(),
-			    lp.num_col_,
-			    highs_solution.row_dual.data(),
-			    lp.num_row_);
+                            highs_solution.col_value.data(), lp.num_col_,
+                            highs_solution.row_dual.data(), lp.num_row_);
   fclose(w->debug_pdlp_log_file_);
   model_status = HighsModelStatus::kUnknown;
   highs_solution.value_valid = value_valid;

highs/pdlp/HiPdlpWrapper.cpp

@@ -11,10 +11,10 @@
  */
 #include "pdlp/HiPdlpWrapper.h"
 
+#include "pdlp/cupdlp/cupdlp_utils.h"
 #include "pdlp/hipdlp/logger.hpp"
 #include "pdlp/hipdlp/pdhg.hpp"
 #include "pdlp/hipdlp/restart.hpp"
-#include "pdlp/cupdlp/cupdlp_utils.h"
 
 HighsStatus solveLpHiPdlp(HighsLpSolverObject& solver_object) {
   return solveLpHiPdlp(solver_object.options_, solver_object.timer_,
@@ -68,7 +68,9 @@ HighsStatus solveLpHiPdlp(const HighsOptions& options, HighsTimer& timer,
   }
 
   // print col_value and row_dual
-  debugPdlpFinalSolutionLog(pdlp.debug_pdlp_log_file_, pdlp_solution.col_value.data(), lp.num_col_, pdlp_solution.row_dual.data(), lp.num_row_);
+  debugPdlpFinalSolutionLog(pdlp.debug_pdlp_log_file_,
+                            pdlp_solution.col_value.data(), lp.num_col_,
+                            pdlp_solution.row_dual.data(), lp.num_row_);
   pdlp.solveReturn();
   // --- Print Summary ---
   pdlp.logSummary();
@@ -118,15 +120,15 @@ HighsStatus solveLpHiPdlp(const HighsOptions& options, HighsTimer& timer,
   }
   assert(termination_status == TerminationStatus::OPTIMAL ||
          termination_status == TerminationStatus::TIMEOUT);
-  //highs_solution.col_value = x;
-  //highs_solution.col_value.resize(lp.num_col_);
-  //highs_solution.row_dual = y;
-  //lp.a_matrix_.product(highs_solution.row_value, highs_solution.col_value);
-  //lp.a_matrix_.productTranspose(highs_solution.col_dual,
-  //                              highs_solution.row_dual);
-  //for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
-  //  highs_solution.col_dual[iCol] =
-  //      lp.col_cost_[iCol] - highs_solution.col_dual[iCol];
+  // highs_solution.col_value = x;
+  // highs_solution.col_value.resize(lp.num_col_);
+  // highs_solution.row_dual = y;
+  // lp.a_matrix_.product(highs_solution.row_value, highs_solution.col_value);
+  // lp.a_matrix_.productTranspose(highs_solution.col_dual,
+  //                               highs_solution.row_dual);
+  // for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+  //   highs_solution.col_dual[iCol] =
+  //       lp.col_cost_[iCol] - highs_solution.col_dual[iCol];
   highs_solution.value_valid = true;
   highs_solution.dual_valid = true;
   return HighsStatus::kOk;

highs/pdlp/cupdlp/cupdlp_solver.c

@@ -847,6 +847,7 @@ void PDHG_Print_Iter(CUPDLPwork *pdhg) {
   const double gap_log_value = resobj->dRelObjGap;
   const double primal_feas_log_value = resobj->dPrimalFeas / (1.0 + scaling->dNormRhs);
   const double dual_feas_log_value = resobj->dDualFeas / (1.0 + scaling->dNormCost);
+
   cupdlp_printf("%9d  %+15.8e  %+15.8e  %+8.2e  %10.2e  %8.2e %7s [L]\n",
                 timers->nIter, resobj->dPrimalObj, resobj->dDualObj,
 		//                resobj->dDualityGap, resobj->dPrimalFeas, resobj->dDualFeas,
@@ -863,7 +864,6 @@ void PDHG_Print_Iter_Average(CUPDLPwork *pdhg) {
     cupdlp_snprintf(timeString, 8, "%6.2fs", timers->dSolvingTime);
   else
     cupdlp_snprintf(timeString, 8, "%6ds", (cupdlp_int)timers->dSolvingTime);
-
   CUPDLPscaling *scaling = pdhg->scaling;
   const double gap_log_value = resobj->dRelObjGapAverage;
   const double primal_feas_log_value = resobj->dPrimalFeasAverage / (1.0 + scaling->dNormRhs);
@@ -966,6 +966,22 @@ cupdlp_retcode PDHG_Solve(const cupdlp_int* has_variables, CUPDLPwork *pdhg) {
 	iter_log_since_header++;
       }
 
+      debugPdlpFeasOptLog(pdhg->debug_pdlp_log_file_,
+			  pdhg->timers->nIter, 
+			  pdhg->resobj->dPrimalObj,
+			  pdhg->resobj->dDualObj,
+			  pdhg->resobj->dRelObjGap,
+			  pdhg->resobj->dPrimalFeas / (1.0 + pdhg->scaling->dNormRhs),
+			  pdhg->resobj->dDualFeas / (1.0 + pdhg->scaling->dNormCost), "[L]");
+      debugPdlpFeasOptLog(pdhg->debug_pdlp_log_file_,
+			  pdhg->timers->nIter, 
+			  pdhg->resobj->dPrimalObjAverage,
+			  pdhg->resobj->dDualObjAverage,
+			  pdhg->resobj->dRelObjGapAverage,
+			  pdhg->resobj->dPrimalFeasAverage / (1.0 + pdhg->scaling->dNormRhs),
+			  pdhg->resobj->dDualFeasAverage / (1.0 + pdhg->scaling->dNormCost), "[A]");
+      debugPdlpIterHeaderLog(pdhg->debug_pdlp_log_file_);
+      
       // Termination check printing is only done when printing is full
       int termination_print = bool_print && full_print;
       if (PDHG_Check_Termination(pdhg, termination_print)) {

highs/pdlp/cupdlp/cupdlp_step.c

@@ -19,9 +19,9 @@ void PDHG_primalGradientStep(CUPDLPwork *work, CUPDLPvec *xUpdate,
   CUPDLPproblem *problem = work->problem;
 
 // print norm of A'y
-double aty_norm = 0.0;
-cupdlp_twoNorm(work, problem->nCols, ATy->data, &aty_norm);
-  work->debug_pdlp_data_.aty_norm = aty_norm;
+double debug_pdlp_data_aty_norm = 0.0;
+cupdlp_twoNorm(work, problem->nCols, ATy->data, &debug_pdlp_data_aty_norm);
+  work->debug_pdlp_data_.aty_norm = debug_pdlp_data_aty_norm;
 
 #if !defined(CUPDLP_CPU) && USE_KERNELS
   cupdlp_pgrad_cuda(xUpdate->data, x->data, problem->cost,
@@ -197,10 +197,10 @@ void PDHG_Update_Iterate_Constant_Step_Size(CUPDLPwork *pdhg) {
 
   Ax(pdhg, ax, x);
   //pint norm of Ax
-  double ax_norm = 0.0;
+  double debug_pdlp_data_ax_norm = 0.0;
   cupdlp_twoNorm(pdhg, problem->nRows, ax->data,
-                  &ax_norm);
-  pdhg->debug_pdlp_data_.ax_norm = ax_norm;
+                  &debug_pdlp_data_ax_norm);
+  pdhg->debug_pdlp_data_.ax_norm = debug_pdlp_data_ax_norm;
   ATy(pdhg, aty, y);
 
   // x^{k+1} = proj_{X}(x^k - dPrimalStep * (c - A'y^k))
@@ -410,9 +410,16 @@ void PDHG_Compute_Average_Iterate(CUPDLPwork *work) {
   Ax(work, iterates->axAverage, iterates->xAverage);
   ATy(work, iterates->atyAverage, iterates->yAverage);
 
-  cupdlp_float ax_norm = 0.0;
-  cupdlp_twoNormSquared(work, lp->nCols, iterates->axAverage->data, &ax_norm);
-  work->debug_pdlp_data_.ax_norm = ax_norm;
+  cupdlp_float debug_pdlp_data_ax_average_norm = 0.0;
+  cupdlp_twoNormSquared(work, lp->nCols, iterates->axAverage->data, &debug_pdlp_data_ax_average_norm);
+  work->debug_pdlp_data_.ax_average_norm = debug_pdlp_data_ax_average_norm;
+
+  /*
+  // Uncomment this once Yanyu is computing aty_average_norm 
+  cupdlp_float debug_pdlp_data_aty_average_norm = 0.0;
+  cupdlp_twoNormSquared(work, lp->nCols, iterates->atyAverage->data, &debug_pdlp_data_aty_average_norm);
+  work->debug_pdlp_data_.aty_average_norm = debug_pdlp_data_aty_average_norm;
+  */
 }
 
 void PDHG_Update_Average(CUPDLPwork *work) {

highs/pdlp/cupdlp/cupdlp_utils.c

@@ -1803,7 +1803,19 @@ void debugPdlpIterLog(FILE* file, const int iter_num, const struct DebugPdlpData
 	  beta);
 }
 
-void debugPdlpFeasOptLog(FILE* file) {
+void debugPdlpFeasOptLog(FILE* file, 
+			 const int iter_num,
+			 const double primal_obj, const double dual_obj,
+			 const double gap, const double primal_feas, const double dual_feas,
+			 const char* type) {
+  fprintf(file,
+	  "%6d Feasibility-optimality %s\n"
+	  "  primal_obj  = %11.4g\n"
+	  "  dual_obj    = %11.4g\n"
+	  "  gap         = %11.4g\n"
+	  "  primal_feas = %11.4g\n"
+	  "  dual_feas   = %11.4g\n",
+	  iter_num, type, primal_obj, dual_obj, gap, primal_feas, dual_feas);
 }
 
 void debugPdlpRestartLog(FILE* file, const int iter_num, const double current_score, const double average_score) {

highs/pdlp/cupdlp/cupdlp_utils.h

@@ -196,7 +196,11 @@ void writeSol(const char *fout, cupdlp_int nCols, cupdlp_int nRows,
   void debugPdlpIterHeaderLog(FILE* file);
   void debugPdlpDataInitialise(struct DebugPdlpData* debug_pdlp);
   void debugPdlpIterLog(FILE* file, const int iter_num, const struct DebugPdlpData* debug_pdlp, const double beta);
-  void debugPdlpFeasOptLog(FILE* file);
+  void debugPdlpFeasOptLog(FILE* file,
+			   const int iter_num,
+			   const double primal_obj, const double dual_obj,
+			   const double gap, const double primal_feas, const double dual_feas,
+			   const char* type);
   void debugPdlpRestartLog(FILE* file, const int iter_num, const double current_score, const double average_score);
   void debugPdlpRestarScoretLog(FILE* file, const double weight_squared, const double primal_feas, const double dual_feas, const double obj_gap);
   void debugPdlpFinalSolutionLog(FILE* file, const double* x, int nCols, const double* y, int nRows);

highs/pdlp/hipdlp/logger.cc

@@ -63,7 +63,6 @@ void Logger::print_header() {
   info("------------------------------------------------------------");
 }
 
-
 void Logger::print_params(const PrimalDualParams& params) {
   info("\nSolver Parameters:");
   std::stringstream ss;