Test failure in `LinearAlgebra/addmul`

[giordano]

On 4633607ce9

The global RNG seed was 0xadd355ff29e71246580cec9b4759733c.

Error in testset LinearAlgebra/addmul:
Test Failed at ~/repo/julia/usr/share/julia/stdlib/v1.12/LinearAlgebra/test/addmul.jl:175
  Expression: ≈(collect(returned_vec), α * Ac * xc + β * yc, rtol = rtol)
   Evaluated: [0.00016286048013114396, 0.0] ≈ [0.0001627827388785974, 0.0] (rtol=0.0003452669770922512)

Offending test is https://github.com/JuliaLang/julia/blob/4633607ce9b9f077f32f89f09a136e04389bbac2/stdlib/LinearAlgebra/test/addmul.jl#L169-L175 I could reliably reproduce this with Base.runtests("LinearAlgebra/addmul"; seed=0xadd355ff29e71246580cec9b4759733c) on both x86_64-linux and aarch64-linux. Note: this is not a test reintroduced by #55775

[jishnub]

Here is a simpler and a more concrete reproducer extracted from the testsets:

julia> using LinearAlgebra

julia> A = Diagonal(Float32[1.2805837, -0.83571815]);

julia> xc = [-10, 0];

julia> yc = [4.479540055608876, 0.0];

julia> α = -0.32607044710547134; β = -0.9321140763694151;

julia> y1 = α * A * xc + β * yc
2-element Vector{Float64}:
 0.0001627827388785974
 0.0

julia> y2 = mul!(copy(yc), A, xc, α, β)
2-element Vector{Float64}:
 0.00016286048013114396
 0.0

julia> norm(y1 - y2)/norm(y1)
0.000477576757106528

[jishnub]

This seems to be an issue with associativity. The operation A * xc * α is computed as

mat_vec_scalar(A, x, γ) = A * (x * γ)

whereas the in-place multiplication operates elementwise, and performs the first multiplication before the second one.

julia> (A * xc * α)[1]
4.175605124232544

julia> A * xc * α == A * (xc * α)
true

julia> ((A * xc) * α)[1]
4.175605201973797

julia> ((A * xc) * α)[1] == A[1,1] * xc[1] * α
true

This is exacerbated by the fact that A * xc * α is very close in magnitude and opposite in sign to yc * β, so the rounding error has a large impact on the rtol to which the results match.

julia> (A * xc * α)[1], (yc * β)[1]
(4.175605124232544, -4.175442341493666)

julia> y1 = (A * xc * α)[1] + (yc * β)[1]
0.0001627827388785974

julia> y2 = ((A * xc) * α)[1] + (yc * β)[1]
0.00016286048013114396

julia> norm(y1 -y2)/norm(y1)
0.000477576757106528

Since we are comparing small numbers, perhaps we need an atol here.