Merge pull request #2410 from fwesselm/fix-2409

Fix 2409

Author: jajhall

check/TestPresolve.cpp

@@ -724,3 +724,14 @@ TEST_CASE("presolve-issue-2402", "[highs_test_presolve]") {
   REQUIRE(highs.presolve() == HighsStatus::kOk);
   REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kInfeasible);
 }
+
+TEST_CASE("presolve-issue-2409", "[highs_test_presolve]") {
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/issue-2409.mps";
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.readModel(model_file);
+  REQUIRE(highs.presolve() == HighsStatus::kOk);
+  REQUIRE(highs.getModelPresolveStatus() ==
+          HighsPresolveStatus::kReducedToEmpty);
+}

check/instances/issue-2409.mps

@@ -0,0 +1,21 @@
+NAME        
+ROWS
+ N  Obj     
+ G  r0      
+ G  r1      
+COLUMNS
+    c0        Obj       -1
+    c0        r0        -1
+    c0        r1        1
+    MARK0000  'MARKER'                 'INTORG'
+    c1        Obj       1
+    c1        r0        1
+    c1        r1        1
+    MARK0001  'MARKER'                 'INTEND'
+RHS
+    RHS_V     r0        0.1
+    RHS_V     r1        0.1
+BOUNDS
+ FR BOUND     c0      
+ FR BOUND     c1      
+ENDATA

highs/presolve/HPresolve.cpp

@@ -4190,14 +4190,21 @@ HPresolve::Result HPresolve::detectDominatedCol(
         HPRESOLVE_CHECKED_CALL(removeRowSingletons(postsolve_stack));
       return checkLimits(postsolve_stack);
     } else if (analysis_.allow_rule_[kPresolveRuleForcingCol]) {
-      // get bound on dual (column) activity
-      HighsCDouble sum = 0;
-      if (direction > 0)
-        sum = impliedDualRowBounds.getSumUpperOrig(col);
-      else
-        sum = impliedDualRowBounds.getSumLowerOrig(col);
-      if (sum == 0.0) {
-        // remove column and rows
+      // get bound on column dual using original bounds on row duals
+      double boundOnColDual = direction > 0
+                                  ? -impliedDualRowBounds.getSumUpperOrig(
+                                        col, -model->col_cost_[col])
+                                  : -impliedDualRowBounds.getSumLowerOrig(
+                                        col, -model->col_cost_[col]);
+      if (boundOnColDual == 0.0) {
+        // 1. column's lower bound is infinite (i.e. column dual has upper bound
+        // of zero) and column dual's lower bound is zero as well
+        // (direction = 1) or
+        // 2. column's upper bound is infinite (i.e. column dual has lower bound
+        // of zero) and column dual's upper bound is zero as well
+        // (direction = -1).
+        // thus, the column dual is zero, and we can remove the column and
+        // all its rows
         if (logging_on) analysis_.startPresolveRuleLog(kPresolveRuleForcingCol);
         postsolve_stack.forcingColumn(
             col, getColumnVector(col), model->col_cost_[col], otherBound,

highs/presolve/HPresolve.h

@@ -390,7 +390,7 @@ class HPresolve {
   Result strengthenInequalities(HighsPostsolveStack& postsolve_stack,
                                 HighsInt& num_strenghtened);
 
-  HPresolve::Result detectImpliedIntegers();
+  Result detectImpliedIntegers();
 
   Result detectParallelRowsAndCols(HighsPostsolveStack& postsolve_stack);
 

highs/util/HighsLinearSumBounds.h

@@ -102,41 +102,36 @@ class HighsLinearSumBounds {
   double getResidualSumUpperOrig(HighsInt sum, HighsInt var,
                                  double coefficient) const;
 
-  double getSumLowerOrig(HighsInt sum) const {
-    return numInfSumLowerOrig[sum] == 0 ? double(sumLowerOrig[sum])
-                                        : -kHighsInf;
+  template <typename T = double>
+  double getSumLowerOrig(HighsInt sum, T offset = T()) const {
+    return numInfSumLowerOrig[sum] == 0
+               ? static_cast<double>(sumLowerOrig[sum] +
+                                     static_cast<HighsCDouble>(offset))
+               : -kHighsInf;
   }
 
-  double getSumUpperOrig(HighsInt sum) const {
-    return numInfSumUpperOrig[sum] == 0 ? double(sumUpperOrig[sum]) : kHighsInf;
+  template <typename T = double>
+  double getSumUpperOrig(HighsInt sum, T offset = T()) const {
+    return numInfSumUpperOrig[sum] == 0
+               ? static_cast<double>(sumUpperOrig[sum] +
+                                     static_cast<HighsCDouble>(offset))
+               : kHighsInf;
   }
 
-  double getSumLower(HighsInt sum) const {
-    return numInfSumLower[sum] == 0 ? double(sumLower[sum]) : -kHighsInf;
+  template <typename T = double>
+  double getSumLower(HighsInt sum, T offset = T()) const {
+    return numInfSumLower[sum] == 0
+               ? static_cast<double>(sumLower[sum] +
+                                     static_cast<HighsCDouble>(offset))
+               : -kHighsInf;
   }
 
-  double getSumUpper(HighsInt sum) const {
-    return numInfSumUpper[sum] == 0 ? double(sumUpper[sum]) : kHighsInf;
-  }
-
-  double getSumLower(HighsInt sum, double offset) const {
-    return numInfSumLower[sum] == 0 ? double(sumLower[sum] + offset)
-                                    : -kHighsInf;
-  }
-
-  double getSumUpper(HighsInt sum, double offset) const {
-    return numInfSumUpper[sum] == 0 ? double(sumUpper[sum] + offset)
-                                    : kHighsInf;
-  }
-
-  double getSumLower(HighsInt sum, HighsCDouble offset) const {
-    return numInfSumLower[sum] == 0 ? double(sumLower[sum] + offset)
-                                    : -kHighsInf;
-  }
-
-  double getSumUpper(HighsInt sum, HighsCDouble offset) const {
-    return numInfSumUpper[sum] == 0 ? double(sumUpper[sum] + offset)
-                                    : kHighsInf;
+  template <typename T = double>
+  double getSumUpper(HighsInt sum, T offset = T()) const {
+    return numInfSumUpper[sum] == 0
+               ? static_cast<double>(sumUpper[sum] +
+                                     static_cast<HighsCDouble>(offset))
+               : kHighsInf;
   }
 
   HighsInt getNumInfSumLower(HighsInt sum) const { return numInfSumLower[sum]; }