Dispatch even more to BLAS

[dkarrasch]

This addresses an issue raised in #29634 (comment), by applying the methodology of #33229 to the entire matmul.jl, (edit: and fixes an issue in the syrk_wrapper!). One drawback, just as in #33229, is that also Bool parameters are promoted to the numeric type of the matrices, which is unnecessary for BLAS.gemv! and BLAS.gemm!. We can think about how to work around that, perhaps by working on the types only and introduce some more if branches, but I guess avoiding the generic fallback in some more cases is already a significant improvement.

[dkarrasch]

Earlier commits were failing a specific test:

julia/stdlib/LinearAlgebra/test/matmul.jl

Line 155 in c7e4b99

@test mul!(C, vf, transpose(vf), 2, 3) == 2vf*vf' .+ 3C0

We are passing it on master by pure luck. 2 and 3 are not of the vectors' eltype, and so the multiplication is performed by the generic_matmatmul. Changing them to 2.0 and 3.0, we start to fail the test, because we pass the task to BLAS.syrk!. The reason is that there is no test of symmetry for C when beta is not zero. So we have two options: switch to generic matmul directly when beta is not zero, or test symmetry. I have currently implemented the latter. The symmetry test is most expensive for symmetric matrices, but for those we benefit the most from saving half of the calculations.

[tkf]

Maybe you can use 2 + 0im to avoid dispatching into BLAS? I used this hack here:

julia/stdlib/LinearAlgebra/test/matmul.jl

Lines 584 to 598 in f806717

	@testset "Issue #33214: tiled generic mul!" begin
	n = 100
	A = rand(n, n)
	B = rand(n, n)
	C = zeros(n, n)
	mul!(C, A, B, -1 + 0im, 0)
	D = -A * B
	@test D ≈ C
	# Just in case dispatching on the surface API `mul!` is changed in the future,
	# let's test the function where the tiled multiplication is defined.
	fill!(C, 0)
	LinearAlgebra._generic_matmatmul!(C, 'N', 'N', A, B, LinearAlgebra.MulAddMul(-1, 0))
	@test D ≈ C
	end

[dkarrasch]

The test I mentioned does not dispatch to syrk!, because the matrix C (that we're doubling and adding the symmetric stuff to) is not symmetric. I have added a symmetry check now, and if beta is not zero and C not symmetric, we fall back to gemm! generic_matmatmul!. But maybe we should give the gemm_wrapper! another chance...

[tkf]

Ah, sorry, I misread your earlier comment.

[dkarrasch]

I have extended the coefficient promotion to sym(v/m)! and her(m/v)!. Unfortunately, I did not see how to do something analogous to the gem(v/m)_wrapper! due to the many combinations of arguments.

[dkarrasch]

I think this is ready for review.

[dkarrasch]

BTW, should this be a 1.3 milestone? I guess we wouldn't really want to ship 5-arg mul!, knowing that you can't call it on some Float matrix with alpha and beta equal to some simple Int and end up with BLAS.

[dkarrasch]

Gentle bump.

[tkf]

FWIW, it looks good to me. It's a performance bug fix that does not change the API at all so it would be nice if it can be backported to 1.3 (especially because it's still in RC).

[ViralBShah]

The failing tests appear not to be relevant. Is this otherwise good to merge? Do we need reviews from anyone else?

[dkarrasch]

Unless there are requests for changes, I think this is good to merge.

[tkf]

How about putting the backport label, so that @KristofferC can check if it's OK for backport?

[dkarrasch]

I added the backport label. Let's see...

[ViralBShah]

Should certainly be 1.3.1 or later.

[Jutho]

Thanks @dkarrasch !

[tkf]

This PR probably introduced recent test failures: #32924 (comment)

[dkarrasch]

Wherever this goes (in terms of backport), #33843 must go with it.