fix(presolve): clamp singleton row bound tolerance to machine epsilon

[HuAliAsu]

Problem

When presolving MILP models with large matrix coefficients, HiGHS can erroneously over-tighten integer variable bounds in singleton rows, cutting off valid optimal integer solutions.

Technical Root Cause

In HPresolve.cpp::removeRowSingletons (line 3273), the bound tolerance for singleton rows is calculated as:

const double boundTol = primal_feastol / std::max(1.0, std::fabs(val));

When a singleton row has a large matrix coefficient (val = 10^12), boundTol underflows to 1e-19, which is far below IEEE 754 machine epsilon (2.22e-16).

This causes over-rounding during integer bound tightening:

if (isIntegral) newColLower = std::ceil(newColLower - boundTol);

If floating-point division noise leaves newColLower at 3.0000000000000004 (instead of 3.0), subtracting 1e-19 has no effect because 1e-19 is smaller than machine epsilon. Consequently, std::ceil(3.0000000000000004) rounds up to 4.0 instead of 3.0, cutting off valid integer solutions.

The Fix

Clamp boundTol to a minimum threshold of 1e-15 (IEEE 754 machine epsilon safeguard):

const double boundTol = std::max(primal_feastol / std::max(1.0, std::fabs(val)), 1e-15);

[Opt-Mucca]

Similar to the other PR I think this change makes sense! @HuAliAsu can you chuck both into a single PR and make sure it respects clang formatting?

[jajhall]

Thanks for these @HuAliAsu

Although we'll accept clamping values to avoid underflow, based on the reproducible LP example, your problems are so badly scaled that HiGHS (and other solvers) are likely to have numerical issues (or fail) when trying to solve them. I've tried the reproducible LP example, and HiGHS (without presolve) finds the origin to be optimal - as does a commercial solver. Both solvers give strong warnings about the large values, and recommending that the model be reformulated.

Specifically, for the reproducible LP example,

• the costs should be scaled down by 1e6
• the equations should be scaled down so the largest coefficient of each is O(1)

Ideally all variables should have an optimal value of order 1 - or zero - and all binding constraints should have a RHS of order 1. Not only will this improve solver robustness, the solver is likely to run faster.

To learn more, look at this JuMP tutorial.