Fix 1966a

check/TestIpm.cpp

@@ -81,3 +81,56 @@ TEST_CASE("test-analytic-centre-box", "[highs_ipm]") {
   REQUIRE(solution_norm < 1e-6);
   if (dev_run) printf("Analytic centre solution norm is %g\n", solution_norm);
 }
+
+TEST_CASE("test-1966", "[highs_ipm]") {
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  const HighsInfo& info = highs.getInfo();
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 2;
+  lp.col_cost_ = {2, -1};
+  lp.col_lower_ = {0, 0};
+  lp.col_upper_ = {kHighsInf, kHighsInf};
+  lp.row_lower_ = {-kHighsInf, 2};
+  lp.row_upper_ = {1, kHighsInf};
+  lp.a_matrix_.start_ = {0, 2, 4};
+  lp.a_matrix_.index_ = {0, 1, 0, 1};
+  lp.a_matrix_.value_ = {1, -1, 1, -1};
+  highs.passModel(lp);
+  highs.setOptionValue("solver", kIpmString);
+  highs.setOptionValue("presolve", kHighsOffString);
+
+  if (dev_run) printf("\nWith default residual tolerances\n");
+  highs.run();
+  if (dev_run) {
+    highs.writeSolution("", kSolutionStylePretty);
+    printf("Num primal infeasibilities = %d\n",
+           int(info.num_primal_infeasibilities));
+    printf("Max primal infeasibility   = %g\n", info.max_primal_infeasibility);
+    printf("Sum primal infeasibilities = %g\n",
+           info.sum_primal_infeasibilities);
+    printf("Num   dual infeasibilities = %d\n",
+           int(info.num_dual_infeasibilities));
+    printf("Max   dual infeasibility   = %g\n", info.max_dual_infeasibility);
+    printf("Sum   dual infeasibilities = %g\n", info.sum_dual_infeasibilities);
+  }
+  highs.clearSolver();
+
+  if (dev_run) printf("\nWith infinite residual tolerances\n");
+  highs.setOptionValue("primal_residual_tolerance", 1e30);
+  highs.setOptionValue("dual_residual_tolerance", 1e30);
+  highs.run();
+  if (dev_run) {
+    highs.writeSolution("", kSolutionStylePretty);
+    printf("Num primal infeasibilities = %d\n",
+           int(info.num_primal_infeasibilities));
+    printf("Max primal infeasibility   = %g\n", info.max_primal_infeasibility);
+    printf("Sum primal infeasibilities = %g\n",
+           info.sum_primal_infeasibilities);
+    printf("Num   dual infeasibilities = %d\n",
+           int(info.num_dual_infeasibilities));
+    printf("Max   dual infeasibility   = %g\n", info.max_dual_infeasibility);
+    printf("Sum   dual infeasibilities = %g\n", info.sum_dual_infeasibilities);
+  }
+}

check/TestLpSolvers.cpp

@@ -438,3 +438,35 @@ TEST_CASE("dual-objective-upper-bound", "[highs_lp_solver]") {
   if (dev_run) printf("\nOptimal objective value error = %g\n", error);
   REQUIRE(error < 1e-10);
 }
+
+TEST_CASE("blending-lp-ipm", "[highs_lp_solver]") {
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 2;
+  lp.col_cost_ = {-8, -10};
+  lp.col_lower_ = {0, 0};
+  lp.col_upper_ = {kHighsInf, kHighsInf};
+  lp.row_lower_ = {-kHighsInf, -kHighsInf};
+  lp.row_upper_ = {80, 120};
+  lp.a_matrix_.start_ = {0, 2, 4};
+  lp.a_matrix_.index_ = {0, 1, 0, 1};
+  lp.a_matrix_.value_ = {1, 1, 2, 4};
+  highs.passModel(lp);
+  highs.setOptionValue("solver", kIpmString);
+  highs.setOptionValue("presolve", kHighsOffString);
+  highs.run();
+  HighsInfo info = highs.getInfo();
+  if (dev_run) {
+    printf("Num primal infeasibilities = %d\n",
+           int(info.num_primal_infeasibilities));
+    printf("Max primal infeasibility   = %g\n", info.max_primal_infeasibility);
+    printf("Sum primal infeasibilities = %g\n",
+           info.sum_primal_infeasibilities);
+    printf("Num   dual infeasibilities = %d\n",
+           int(info.num_dual_infeasibilities));
+    printf("Max   dual infeasibility   = %g\n", info.max_dual_infeasibility);
+    printf("Sum   dual infeasibilities = %g\n", info.sum_dual_infeasibilities);
+  }
+}

