Overzealous `fmul contract` marking on `*` outside of `muladd`?

[mbauman]

I think Julia is being overzealous in marking the first fmul as contact-able here:

julia> smat = (0.3759580090201759539070280879968777298927307128906250000000000000, -0.2925452398379911889136906211206223815679550170898437500000000000, -0.8792456187560546698733787707169540226459503173828125000000000000, -0.6857980393775208183271274720027577131986618041992187500000000000, 0.5502662724003037908460100879892706871032714843750000000000000000, -0.4763277008999405315314845665852772071957588195800781250000000000, 0.6231666106584448083793859041179530322551727294921875000000000000, 0.7820641355456769971965513832401484251022338867187500000000000000, 0.0062500602924510650915124188031768426299095153808593750000000000);

julia> svec = (0.5513258229231927654012679340667091310024261474609375000000000000, -0.8283034402416151742443162220297381281852722167968750000000000000, -0.0997659654489912867125767093057220336049795150756835937500000000)
(0.5513258229231928, -0.8283034402416152, -0.09976596544899129)

julia> f3(a, b) = muladd(a[9], b[3], muladd(a[6], b[2], a[3] * b[1]))
f3 (generic function with 1 method)

julia> @code_llvm f3(smat, svec)

; Function Signature: f3(NTuple{9, Float64}, Tuple{Float64, Float64, Float64})
;  @ REPL[4]:1 within `f3`
define double @julia_f3_1940(ptr nocapture noundef nonnull readonly align 8 dereferenceable(72) %"a::Tuple", ptr nocapture noundef nonnull readonly align 8 dereferenceable(24) %"b::Tuple") #0 {
top:
; ┌ @ tuple.jl:31 within `getindex`
   %"a::Tuple[9]_ptr" = getelementptr inbounds i8, ptr %"a::Tuple", i64 64
   %"b::Tuple[3]_ptr" = getelementptr inbounds i8, ptr %"b::Tuple", i64 16
   %"a::Tuple[6]_ptr" = getelementptr inbounds i8, ptr %"a::Tuple", i64 40
   %"b::Tuple[2]_ptr" = getelementptr inbounds i8, ptr %"b::Tuple", i64 8
   %"a::Tuple[3]_ptr" = getelementptr inbounds i8, ptr %"a::Tuple", i64 16
; └
; ┌ @ float.jl:493 within `*`
   %"a::Tuple[3]_ptr.unbox" = load double, ptr %"a::Tuple[3]_ptr", align 8
   %"b::Tuple.unbox" = load double, ptr %"b::Tuple", align 8
   %0 = fmul contract double %"a::Tuple[3]_ptr.unbox", %"b::Tuple.unbox"
; └
; ┌ @ float.jl:496 within `muladd`
   %"a::Tuple[6]_ptr.unbox" = load double, ptr %"a::Tuple[6]_ptr", align 8
   %"b::Tuple[2]_ptr.unbox" = load double, ptr %"b::Tuple[2]_ptr", align 8
   %1 = fmul contract double %"a::Tuple[6]_ptr.unbox", %"b::Tuple[2]_ptr.unbox"
   %2 = fadd contract double %0, %1
   %"a::Tuple[9]_ptr.unbox" = load double, ptr %"a::Tuple[9]_ptr", align 8
   %"b::Tuple[3]_ptr.unbox" = load double, ptr %"b::Tuple[3]_ptr", align 8
   %3 = fmul contract double %"a::Tuple[9]_ptr.unbox", %"b::Tuple[3]_ptr.unbox"
   %4 = fadd contract double %2, %3
   ret double %4
; └
}

Why is that first %0 = fmul contract? That's the * I wrote — and it is unexpected that passing that result to muladd would induce a contract without @fastmath.

We noticed this because having that * marked contract means that LLVM rearranges this operation on my ARM64 CPU. Specifically, f3 gives different answers in a function than the same code does at top-level.

julia> f3(smat, svec)
-0.0908304842736846

julia> (a, b) = (smat, svec);

julia> muladd(a[9], b[3], muladd(a[6], b[2], a[3] * b[1]))
-0.09083048427368462

julia> # it's actually doing this:
       muladd(a[9], b[3], muladd(a[3], b[1], a[6] * b[2]))
-0.0908304842736846

I see this behavior both on v1.10.5 and today's nightly (1.12.0-DEV.1209 @ 55c40ce).

julia> versioninfo()
Julia Version 1.12.0-DEV.1209
Commit 55c40ce52eb (2024-09-16 17:25 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M1 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 6 virtual cores)

[mbauman]

