Inaccurate Matrix Multiplication

fma is a red herring here. Any way of computing a matmul will have cases with arbitrarily high forward error. On average, fma lowers forward error, but not always.
For example, in this case

julia> prevfloat(exp2(100))*nextfloat(exp2(100))-exp2(100)^2
0.0
julia> fma(prevfloat(exp2(100)),nextfloat(exp2(100)),-exp2(100)^2)
1.7840596158824495e44
julia> big(prevfloat(exp2(100)))*big(nextfloat(exp2(100)))-big(exp2(100))^2
1.7840596158824494551820448904901809527586816e+44

The fma based answer gives a correctly rounded answer, while the version without fma is completely wrong.

3 Likes