Fix #3557

check/TestPresolve.cpp

@@ -2034,6 +2034,43 @@ TEST_CASE("dominated-cols-container-overflow", "[highs_test_presolve]") {
   h.resetThreadScheduler(true);
 }
 
+TEST_CASE("doubleton-equations-reinserted-row", "[highs_test_presolve]") {
+  // model in issue 3557. with the rules for forcing rows and doubleton
+  // equations switched off, removeDoubletonEquations calls rowPresolve on the
+  // equation, which fixes a column and re-keys the equation in the set of
+  // equations, so the iterator was incremented after its entry was erased
+  //
+  // min -8 x1 + 10 x2
+  // s.t. 2 x0 - 5 x1 + 5 x2 = -16
+  //      -1 <= x0 <= 2, 0 <= x1 <= 3, -1 <= x2 <= 0, x0, x1, x2 integer
+  HighsLp lp;
+  lp.num_col_ = 3;
+  lp.num_row_ = 1;
+  lp.col_cost_ = {0, -8, 10};
+  lp.col_lower_ = {-1, 0, -1};
+  lp.col_upper_ = {2, 3, 0};
+  lp.row_lower_ = {-16};
+  lp.row_upper_ = {-16};
+  lp.a_matrix_.format_ = MatrixFormat::kColwise;
+  lp.a_matrix_.start_ = {0, 1, 2, 3};
+  lp.a_matrix_.index_ = {0, 0, 0};
+  lp.a_matrix_.value_ = {2, -5, 5};
+  lp.integrality_ = {HighsVarType::kInteger, HighsVarType::kInteger,
+                     HighsVarType::kInteger};
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  h.setOptionValue(
+      "presolve_rule_off",
+      (1 << kPresolveRuleForcingRow) | (1 << kPresolveRuleDoubletonEquation));
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+  REQUIRE(h.getInfo().objective_function_value == -34);
+
+  h.resetThreadScheduler(true);
+}
+
 #ifdef HIPO
 TEST_CASE("hipo-ipx-pdlp-no-crossover", "[highs_test_presolve]") {
   auto setPresolveRuleOff = [&]() {

highs/presolve/HPresolve.cpp

@@ -4907,11 +4907,11 @@ HPresolve::Result HPresolve::removeDoubletonEquations(
     assert(eq->first == rowsize[eqrow]);
     assert(isEquation(eqrow));
     if (rowsize[eqrow] > 2) return Result::kOk;
+    // rowPresolve may erase or re-key the entry of eqrow (e.g. removeFixedCol
+    // calls reinsertEquation), so continue from a copy of its key
+    const auto key = *eq;
     HPRESOLVE_CHECKED_CALL(rowPresolve(postsolve_stack, eqrow));
-    if (rowDeleted[eqrow])
-      eq = equations.begin();
-    else
-      ++eq;
+    eq = rowDeleted[eqrow] ? equations.begin() : equations.upper_bound(key);
   }
 
   return Result::kOk;