It feels even spookier if I explore a few more function definitions:

julia> f3(smat, svec)
-0.0908304842736846

julia> f3a(a, b) = (x = a[3] * b[1]; muladd(a[9], b[3], muladd(a[6], b[2], x)))
f3a (generic function with 4 methods)

julia> f3a(smat, svec)
-0.0908304842736846

julia> f3b(a, b) = (x = a[3] * b[1]; string(x); muladd(a[9], b[3], muladd(a[6], b[2], x)))
f3b (generic function with 2 methods)

julia> f3b(smat, svec)
-0.09083048427368462

[mbauman]

Julia v1.9 did not mark this first fmul as contract:

▶ julia +1.9 -q

julia> smat = (0.3759580090201759539070280879968777298927307128906250000000000000, -0.2925452398379911889136906211206223815679550170898437500000000000, -0.8792456187560546698733787707169540226459503173828125000000000000, -0.6857980393775208183271274720027577131986618041992187500000000000, 0.5502662724003037908460100879892706871032714843750000000000000000, -0.4763277008999405315314845665852772071957588195800781250000000000, 0.6231666106584448083793859041179530322551727294921875000000000000, 0.7820641355456769971965513832401484251022338867187500000000000000, 0.0062500602924510650915124188031768426299095153808593750000000000);

julia> svec = (0.5513258229231927654012679340667091310024261474609375000000000000, -0.8283034402416151742443162220297381281852722167968750000000000000, -0.0997659654489912867125767093057220336049795150756835937500000000)
(0.5513258229231928, -0.8283034402416152, -0.09976596544899129)

julia> f3(a, b) = muladd(a[9], b[3], muladd(a[6], b[2], a[3] * b[1]))
f3 (generic function with 1 method)

julia> @code_llvm f3(smat, svec)

;  @ REPL[3]:1 within `f3`
define double @julia_f3_108([9 x double]* nocapture noundef nonnull readonly align 8 dereferenceable(72) %0, [3 x double]* nocapture noundef nonnull readonly align 8 dereferenceable(24) %1) #0 {
top:
; ┌ @ tuple.jl:29 within `getindex`
   %2 = getelementptr inbounds [9 x double], [9 x double]* %0, i64 0, i64 8
   %3 = getelementptr inbounds [3 x double], [3 x double]* %1, i64 0, i64 2
   %4 = getelementptr inbounds [9 x double], [9 x double]* %0, i64 0, i64 5
   %5 = getelementptr inbounds [3 x double], [3 x double]* %1, i64 0, i64 1
   %6 = getelementptr inbounds [9 x double], [9 x double]* %0, i64 0, i64 2
   %7 = getelementptr inbounds [3 x double], [3 x double]* %1, i64 0, i64 0
; └
; ┌ @ float.jl:410 within `*`
   %8 = load double, double* %6, align 8
   %9 = load double, double* %7, align 8
   %10 = fmul double %8, %9
; └
; ┌ @ float.jl:413 within `muladd`
   %11 = load double, double* %4, align 8
   %12 = load double, double* %5, align 8
   %13 = fmul contract double %11, %12
   %14 = fadd contract double %10, %13
   %15 = load double, double* %2, align 8
   %16 = load double, double* %3, align 8
   %17 = fmul contract double %15, %16
   %18 = fadd contract double %14, %17
; └
  ret double %18
}

[giordano]

I wonder if this is related to #49405 just because it touches the llvm-muladd file and first appeared in 1.10, but the point of that PR was to reduce automatic unsafe operations, so I'd hope it's not that.

[gbaraldi]

It's technically that PR but I'm not super sure what that pass is supposed to do. Basically it looks for floating point operations that are only used in other floating point operations with contract and adds contract to them. So it's intentionally doing this right or wrong. The change is that it used to look for fast and now looks for contract which is what muladd does.
Specifically this

/**
 * Combine
 * ```
 * %v0 = fmul ... %a, %b
 * %v = fadd contract ... %v0, %c
 * ```
 * to
 * %v0 = fmul contract ... %a, %b
 * %v = fadd contract ... %v0, %c
 * when `%v0` has no other use
 */

[mbauman]

Yeah, I wasn't initially sure this was wrong (hence the ? in the title), but I'm slowly becoming convinced. This seems not great. What happens if we completely disable that pass? Looks like it dates way back to #22262 and is explicitly tested — which is surely how/why @vchuravy updated that pass in his change there.

But the tests themselves are quite contrived LLVM, not higher-level Julia tests. Perhaps at one point the muladd_fast intrinsic only marked the add as fast? It looks like it now marks both directly (the boolean args to math_builder are fast and contract flags, respectively).

