Add additional utilities to presolve

src/presolve/HPresolve.cpp

@@ -3089,6 +3089,20 @@ HPresolve::Result HPresolve::rowPresolve(HighsPostsolveStack& postsolve_stack,
               continue;
             }
 
+            auto remDoubletonEq = [&](HighsInt col, HighsInt binCol,
+                                      HighsInt direction) {
+              double bound = direction >= 0 ? model->col_upper_[col]
+                                            : model->col_lower_[col];
+              double scale =
+                  direction * (model->col_lower_[col] - model->col_upper_[col]);
+              double offset = bound - model->col_lower_[binCol] * scale;
+              postsolve_stack.doubletonEquation(
+                  -1, col, binCol, 1.0, -scale, offset, model->col_lower_[col],
+                  model->col_upper_[col], 0.0, false, false,
+                  HighsPostsolveStack::RowType::kEq, HighsEmptySlice());
+              substitute(col, binCol, offset, scale);
+            };
+
             if (std::signbit(binCoef) == std::signbit(nonz.value())) {
               // binary coefficient is positive:
               // setting the binary to its upper bound
@@ -3104,16 +3118,7 @@ HPresolve::Result HPresolve::rowPresolve(HighsPostsolveStack& postsolve_stack,
               // nonzCol = colUb - (colUb - colLb)(binCol - binLb)
               // nonzCol = colUb + binLb * (colUb - colLb) - (colUb - colLb) *
               // binCol
-              double scale = model->col_lower_[nonz.index()] -
-                             model->col_upper_[nonz.index()];
-              double offset = model->col_upper_[nonz.index()] -
-                              model->col_lower_[binCol] * scale;
-              postsolve_stack.doubletonEquation(
-                  -1, nonz.index(), binCol, 1.0, -scale, offset,
-                  model->col_lower_[nonz.index()],
-                  model->col_upper_[nonz.index()], 0.0, false, false,
-                  HighsPostsolveStack::RowType::kEq, HighsEmptySlice());
-              substitute(nonz.index(), binCol, offset, scale);
+              remDoubletonEq(nonz.index(), binCol, HighsInt{1});
             } else {
               // This case yields the following implications:
               // binCol = lb -> nonzCol = lb
@@ -3122,16 +3127,7 @@ HPresolve::Result HPresolve::rowPresolve(HighsPostsolveStack& postsolve_stack,
               // nonzCol = colLb + (colUb - colLb)(binCol - binLb)
               // nonzCol =
               //    colLb - binLb*(colUb - colLb) + (colUb - colLb)*binCol
-              double scale = model->col_upper_[nonz.index()] -
-                             model->col_lower_[nonz.index()];
-              double offset = model->col_lower_[nonz.index()] -
-                              model->col_lower_[binCol] * scale;
-              postsolve_stack.doubletonEquation(
-                  -1, nonz.index(), binCol, 1.0, -scale, offset,
-                  model->col_lower_[nonz.index()],
-                  model->col_upper_[nonz.index()], 0.0, false, false,
-                  HighsPostsolveStack::RowType::kEq, HighsEmptySlice());
-              substitute(nonz.index(), binCol, offset, scale);
+              remDoubletonEq(nonz.index(), binCol, HighsInt{-1});
             }
           }
 
@@ -5915,8 +5911,12 @@ HPresolve::Result HPresolve::detectParallelRowsAndCols(
     auto it = buckets.find(rowHashes[i]);
     decltype(it) last = it;
 
-    const HighsInt* numSingletonPtr = numRowSingletons.find(i);
-    HighsInt numSingleton = numSingletonPtr ? *numSingletonPtr : 0;
+    auto getNumSingletons = [&](HighsInt row) {
+      const HighsInt* numSingletonPtr = numRowSingletons.find(row);
+      return (numSingletonPtr ? *numSingletonPtr : 0);
+    };
+
+    const HighsInt numSingleton = getNumSingletons(i);
 
 #if !ENABLE_SPARSIFY_FOR_LP
     if (mipsolver == nullptr && options->lp_presolve_requires_basis_postsolve &&
@@ -5929,9 +5929,8 @@ HPresolve::Result HPresolve::detectParallelRowsAndCols(
       HighsInt parallelRowCand = it->second;
       last = it++;
 
-      numSingletonPtr = numRowSingletons.find(parallelRowCand);
-      const HighsInt numSingletonCandidate =
-          numSingletonPtr ? *numSingletonPtr : 0;
+      const HighsInt numSingletonCandidate = getNumSingletons(parallelRowCand);
+
 #if !ENABLE_SPARSIFY_FOR_LP
       if (mipsolver == nullptr &&
           options->lp_presolve_requires_basis_postsolve &&
@@ -6072,31 +6071,8 @@ HPresolve::Result HPresolve::detectParallelRowsAndCols(
         //    HIGHSINT_FORMAT ")\n", numSingleton, numSingletonCandidate,
         //    model->row_lower_[parallelRowCand] ==
         //        model->row_upper_[parallelRowCand]);
-        postsolve_stack.equalityRowAddition(parallelRowCand, i, -rowScale,
-                                            getStoredRow());
-        for (const HighsSliceNonzero& rowNz : getStoredRow()) {
-          HighsInt pos = findNonzero(parallelRowCand, rowNz.index());
-          if (pos != -1)
-            unlink(pos);  // all common nonzeros are cancelled, as the rows are
-                          // parallel
-          else            // might introduce a singleton
-            addToMatrix(parallelRowCand, rowNz.index(),
-                        -rowScale * rowNz.value());
-        }
-
-        if (model->row_upper_[parallelRowCand] != kHighsInf)
-          model->row_upper_[parallelRowCand] =
-              double(model->row_upper_[parallelRowCand] -
-                     HighsCDouble(rowScale) * model->row_upper_[i]);
-        if (model->row_lower_[parallelRowCand] != -kHighsInf)
-          model->row_lower_[parallelRowCand] =
-              double(model->row_lower_[parallelRowCand] -
-                     HighsCDouble(rowScale) * model->row_upper_[i]);
-
-        // parallelRowCand is now a singleton row, doubleton equation, or a row
-        // that contains only singletons and we let the normal row presolve
-        // handle the cases
-        HPRESOLVE_CHECKED_CALL(rowPresolve(postsolve_stack, parallelRowCand));
+        HPRESOLVE_CHECKED_CALL(equalityRowAddition(
+            postsolve_stack, i, parallelRowCand, -rowScale, getStoredRow()));
         delRow = parallelRowCand;
       } else if (model->row_lower_[parallelRowCand] ==
                  model->row_upper_[parallelRowCand]) {
@@ -6105,28 +6081,10 @@ HPresolve::Result HPresolve::detectParallelRowsAndCols(
         //    row and %" HIGHSINT_FORMAT " " "singletons in other inequality
         //    row\n", numSingletonCandidate, numSingleton);
         // the row parallelRowCand is an equation; add it to the other row
-        double scale = -rowMax[i].first / rowMax[parallelRowCand].first;
-        postsolve_stack.equalityRowAddition(i, parallelRowCand, scale,
-                                            getRowVector(parallelRowCand));
-        for (const HighsSliceNonzero& rowNz : getRowVector(parallelRowCand)) {
-          HighsInt pos = findNonzero(i, rowNz.index());
-          if (pos != -1)
-            unlink(pos);  // all common nonzeros are cancelled, as the rows are
-                          // parallel
-          else            // might introduce a singleton
-            addToMatrix(i, rowNz.index(), scale * rowNz.value());
-        }
-
-        if (model->row_upper_[i] != kHighsInf)
-          model->row_upper_[i] =
-              double(model->row_upper_[i] +
-                     HighsCDouble(scale) * model->row_upper_[parallelRowCand]);
-        if (model->row_lower_[i] != -kHighsInf)
-          model->row_lower_[i] =
-              double(model->row_lower_[i] +
-                     HighsCDouble(scale) * model->row_upper_[parallelRowCand]);
-
-        HPRESOLVE_CHECKED_CALL(rowPresolve(postsolve_stack, i));
+        HPRESOLVE_CHECKED_CALL(equalityRowAddition(
+            postsolve_stack, parallelRowCand, i,
+            -rowMax[i].first / rowMax[parallelRowCand].first,
+            getRowVector(parallelRowCand)));
         delRow = i;
       } else {
         assert(numSingleton == 1);
@@ -6170,6 +6128,36 @@ HPresolve::Result HPresolve::detectParallelRowsAndCols(
   return Result::kOk;
 }
 
+template <typename RowStorageFormat>
+HPresolve::Result HPresolve::equalityRowAddition(
+    HighsPostsolveStack& postsolve_stack, HighsInt stayrow, HighsInt removerow,
+    double scale, const HighsMatrixSlice<RowStorageFormat>& vector) {
+  postsolve_stack.equalityRowAddition(removerow, stayrow, scale, vector);
+  for (const auto& rowNz : vector) {
+    HighsInt pos = findNonzero(removerow, rowNz.index());
+    if (pos != -1)
+      unlink(pos);  // all common nonzeros are cancelled, as the rows are
+                    // parallel
+    else            // might introduce a singleton
+      addToMatrix(removerow, rowNz.index(), scale * rowNz.value());
+  }
+
+  if (model->row_upper_[removerow] != kHighsInf)
+    model->row_upper_[removerow] =
+        double(model->row_upper_[removerow] +
+               HighsCDouble(scale) * model->row_upper_[stayrow]);
+  if (model->row_lower_[removerow] != -kHighsInf)
+    model->row_lower_[removerow] =
+        double(model->row_lower_[removerow] +
+               HighsCDouble(scale) * model->row_upper_[stayrow]);
+
+  // row is now a singleton row, doubleton equation, or a row
+  // that contains only singletons and we let the normal row presolve
+  // handle the cases
+  HPRESOLVE_CHECKED_CALL(rowPresolve(postsolve_stack, removerow));
+  return Result::kOk;
+}
+
 void HPresolve::setRelaxedImpliedBounds() {
   double hugeBound = primal_feastol / kHighsTiny;
   for (HighsInt i = 0; i != model->num_col_; ++i) {

src/presolve/HPresolve.h

@@ -362,6 +362,11 @@ class HPresolve {
 
   Result detectParallelRowsAndCols(HighsPostsolveStack& postsolve_stack);
 
+  template <typename RowStorageFormat>
+  Result equalityRowAddition(HighsPostsolveStack& postsolve_stack,
+                             HighsInt stayrow, HighsInt removerow, double scale,
+                             const HighsMatrixSlice<RowStorageFormat>& vector);
+
   Result sparsify(HighsPostsolveStack& postsolve_stack);
 
   void setRelaxedImpliedBounds();