Added test code in 1.9.0

Author: jajhall

src/lp_data/Highs.cpp

@@ -4607,6 +4607,8 @@ void Highs::reportSolvedLpQpStats() {
     const double relative_primal_dual_gap =
         std::fabs(info_.objective_function_value - dual_objective_value) /
         std::max(1.0, std::fabs(info_.objective_function_value));
+    printf("P-D: primal objective = %g\n", info_.objective_function_value);
+    printf("P-D: dual   objective = %g\n", dual_objective_value);
     highsLogUser(log_options, HighsLogType::kInfo,
                  "Relative P-D gap    : %17.10e\n", relative_primal_dual_gap);
   }

src/lp_data/HighsSolution.cpp

@@ -625,6 +625,10 @@ bool computeDualObjectiveValue(const HighsLp& lp, const HighsSolution& solution,
 
   dual_objective_value = lp.offset_;
   double bound = 0;
+  assert(HighsInt(lp.col_lower_.size()) >= lp.num_col_);
+  assert(HighsInt(lp.col_upper_.size()) >= lp.num_col_);
+  assert(HighsInt(lp.row_lower_.size()) >= lp.num_row_);
+  assert(HighsInt(lp.row_upper_.size()) >= lp.num_row_);
   for (HighsInt iVar = 0; iVar < lp.num_col_ + lp.num_row_; iVar++) {
     const bool is_col = iVar < lp.num_col_;
     const HighsInt iRow = iVar - lp.num_col_;
@@ -634,6 +638,7 @@ bool computeDualObjectiveValue(const HighsLp& lp, const HighsSolution& solution,
         is_col ? solution.col_dual[iVar] : solution.row_dual[iRow];
     const double lower = is_col ? lp.col_lower_[iVar] : lp.row_lower_[iRow];
     const double upper = is_col ? lp.col_upper_[iVar] : lp.row_upper_[iRow];
+    if (!is_col) assert(iRow < lp.num_row_);
     if (lower <= -kHighsInf && upper >= kHighsInf) {
       // Free
       bound = 1;
@@ -642,6 +647,12 @@ bool computeDualObjectiveValue(const HighsLp& lp, const HighsSolution& solution,
       bound = primal < mid ? lower : upper;
     }
     dual_objective_value += bound * dual;
+    if (bound * dual)
+      printf(
+          "computeDualObjectiveValue: %s %6d:  bound = %11.5g; dual = %11.5g; "
+          "objective = %11.5g\n",
+          is_col ? "col" : "row", is_col ? int(iVar) : int(iRow), bound, dual,
+          dual_objective_value);
   }
   return true;
 }

src/simplex/HEkk.cpp

@@ -1141,6 +1141,13 @@ HighsStatus HEkk::solve(const bool force_phase2) {
   reportSimplexPhaseIterations(options_->log_options, iteration_count_, info_);
   if (return_status == HighsStatus::kError)
     return returnFromEkkSolve(return_status);
+  printf("EKK max_primal_infeasibility   = %g\n",
+         info_.max_primal_infeasibility);
+  printf("EKK sum_primal_infeasibilities = %g\n",
+         info_.sum_primal_infeasibilities);
+  printf("EKK max_  dual_infeasibility   = %g\n", info_.max_dual_infeasibility);
+  printf("EKK sum_  dual_infeasibilities = %g\n",
+         info_.sum_dual_infeasibilities);
   highsLogDev(options_->log_options, HighsLogType::kInfo,
               "EKK %s simplex solver returns %" HIGHSINT_FORMAT
               " primal and %" HIGHSINT_FORMAT