Merge pull request ERGO-Code#2321 from fwesselm/coeffStren

Fix 1517

Author: jajhall

highs/presolve/HPresolve.cpp

@@ -3774,28 +3774,41 @@ HPresolve::Result HPresolve::rowPresolve(HighsPostsolveStack& postsolve_stack,
       }
 
       auto strengthenCoefs = [&](HighsCDouble& rhs, HighsInt direction,
-                                 double maxCoefValue) {
+                                 HighsCDouble maxAbsCoefValue) {
         // iterate over non-zero positions instead of iterating over the
         // HighsMatrixSlice (provided by HPresolve::getStoredRow) because the
         // latter contains pointers to Acol and Avalue that may be invalidated
         // if these vectors are reallocated (see std::vector::push_back
         // performed in HPresolve::addToMatrix).
         for (HighsInt rowiter : rowpositions) {
+          // max. absolute coefficient should not be negative
+          assert(maxAbsCoefValue >= 0);
+
           // get column index and coefficient
           HighsInt col = Acol[rowiter];
-          double val = Avalue[rowiter];
+          double val = direction * Avalue[rowiter];
+
+          // get lower and upper bounds
+          double col_lower = impliedRowBounds.getImplVarLower(row, col);
+          double col_upper = impliedRowBounds.getImplVarUpper(row, col);
 
           // skip continuous variables
           if (model->integrality_[col] == HighsVarType::kContinuous) continue;
 
-          if (direction * val > maxCoefValue + primal_feastol) {
-            double delta = direction * maxCoefValue - val;
-            addToMatrix(row, col, delta);
-            rhs += delta * model->col_upper_[col];
-          } else if (direction * val < -maxCoefValue - primal_feastol) {
-            double delta = -direction * maxCoefValue - val;
-            addToMatrix(row, col, delta);
-            rhs += delta * model->col_lower_[col];
+          if (val > maxAbsCoefValue + primal_feastol) {
+            assert(col_upper != kHighsInf);
+            // new matrix coefficient is direction * maxAbsCoefValue; subtract
+            // existing matrix coefficient to get delta
+            HighsCDouble delta = direction * (maxAbsCoefValue - val);
+            addToMatrix(row, col, static_cast<double>(delta));
+            rhs += delta * col_upper;
+          } else if (val < -maxAbsCoefValue - primal_feastol) {
+            assert(col_lower != -kHighsInf);
+            // new matrix coefficient is (-direction) * maxAbsCoefValue;
+            // subtract existing matrix coefficient to get delta
+            HighsCDouble delta = -direction * (maxAbsCoefValue + val);
+            addToMatrix(row, col, static_cast<double>(delta));
+            rhs += delta * col_lower;
           }
         }
       };
@@ -3805,7 +3818,8 @@ HPresolve::Result HPresolve::rowPresolve(HighsPostsolveStack& postsolve_stack,
         // <= constraint: try to strengthen coefficients
         HighsCDouble rhs = model->row_upper_[row];
         strengthenCoefs(rhs, HighsInt{1},
-                        impliedRowUpper - model->row_upper_[row]);
+                        static_cast<HighsCDouble>(impliedRowUpper) -
+                            model->row_upper_[row]);
         model->row_upper_[row] = static_cast<double>(rhs);
       }
 
@@ -3814,7 +3828,8 @@ HPresolve::Result HPresolve::rowPresolve(HighsPostsolveStack& postsolve_stack,
         // >= constraint: try to strengthen coefficients
         HighsCDouble rhs = model->row_lower_[row];
         strengthenCoefs(rhs, HighsInt{-1},
-                        model->row_lower_[row] - impliedRowLower);
+                        model->row_lower_[row] -
+                            static_cast<HighsCDouble>(impliedRowLower));
         model->row_lower_[row] = static_cast<double>(rhs);
       }
     }