Why doesn’t Julia fma automatically?

[giordano]

@ParadaCarleton

I think your post is the first one to bring up an actual case of these kinds of bugs causing problems in real life.

There were lots of concrete bogus examples in the other thread, I think this was one of the nicest ones:

X * y + z does not automatically use FMA instruction

Using FMA automatically can cause surprises. Consider the following: function f(a,b,c) @assert a*b ≥ c return sqrt(a*b-c) end a = 1.0 + 0.5^27 b = 1.0 - 0.5^27 c = 1.0 f(a,b,c) It works as expected. Now suppose the compiler automatically uses FMA inside f, like this: function f(a,b,c) @assert a*b ≥ c return sqrt(fma(a,b,-c)) end Now f(a,b,c) throws and exception, because fma(a,b,-c) is negative, despite the assertion a*b ≥ c.

And this was your “basic/case easy win a*x + b”

[vchuravy]

@ParadaCarleton

Were the bugs all caused by people writing things like

function kernel(rho, T)
  P = rho[1] * T[1]

  if (abs(P - P) > 1e-16)
    @cuprintf("diff = %.16e\n", P - P)
  end
  nothing
end

rho = CuArray([1e-1])
T = CuArray([300.0])
@launch CUDA() kernel(rho, T, threads=1, blocks=1)

with the output

diff = 1.6653345369377348e-15

So abs(P-P) being different from P-P

	mul.f64 	%fd4, %fd1, %fd3;
	neg.f64 	%fd5, %fd4;
	fma.rn.f64 	%fd2, %fd1, %fd3, %fd5;
	abs.f64 	%fd6, %fd2;

[stevengj]

@Oscar_Smith

the problem is that compilers don’t work one line at a time. You basically only have 2 options: allow re-association and contraction or prevent them.

clang allows finer-grained control than this — -ffp-contract=on will only fuse within a single statement in the source code, as allowed by the C99 standard, whereas -ffp-contract=fast will fuse across statements.

[mikmoore]

Plugging this tangentially related issue for selectable @fastmath flags for those wishing to exercise restraint in their fastmath. We already have @simd with less @fastmath slated for v1.10 so that it isn’t such a hazard.