Fix type instability in generic matmul

[simonster]

For C::AbstractVecOrMat{R}, A::AbstractVecOrMat{T}, B::AbstractVecOrMat{S}, we were sometimes initializing the intermediate sum as zero(R), but it wouldn't stay that way if R was narrower than T and S.

This promotes the intermediate sums (s and Ctmp in the code) to the promotion of R and zero(T)*zero(S) + zero(T)*zero(S). Thus, for A_mul_B!(C::Matrix{Float32}, A::Matrix{Float64}, B::Matrix{Float64}, the intermediate sums are Float64. The other possibility would be to enforce that the intermediate sums are the same type as R (i.e. Float32 in this example) but it seemed more conservative to potentially use higher precision than lower precision.

Before:

julia> C = Array(Float32, 100, 100);
       A = rand(Float64, 100, 100);
       B = rand(Float64, 100, 100);

julia> @time A_mul_B!(C, A, B);
  0.031763 seconds (2.11 M allocations: 32.154 MB, 9.28% gc time)

After:

julia> @time A_mul_B!(C, A, B);
  0.000870 seconds (10 allocations: 448 bytes)

[andreasnoack]

LGTM

[simonster]

The RootInt type in one of the tests doesn't have zero defined, so that test is failing. In principle I guess we should be setting up the zero values for the blocked algorithm the same way as we currently do for the naive algorithm.

[simonster]

If A or B is empty, maybe we could just skip the loop entirely and fill!(C, zero(R))?

[andreasnoack]

Yes. I think this check can then be moved outside the transpose branches which will make the code a bit shorter and simpler.

[andreasnoack]

In principle I guess we should be setting up the zero values for the blocked algorithm the same way as we currently do for the naive algorithm.

I.e. check for non-emptiness and then use the elements of the arrays instead of the array element types?

[tkelman]

can this be tested against via @inferred or looking at the @code_warntype output to prevent it from regressing? cc @timholy on the zero issue

[timholy]

On the zero issue, there's no mathematical reason to prevent defining zero(RootInt). OTOH maybe it's a useful test case for such situations.

[simonster]

Updated to:

• Bail out early if A or B is empty and just fill C with zeros, instead of putting a bunch of branches everywhere.
• Use zero based on elements themselves rather than types for the blocked case. (I just use the first elements in both arrays, since if they're isbits the types will be the same. If you tried hard enough, I think you could construct a type where this is not mathematically correct because zero does different things for different elements, but then the current code wouldn't work either.)
• Add a test that allocation doesn't scale with input size for A_mul_B!(::Matrix{Float32}, ::Matrix{Float64}, ::Matrix{Float64}). It's possibly more fragile than is ideal, but it was the easiest approach I could think of that would catch a regression here.

[tkelman]

this fails when inlining is off https://build.julialang.org/builders/coverage_ubuntu14.04-x64/builds/404/steps/Run%20non-inlined%20tests/logs/stdio

[tkelman]

Bump! This is still blocking coverage. https://build.julialang.org/builders/coverage_ubuntu14.04-x64/builds/461/steps/Run%20non-inlined%20tests/logs/stdio

[simonster]

Sorry, I totally missed this the first time around. Filed #17327 to revert the test. Will look into adding a benchmark instead tomorrow.