Drop @fastmath from L2_NORM/Linf_NORM so isfinite guards survive

[ChrisRackauckas-Claude]

Please ignore until reviewed by @ChrisRackauckas.

Released NonlinearSolve can return a successful retcode on a diverged solve. Reproduced against the registry, not this branch — NonlinearSolve v4.27.0 / NonlinearSolveBase v2.46.0:

using NonlinearSolve, SciMLBase
function nanres!(du, u, p)          # Inf - Inf = NaN, e.g. an Arrhenius overflow
    du[1] = exp(u[1]) - exp(2*u[1]) + 1.0
    du[2] = u[2] - 1.0
end
prob = NonlinearLeastSquaresProblem(
    NonlinearFunction(nanres!; resid_prototype = zeros(2)), [1000.0, 0.0])
sol = solve(prob, TrustRegion(); maxiters = 100)

sol.retcode                       # ReturnCode.StalledSuccess
SciMLBase.successful_retcode(sol) # true
sol.resid                         # [NaN, -1.0]

Cause

L2_NORM and Linf_NORM are defined with @fastmath (lib/NonlinearSolveBase/src/common_defaults.jl, the only @fastmath in the repo). @fastmath sqrt emits call fast double @llvm.sqrt.f64, and fast implies nnan ninf, so LLVM proves any isfinite on the result is true and deletes it. With the released code,

g(du::Vector{Float64}) = !isfinite(NonlinearSolveBase.L2_NORM(du))

compiles to, in full:

define i8 @julia_g(ptr noundef nonnull align 8 dereferenceable(24) %"du::Array") #0 {
top:
  ret i8 0
}

The safe-mode protective break is written exactly that way (termination_conditions.jl:256, if !isfinite(objective)), and the least-squares default internalnorm is Base.Fix2(norm, 2), which standardize_norms to L2_NORM. So for least-squares the guard is dead. Square NonlinearProblem defaults to Base.Fix1(maximum, abs) → Linf_NORM (maximum(abs, u), no @fastmath on the reduction) and is unaffected.

From there: the guard misses → TrustRegion rejects every step (ρ = NaN, so NaN > threshold is false) → u == uprev exactly → L2_NORM(u - uprev) = 0 ≤ abstol → the stall branch fires with leastsq == true → StalledSuccess, which successful_retcode reports as true.

Causal isolation — same mode, same input, only the norm differs:

internalnorm	du = [NaN, 1.0]	[-Inf, Inf]
Base.Fix2(norm, 2) → L2_NORM	Failure	Failure
Base.Fix1(maximum, abs) → Linf_NORM	Unstable	Unstable
the same L2_NORM behind @noinline	Unstable	Unstable

The @noinline row rules out an ordinary logic bug: the guard's input is correct (L2_NORM([NaN, 1.0]) is NaN at top level), and only inlining-plus-fold explains the miss. code_typed still contains the check; it is an LLVM fold, and it happens at -O1 and above on every Julia version tested (1.10, 1.11, 1.12, 1.13-rc2).

It is broader than NaN, and broader than isfinite

• A fully finite residual is enough. L2_NORM computes an unscaled sum(abs2), so any ‖F‖ ≳ 1.34e154 overflows to Inf internally. Starting from u0 = [200.0, 0.0] with r(u0) = [-5.22e173, -1.0] — every entry finite — master returns StalledSuccess/successful; this branch returns Unstable.
• isnan is folded too, though isinf is not.
• max/min swallow NaN, because Julia builds them on an internal isnan test which is also folded: on master max(L2_NORM([NaN,1.0]), 1.0) == 1.0. That launders NaN into the initial trust radius (trust_region.jl:336), utils.jl:313, homotopy_sweep.jl:839 and two SimpleNonlinearSolve sites. Removing the taint at source fixes all of them — verified NaN after the change.

Ordered comparisons (fcmp ole/ogt) are not affected and survive verbatim; every accept/reject site consuming a tainted norm was checked and behaves identically before and after, always in the reject direction.

Fix

Drop @fastmath from the two norms rather than rewriting the guards. Magnitude comparisons do not work here — NaN > 1e300 is false, so they catch Inf but not NaN — and !(x <= c) formulations are themselves fast-math-dependent. Screening the raw residual costs an extra O(n) pass per iteration.

The cost is nil, and the reason is that @fastmath was not doing the work. @simd already supplies the reassoc+contract that vectorizes the reduction. Codegen for L2_NORM(::Vector{Float64}), counted in-process so the shared machine cannot distort it:

	before	after
<4 x double> / fmul / fadd	20 / 5 / 9	20 / 5 / 9
vfmadd / vaddpd (native)	5 / 4	5 / 4
reduction flags	fast	reassoc, contract
allocations, n = 8…65536	0 B	0 B

Identical vector width, identical flops, identical allocations. The only added work is sqrt's negative-domain check. Interleaved wall clock is +1.0% at n=8, +0.7% at n=65536 — and the norm is O(n) against an O(n²)–O(n³) linear solve per iteration.

Effect

42 (problem × algorithm) cells: successful retcode on a non-finite objective 2 → 0. Both were TrustRegion on least-squares. Separately, the same dead guard was turning Unstable into MaxIters for NewtonRaphson, GaussNewton, PseudoTransient, Broyden, Klement and the default polyalgorithm on least-squares — 7 algorithms × 4 problems now report Unstable instead of burning 100 iterations.

Also fixed: Base.FastMath.abs_fast(::Complex) skips hypot scaling, so L2_NORM(1e200 + 1e200im) returned Inf instead of 1.414e200.

Tests

Testing this is awkward, because the fold depends on the optimizer. A naive @test !isfinite(L2_NORM(x)) passes on unfixed master — the assertion is itself folded. The tests therefore assert observable behaviour (cache.retcode == Unstable, !successful_retcode(sol)), never the predicate.

Negative control, src/common_defaults.jl alone reverted, full GROUP=Core:

NonlinearSolveBase       Non-finite objective protective break | 36 pass 14 fail 50
                         failures are exactly the Base.Fix2(norm,2) rows,
                         Evaluated: ReturnCode.Failure == ReturnCode.Unstable
NonlinearSolveFirstOrder exactly 1 failure in the whole Core suite:
                         TrustRegion: !(successful_retcode(sol))

GROUP=Core passes for both packages on the branch. Runic clean.

Known limitation, stated plainly: the test only bites when the optimizer folds. At -O0, or if LLVM changes, it would pass on unfixed code. It can never false-pass on fixed code, which is the property that matters, but nothing here locks in "no @fastmath in these norms" at source level. A lint or a @code_llvm assertion would be a stronger guarantee if you want one.

Related, not fixed here

DiffEqBase.ODE_DEFAULT_NORM has the same hazard — isnan(ODE_DEFAULT_NORM(u, t)) also compiles to ret i8 0. It appears latent today (the initdt path still reports Unstable by another route), but it is the same trap. Worth its own issue.

🤖 Generated with Claude Code