check/TestNames.cpp

@@ -120,13 +120,12 @@ TEST_CASE("highs-names", "[highs_names]") {
 }
 
 TEST_CASE("highs-model-name", "[model_names]") {
-
   Highs highs;
   const HighsLp& lp = highs.getLp();
 
   std::string name = lp.model_name_;
   REQUIRE(name == "");
-  
+
   highs.passModelName("new_name");
   name = lp.model_name_;
   REQUIRE(name == "new_name");

src/lp_data/Highs.cpp

@@ -1959,10 +1959,6 @@ HighsStatus Highs::getReducedColumn(const HighsInt col, double* col_vector,
 
 HighsStatus Highs::getKappa(double& kappa, const bool exact,
                             const bool report) {
-  printf(
-      "Highs::getKappa basis_.valid = %d, ekk_instance_.status_.has_invert = "
-      "%d\n",
-      int(basis_.valid), int(ekk_instance_.status_.has_invert));
   if (!ekk_instance_.status_.has_invert)
     return invertRequirementError("getBasisInverseRow");
   kappa = ekk_instance_.computeBasisCondition(this->model_.lp_, exact, report);

src/lp_data/HighsOptions.h

@@ -389,6 +389,8 @@ struct HighsOptionsStruct {
   bool less_infeasible_DSE_choose_row;
   bool use_original_HFactor_logic;
   bool run_centring;
+  double primal_residual_tolerance;
+  double dual_residual_tolerance;
   HighsInt max_centring_steps;
   double centring_ratio_tolerance;
 
@@ -524,6 +526,8 @@ struct HighsOptionsStruct {
         run_centring(false),
         max_centring_steps(0),
         centring_ratio_tolerance(0.0),
+        primal_residual_tolerance(0.0),
+        dual_residual_tolerance(0.0),
         icrash(false),
         icrash_dualize(false),
         icrash_strategy(""),
@@ -1378,6 +1382,16 @@ class HighsOptions : public HighsOptionsStruct {
         advanced, &centring_ratio_tolerance, 0, 100, kHighsInf);
     records.push_back(record_double);
 
+    record_double = new OptionRecordDouble(
+        "primal_residual_tolerance", "Primal residual tolerance", advanced,
+        &primal_residual_tolerance, 1e-10, 1e-7, kHighsInf);
+    records.push_back(record_double);
+
+    record_double = new OptionRecordDouble(
+        "dual_residual_tolerance", "Dual residual tolerance", advanced,
+        &dual_residual_tolerance, 1e-10, 1e-7, kHighsInf);
+    records.push_back(record_double);
+
     // Set up the log_options aliases
     log_options.clear();
     log_options.log_stream =

src/lp_data/HighsSolution.cpp

@@ -282,7 +282,8 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp,
         relative_primal_infeasibility, dual_infeasibility, value_residual);
     if (!status_value_ok)
       highsLogUser(options.log_options, HighsLogType::kError,
-                   "getKktFailures: %s %d status-value error: [%g; %g; %g] has "
+                   "getKktFailures: %s %d status-value error: [%23.18g; "
+                   "%23.18g; %23.18g] has "
                    "residual %g\n",
                    iVar < lp.num_col_ ? "Column" : "Row   ",
                    iVar < lp.num_col_ ? int(iVar) : int(iVar - lp.num_col_),
@@ -531,14 +532,20 @@ bool getVariableKktFailures(const double primal_feasibility_tolerance,
     // If the variable is basic, then consider it not to be at a bound
     // so that any dual value yields an infeasibility value
     if (*status_pointer == HighsBasisStatus::kBasic) at_a_bound = false;
+    // With very large values, accuracy is lost in adding/subtracting
+    // the feasibility tolerance from the bounds, so skip if this may
+    // occur
+    //
     // Check that kLower and kUpper are consistent with value and
     // bounds - for debugging QP basis errors
     if (*status_pointer == HighsBasisStatus::kLower) {
-      status_value_ok = value >= lower - primal_feasibility_tolerance &&
-                        value <= lower + primal_feasibility_tolerance;
+      if (std::fabs(lower) / primal_feasibility_tolerance < 1e-16)
+        status_value_ok = value >= lower - primal_feasibility_tolerance &&
+                          value <= lower + primal_feasibility_tolerance;
     } else if (*status_pointer == HighsBasisStatus::kUpper) {
-      status_value_ok = value >= upper - primal_feasibility_tolerance &&
-                        value <= upper + primal_feasibility_tolerance;
+      if (std::fabs(upper) / primal_feasibility_tolerance < 1e-16)
+        status_value_ok = value >= upper - primal_feasibility_tolerance &&
+                          value <= upper + primal_feasibility_tolerance;
     }
   }
   if (at_a_bound) {
@@ -1369,6 +1376,73 @@ void accommodateAlienBasis(HighsLpSolverObject& solver_object) {
   assert(num_basic_variables == num_row);
 }
 
+void correctResiduals(HighsLpSolverObject& solver_object) {
+  HighsOptions& options = solver_object.options_;
+  HighsSolution& solution = solver_object.solution_;
+  const bool& have_primal_solution = solution.value_valid;
+  const bool& have_dual_solution = solution.dual_valid;
+  assert(have_primal_solution);
+  assert(have_dual_solution);
+  HighsLp& lp = solver_object.lp_;
+  std::vector<double> check_row_value;
+  std::vector<double> check_col_dual;
+  lp.a_matrix_.productQuad(check_row_value, solution.col_value);
+  if (have_dual_solution) {
+    lp.a_matrix_.productTransposeQuad(check_col_dual, solution.row_dual);
+    for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+      check_col_dual[iCol] -= lp.col_cost_[iCol];
+  }
+  double norm_primal_residual = 0;
+  HighsInt num_primal_correction = 0;
+  double max_primal_correction = 0;
+  double sum_primal_correction = 0;
+  double use_primal_residual_tolerance = options.primal_residual_tolerance;
+
+  double norm_dual_residual = 0;
+  HighsInt num_dual_correction = 0;
+  double max_dual_correction = 0;
+  double sum_dual_correction = 0;
+  double use_dual_residual_tolerance = options.dual_residual_tolerance;
+
+  for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) {
+    double primal_residual = check_row_value[iRow] - solution.row_value[iRow];
+    double abs_primal_residual = std::fabs(primal_residual);
+    if (abs_primal_residual > use_primal_residual_tolerance) {
+      solution.row_value[iRow] += primal_residual;
+      double check_primal_residual =
+          std::fabs(check_row_value[iRow] - solution.row_value[iRow]);
+      assert(check_primal_residual < 1e-10);
+      num_primal_correction++;
+      max_primal_correction =
+          std::max(abs_primal_residual, max_primal_correction);
+      sum_primal_correction += abs_primal_residual;
+    }
+    norm_primal_residual = std::max(abs_primal_residual, norm_primal_residual);
+  }
+  if (have_dual_solution) {
+    for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) {
+      double dual_residual = check_col_dual[iCol] + solution.col_dual[iCol];
+      double abs_dual_residual = std::fabs(dual_residual);
+      if (abs_dual_residual > use_dual_residual_tolerance) {
+        solution.col_dual[iCol] -= dual_residual;
+        num_dual_correction++;
+        max_dual_correction = std::max(abs_dual_residual, max_dual_correction);
+        sum_dual_correction += abs_dual_residual;
+      }
+      norm_dual_residual = std::max(abs_dual_residual, norm_dual_residual);
+    }
+  }
+
+  if (num_primal_correction > 0 || num_dual_correction > 0)
+    highsLogUser(
+        options.log_options, HighsLogType::kWarning,
+        "LP solver residuals: primal = %g; dual = %g yield num/max/sum primal "
+        "(%d/%g/%g) and dual (%d/%g/%g) corrections\n",
+        norm_primal_residual, norm_dual_residual, int(num_primal_correction),
+        max_primal_correction, sum_primal_correction, int(num_dual_correction),
+        max_dual_correction, sum_dual_correction);
+}
+
 void resetModelStatusAndHighsInfo(HighsLpSolverObject& solver_object) {
   solver_object.model_status_ = HighsModelStatus::kNotset;
   solver_object.highs_info_.objective_function_value = 0;

src/lp_data/HighsSolution.h

@@ -123,6 +123,8 @@ HighsStatus formSimplexLpBasisAndFactor(
 
 void accommodateAlienBasis(HighsLpSolverObject& solver_object);
 
+void correctResiduals(HighsLpSolverObject& solver_object);
+
 void resetModelStatusAndHighsInfo(HighsLpSolverObject& solver_object);
 void resetModelStatusAndHighsInfo(HighsModelStatus& model_status,
                                   HighsInfo& highs_info);

src/lp_data/HighsSolve.cpp

@@ -72,13 +72,22 @@ HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) {
                                           return_status, "solveLpCupdlp");
     }
     if (return_status == HighsStatus::kError) return return_status;
+    if (solver_object.model_status_ == HighsModelStatus::kOptimal)
+      // IPM (and PDLP?) can claim optimality with large primal and/or
+      // dual residual errors, so correct any residual errors that
+      // exceed the tolerance
+      correctResiduals(solver_object);
     // Non-error return requires a primal solution
     assert(solver_object.solution_.value_valid);
     // Get the objective and any KKT failures
     solver_object.highs_info_.objective_function_value =
         solver_object.lp_.objectiveValue(solver_object.solution_.col_value);
     getLpKktFailures(options, solver_object.lp_, solver_object.solution_,
                      solver_object.basis_, solver_object.highs_info_);
+    if (solver_object.model_status_ == HighsModelStatus::kOptimal &&
+        (solver_object.highs_info_.num_primal_infeasibilities > 0 ||
+         solver_object.highs_info_.num_dual_infeasibilities))
+      solver_object.model_status_ = HighsModelStatus::kUnknown;
     if (options.solver == kIpmString || options.run_centring) {
       // Setting the IPM-specific values of (highs_)info_ has been done in
       // solveLpIpx

src/pdlp/CupdlpWrapper.cpp

@@ -207,7 +207,7 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer,
 #if CUPDLP_DEBUG
   analysePdlpSolution(options, lp, highs_solution);
 #endif
-  
+
   free(cost);
   free(lower);
   free(upper);
@@ -225,15 +225,15 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer,
   free(prob->lower);
   free(prob->upper);
   free(prob->rhs);
-   
+
   free(prob->hasLower);
   free(prob->hasUpper);
 
   free(prob->data->csr_matrix->rowMatBeg);
   free(prob->data->csr_matrix->rowMatIdx);
   free(prob->data->csr_matrix->rowMatElem);
   free(prob->data->csr_matrix);
-   
+
   free(prob->data->csc_matrix->colMatBeg);
   free(prob->data->csc_matrix->colMatIdx);
   free(prob->data->csc_matrix->colMatElem);
@@ -246,13 +246,11 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer,
   free(csc_cpu->colMatBeg);
   free(csc_cpu->colMatIdx);
   free(csc_cpu->colMatElem);
-  
+
   free(csc_cpu);
 
-  if (scaling->rowScale != nullptr)
-    free(scaling->rowScale);
-  if (scaling->colScale != nullptr)
-    free(scaling->colScale);
+  if (scaling->rowScale != nullptr) free(scaling->rowScale);
+  if (scaling->colScale != nullptr) free(scaling->colScale);
   free(scaling);
 
   return HighsStatus::kOk;

src/util/HighsSparseMatrix.cpp

@@ -1174,6 +1174,32 @@ void HighsSparseMatrix::productQuad(vector<double>& result,
   }
 }
 
+void HighsSparseMatrix::productTransposeQuad(
+    vector<double>& result, const vector<double>& row,
+    const HighsInt debug_report) const {
+  assert(this->formatOk());
+  assert((int)row.size() >= this->num_row_);
+  result.assign(this->num_col_, 0.0);
+  if (this->isColwise()) {
+    for (HighsInt iCol = 0; iCol < this->num_col_; iCol++) {
+      HighsCDouble value = 0.0;
+      for (HighsInt iEl = this->start_[iCol]; iEl < this->start_[iCol + 1];
+           iEl++)
+        value += row[this->index_[iEl]] * this->value_[iEl];
+      result[iCol] = double(value);
+    }
+  } else {
+    std::vector<HighsCDouble> value(this->num_col_, 0);
+    for (HighsInt iRow = 0; iRow < this->num_row_; iRow++) {
+      for (HighsInt iEl = this->start_[iRow]; iEl < this->start_[iRow + 1];
+           iEl++)
+        value[this->index_[iEl]] += row[iRow] * this->value_[iEl];
+    }
+    for (HighsInt iCol = 0; iCol < this->num_col_; iCol++)
+      result[iCol] = double(value[iCol]);
+  }
+}
+
 void HighsSparseMatrix::productTransposeQuad(
     vector<double>& result_value, vector<HighsInt>& result_index,
     const HVector& column, const HighsInt debug_report) const {

src/util/HighsSparseMatrix.h

@@ -94,6 +94,9 @@ class HighsSparseMatrix {
   void productTranspose(vector<double>& result, const vector<double>& x) const;
   void productQuad(vector<double>& result, const vector<double>& x,
                    const HighsInt debug_report = kDebugReportOff) const;
+  void productTransposeQuad(
+      vector<double>& result_value, const vector<double>& x,
+      const HighsInt debug_report = kDebugReportOff) const;
   void productTransposeQuad(
       vector<double>& result_value, vector<HighsInt>& result_index,
       const HVector& x, const HighsInt debug_report = kDebugReportOff) const;