`mul` dispatch to BLAS incomplete?

[goerch]

Based on the discussion here (https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831) I'd like to note the following observations:

• mul!(C,temp,temp',1.,1.) with temp being a vector seems to compute the same as BLAS.ger!(1.,temp,temp,C) or BLAS.syr!(1.,temp,C) but does not dispatch to them AFAIU.
• mul!(Cx,temp,x,1.,1.) with temp being a vector and x being a scalar seems to work (which I like very much), although this is not documented?
• Runtime behaviour for mul!(C,temp,temp',1.,1.) could be surprising for users because of dispatch to BLAS.syrk! AFAIU, so should this be documented?

[N5N3]

My local bench with MKL shows that ger! is faster than gemm! only in limited circumstances. For small inputs (<64*64), ger! is generally slower. While for large inputs(>1024*1024), their time costs are almost the same.
Benchmark code:

using BenchmarkTools, Random
elt = Float64
for n in (32, 64, 128, 256, 512, 1024, 2048)
    Random.seed!(1234)
    A = rand(elt, n, n)
    b = rand(elt, n)
    c = rand(elt, n)
    t1 = @benchmark BLAS.ger!(elt(1), $b, $c, a) setup = (a = copy($A))
    t2 = @benchmark BLAS.gemm!('N', 'C', true, $b, $c, true, a) setup = (a = copy($A))
    println(elt, '|', n, ':', ratio(median(t1), median(t2)))
end

Result:

ComplexF32|32:TrialRatio(154.18%)
ComplexF32|64:TrialRatio(121.27%)
ComplexF32|128:TrialRatio(108.24%)
ComplexF32|256:TrialRatio(94.21%)
ComplexF32|512:TrialRatio(245.15%)
ComplexF32|1024:TrialRatio(100.99%)
ComplexF32|2048:TrialRatio(99.74%)

ComplexF64|32:TrialRatio(125.49%)
ComplexF64|64:TrialRatio(97.30%)
ComplexF64|128:TrialRatio(95.75%)
ComplexF64|256:TrialRatio(85.43%)
ComplexF64|512:TrialRatio(100.64%)
ComplexF64|1024:TrialRatio(101.00%)
ComplexF64|2048:TrialRatio(101.02%)

Float32|32:TrialRatio(174.11%)
Float32|64:TrialRatio(163.56%)
Float32|128:TrialRatio(128.68%)
Float32|256:TrialRatio(110.58%)
Float32|512:TrialRatio(96.79%)
Float32|1024:TrialRatio(100.00%)
Float32|2048:TrialRatio(100.38%)

Float64|32:TrialRatio(156.26%)
Float64|64:TrialRatio(125.62%)
Float64|128:TrialRatio(106.88%)
Float64|256:TrialRatio(100.89%)
Float64|512:TrialRatio(100.37%)
Float64|1024:TrialRatio(100.20%)
Float64|2048:TrialRatio(100.18%)

[goerch]

@N5N3: your benchmark indeed doesn't reproduce the extent of the problem. Therefore a somewhat reduced test case from OP:

using LinearAlgebra
using BenchmarkTools

function _kron!(C::T1,A::T2,B::T3) where {T1 <: AbstractVector,T2 <: AbstractVector,T3 <: AbstractVector}
    @views for i = 1:size(A,1)
        mul!(C[(i-1)*size(B,1)+1:i*size(B,1)],A[i],B)
    end
    return C
end

function _kron(A::T1,B::T2) where {T1 <: AbstractVector,T2 <: AbstractVector}
    C = zeros(size(A,1) * size(B,1))
    _kron!(C,A,B)
    return C
end

function normal1!(C,Cx,x,A,B)
    temp = Vector{Float64}(undef,size(A,2)*size(B,2))
    @views for n = 1:size(A,1)
        _kron!(temp,A[n,:],B[n,:])
        BLAS.ger!(1.,temp,temp,C)
        mul!(Cx,temp,x[n],1.,1.)
    end
    return Cx
end

function normal2!(C,Cx,x,A,B)
    temp = Vector{Float64}(undef,size(A,2)*size(B,2))
    @views for n = 1:size(A,1)
        _kron!(temp,A[n,:],B[n,:])
        mul!(C,temp,temp',1.,1.)
        mul!(Cx,temp,x[n],1.,1.)
    end
    return Cx
end

N = 10000
R = 20
M = 40

A = rand(N,M)
B = rand(N,R)
x = rand(N,)

C = zeros(M*R,M*R)
Cx = zeros(M*R,)

normal1!(C,Cx,x,A,B)
CNormal = copy(C)
CxNormal = copy(Cx)

C = zeros(M*R,M*R)
Cx = zeros(M*R,)
normal2!(C,Cx,x,A,B)
@assert C ≈ CNormal
@assert Cx ≈ CxNormal

println(BLAS.get_config())
@btime normal1!(copy($C),copy($Cx),$x,$A,$B)
@btime normal2!(copy($C),copy($Cx),$x,$A,$B)

with results

LBTConfig([ILP64] libopenblas64_.dll)
  774.735 ms (4 allocations: 4.90 MiB)
  10.895 s (4 allocations: 4.90 MiB)

and

LBTConfig([ILP64] mkl_rt.1.dll)
  323.367 ms (4 allocations: 4.90 MiB)
  8.694 s (4 allocations: 4.90 MiB)

[N5N3]

Yes, this example makes much more sense.
If temp<:AbstractVector, ishermitian and copytri! in syrk_wrapper! seems too heavy and we'd better call LinearAlgebra.gemm_wrapper directly in this case, i.e. add a size check in mul!:

@inline function mul!(C::StridedMatrix{T}, A::StridedVecOrMat{T}, tB::Transpose{<:Any,<:StridedVecOrMat{T}},
    alpha::Number, beta::Number) where {T<:BlasFloat}
    B = tB.parent
    if A === B && size(A, 2) >= size(A, 1) ÷ 5 # need to be tuned
        return syrk_wrapper!(C, 'N', A, MulAddMul(alpha, beta))
    else
        return gemm_wrapper!(C, 'N', 'T', A, B, MulAddMul(alpha, beta))
    end
end

after the above Revise, your Benchmark example gives

LBTConfig([ILP64] mkl_rt.1.dll)
  354.466 ms (4 allocations: 4.90 MiB)
  312.852 ms (4 allocations: 4.90 MiB)

on my PC.

Of cource making mul! call ger! is another solution, but the above change is simpler and covers more situations I think.

Edit: Also notice that both BLAS.ger! and gemm_wrapper! is less accurate than syrk_wrapper!, as syrk_wrapper! ensures the output's symmetry, while the other 2 not. On the other hand, syr! has no performance difference with syrk! as the bottlenecks are ishermitian and copytri!.

[goerch]

Thank you very much @N5N3! Just a minor note: when I try to modify your benchmark like so

using BenchmarkTools, Random, LinearAlgebra
println(BLAS.get_config())
elt = Float64
for n in (32, 64, 128, 256, 512, 1024, 2048)
    Random.seed!(1234)
    A = rand(elt, n, n)
    b = rand(elt, n)
    t1 = @benchmark BLAS.ger!(elt(1), $b, $b, a) setup = (a = copy($A))
    # t2 = @benchmark BLAS.gemm!('N', 'C', true, $b, $b, true, a) setup = (a = copy($A))
    t2 = @benchmark mul!(a, $b, $b', elt(1), elt(1)) setup = (a = copy($A))
    println(elt, '|', n, ':', ratio(median(t1), median(t2)))
end
using MKL
println(BLAS.get_config())
elt = Float64
for n in (32, 64, 128, 256, 512, 1024, 2048)
    Random.seed!(1234)
    A = rand(elt, n, n)
    b = rand(elt, n)
    t1 = @benchmark BLAS.ger!(elt(1), $b, $b, a) setup = (a = copy($A))
    # t2 = @benchmark BLAS.gemm!('N', 'C', true, $b, $b, true, a) setup = (a = copy($A))
    t2 = @benchmark mul!(a, $b, $b', elt(1), elt(1)) setup = (a = copy($A))
    println(elt, '|', n, ':', ratio(median(t1), median(t2)))
end

I see

LBTConfig([ILP64] libopenblas64_.dll)
Float64|32:TrialRatio(62.27%)
Float64|64:TrialRatio(66.79%)
Float64|128:TrialRatio(1019.51%)
Float64|256:TrialRatio(125.37%)
Float64|512:TrialRatio(94.07%)
Float64|1024:TrialRatio(78.71%)
Float64|2048:TrialRatio(77.39%)
LBTConfig([ILP64] mkl_rt.1.dll)
Float64|32:TrialRatio(65.50%)
Float64|64:TrialRatio(77.40%)
Float64|128:TrialRatio(84.44%)
Float64|256:TrialRatio(83.91%)
Float64|512:TrialRatio(98.47%)
Float64|1024:TrialRatio(98.57%)
Float64|2048:TrialRatio(99.22%)

and in this form

using BenchmarkTools, Random, LinearAlgebra
println(BLAS.get_config())
elt = Float64
for n in (32, 64, 128, 256, 512, 1024, 2048)
    Random.seed!(1234)
    A = rand(elt, n, n)
    b = rand(elt, n)
    # t1 = @benchmark BLAS.ger!(elt(1), $b, $b, a) setup = (a = copy($A))
    t1 = @benchmark BLAS.syr!('U', elt(1), $b, a) setup = (a = copy($A))
    # t2 = @benchmark BLAS.gemm!('N', 'C', true, $b, $b, true, a) setup = (a = copy($A))
    t2 = @benchmark mul!(a, $b, $b', elt(1), elt(1)) setup = (a = copy($A))
    println(elt, '|', n, ':', ratio(median(t1), median(t2)))
end
using MKL
println(BLAS.get_config())
elt = Float64
for n in (32, 64, 128, 256, 512, 1024, 2048)
    Random.seed!(1234)
    A = rand(elt, n, n)
    b = rand(elt, n)
    # t1 = @benchmark BLAS.ger!(elt(1), $b, $b, a) setup = (a = copy($A))
    t1 = @benchmark BLAS.syr!('U', elt(1), $b, a) setup = (a = copy($A))
    # t2 = @benchmark BLAS.gemm!('N', 'C', true, $b, $b, true, a) setup = (a = copy($A))
    t2 = @benchmark mul!(a, $b, $b', elt(1), elt(1)) setup = (a = copy($A))
    println(elt, '|', n, ':', ratio(median(t1), median(t2)))
end

I see

BTConfig([ILP64] libopenblas64_.dll)
Float64|32:TrialRatio(1858.98%)
Float64|64:TrialRatio(1595.74%)
Float64|128:TrialRatio(780.61%)
Float64|256:TrialRatio(120.92%)
Float64|512:TrialRatio(84.84%)
Float64|1024:TrialRatio(48.80%)
Float64|2048:TrialRatio(40.66%)
LBTConfig([ILP64] mkl_rt.1.dll)
Float64|32:TrialRatio(49.25%)
Float64|64:TrialRatio(50.02%)
Float64|128:TrialRatio(65.62%)
Float64|256:TrialRatio(133.55%)
Float64|512:TrialRatio(68.63%)
Float64|1024:TrialRatio(67.26%)
Float64|2048:TrialRatio(51.94%)

which look slightly off for OpenBLAS. Besides performance what do you think about the numerical semantics of using symm vs gemm in this case?

Edit: just noticed:

Edit: Also notice that both BLAS.ger! and gemm_wrapper! is less accurate than syrk_wrapper!, as syrk_wrapper! ensures the output's symmetry, while the other 2 not. On the other hand, syr! has no performance difference with syrk! as the bottlenecks are ishermitian and copytri!.

We agree. That is why I don't like the proposed change

@inline function mul!(C::StridedMatrix{T}, A::StridedVecOrMat{T}, tB::Transpose{<:Any,<:StridedVecOrMat{T}},
    alpha::Number, beta::Number) where {T<:BlasFloat}
    B = tB.parent
    if A === B && size(A, 2) >= size(A, 1) ÷ 5 # need to be tuned
        return syrk_wrapper!(C, 'N', A, MulAddMul(alpha, beta))
    else
        return gemm_wrapper!(C, 'N', 'T', A, B, MulAddMul(alpha, beta))
    end
end

too much...

[N5N3]

As for your modifition, since A = rand(elt, n, n) should not be symetric, so mul!(a, $b, $b', elt(1), elt(1)) always fallbacks to gemm_wrapper!. The extra cost should come from issymetric. (I can't verify, as adding @profview before @benchmark changes the benchmark result mysteriously)
As for symmetric accuracy, under my own use cases, I think we'd better do such optimization to C <: Symmetric{T, <:StridedMatrix{T}} rather than C <: StridedMatrix. But there might be another trade-off.

Edit: I think replacing size(A, 2) >= size(A, 1) ÷ 5 with size(A, 2) >= 2 should be as accurate as the current status, although I didn't do much test.

[goerch]

Thanks for your feedback, @N5N3. Should something like

    t2 = @benchmark mul!(a, $b, $b', elt(1), elt(1)) setup = (a = Symmetric(copy($A)))

work then? I get

LBTConfig([ILP64] libopenblas64_.dll)
ERROR: LoadError: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix

and

LBTConfig([ILP64] mkl_rt.1.dll)
ERROR: LoadError: MethodError: no method matching strides(::Symmetric{Float64, Matrix{Float64}}

with that change.

Edit:
Just back from studying matmul tests. This is the only reference to Symmetric I could find.

Edit:
Nevermind, just found the symmetric tests. And these are the only mul! tests for Symmetric I could find.

[N5N3]

Yes, they haven’t been implemented yet. The reason might be it’s limited usage. And the worry about new ambiguity problem. One might argue that user could use C .+= b .* b’ instead (although this have some performance issue currently)

[goerch]

One might argue that user could use C .+= b .* b’ instead (although this have some performance issue currently)

This is the main reason why I'm here to bother you;)

[N5N3]

You could refer to JuliaLang/julia#43153 for more information. (Unfortunately, the easiest fix is not safe for general usage.)

[goerch]

Going back to general discussion at https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831 for now to get more opinions...

[goerch]

The discussion on discourse showed a sensible way out.