julia/src/intrinsics.cpp

Lines 1507 to 1513 in 02549d5

	case muladd_float: {
	// LLVM 5.0 can create FMA in the backend for contractible fmul and fadd
	// Emitting fmul and fadd here since they are easier for other LLVM passes to
	// optimize.
	auto mathb = math_builder(ctx, false, true);
	return mathb().CreateFAdd(mathb().CreateFMul(x, y), z);
	}

[mbauman]

Aha, I see, it's the @simd loop that requires it. Without that pass, Julia does this:

julia> function f(a, b)
           s = zero(eltype(a))
           @inbounds @simd for i in 1:length(a)
               s += a[i] * b[i]
           end
           return s
       end
f (generic function with 1 method)

julia> @code_llvm f(Float64[], Float64[])
; Function Signature: f(Array{Float64, 1}, Array{Float64, 1})
# ...
vector.body:                                      ; preds = %vector.body, %vector.ph
# ...
; │┌ @ float.jl:493 within `*`
    %10 = fmul <2 x double> %wide.load, %wide.load34
    %11 = fmul <2 x double> %wide.load31, %wide.load35
    %12 = fmul <2 x double> %wide.load32, %wide.load36
    %13 = fmul <2 x double> %wide.load33, %wide.load37
; │└
; │┌ @ float.jl:491 within `+`
    %14 = fadd reassoc contract <2 x double> %vec.phi, %10
    %15 = fadd reassoc contract <2 x double> %vec.phi28, %11
    %16 = fadd reassoc contract <2 x double> %vec.phi29, %12
    %17 = fadd reassoc contract <2 x double> %vec.phi30, %13
; │└

So that's what this is doing when it's working as expected. I see.

[giordano]

Yeah, that was the observation in #49405 (comment) right before I modified the muladd pass.

[mbauman]

Hah, sorry Mosè! Should've read further than the first post in that PR. I'm slowly catching up here. :)

[giordano]

As I understand this, in simd context the user explicitly allows for non-standard behaviour (realloc/contract after the PR, full fastmath before it) so automatic muladd fusion should be acceptable also without the muladd function. The problem is that this is happening also outside of simd context. Probably the issue is that we don't define very well this context?

[mbauman]

Yeah, I think that's just it — I think I do want/expect fmas inside @simd. But I don't think this is working as I'd expect, even in the @simd context. There are just so many non-localities going on here, even within @simd. For example, here's a fun trio of functions:

julia> function f1(a, b)
           s1 = zero(eltype(a))
           s2 = zero(eltype(a))
           @inbounds @simd for i in 1:length(a)
               x = a[i] * b[i]
               s1 += x
               s2 += 2.0
           end
           return s1, s2
       end
f1 (generic function with 1 method)

julia> function f2(a, b)
           s1 = zero(eltype(a))
           s2 = zero(eltype(a))
           @inbounds @simd for i in 1:length(a)
               x = a[i] * b[i]
               s1 += x
               s2 += 2.0x
           end
           return s1, s2
       end
f2 (generic function with 1 method)

julia> (x, y) = ([0.7932039553981265, 0.6781862407025683, 0.5496472683398627, 0.6088071625708916], [0.6045129003064795, 0.1115602608320766, 0.493350882104838, 0.5285717302662082]);

julia> f1(x,y)[1]
1.1481278775183397

julia> f2(x,y)[1]
1.1481278775183394

julia> f1(x,y)[1] == f2(x,y)[1]
false

julia> test_fs(x, y) = f1(x,y)[1] == f2(x,y)[1]
test_fs (generic function with 1 method)

julia> test_fs(x, y)
true

It's perhaps worth noting that I intentionally contrived f1 and f2 and was still caught off guard by the infinite while f1(x,y)[1] == f2(x,y)[1] loop to find the problematic x and y.

[mikmoore]

That contract is being added outside of muladd is an issue, of course, and a correctness bug at that.

Specifically, f3 gives different answers in a function than the same code does at top-level.

It feels even spookier if I explore a few more function definitions:

I know you know this, but it seems that a reminder might be in order. The compiler can choose not to act on contract flags, so you still cannot expect reproducible results while using muladd (you must instead explicitly use fma or unfused * +). Even with this issue fixed, you should still expect inconsistent results from these examples.

[vchuravy]

 // Emitting fmul and fadd here since they are easier for other LLVM passes to 

This statement may no longer be true and maybe worthwile to check...

I am unsure why the @simd pass doesn't touch the fmul and instead relies on the MulAdd pass