# \`mul\` dispatch to BLAS incomplete?

**URL:** <https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831>\
**Category:** Performance\
**Tags:** linearalgebra, numerics\
**Created:** [November 20, 2021, 3:22pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831 "2021-11-20T15:22:45Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 20, 2021, 3:22pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/1 "2021-11-20T15:22:45Z")

</div>

When translating the benchmark from [Outperformed by Matlab](https://discourse.julialang.org/t/outperformed-by-matlab/71603)

```julia
using LinearAlgebra
using MKL
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)
        BLAS.axpy!(x[n],temp,Cx)
    end
    return Cx
end

function kr!(C,A,B)
    @views for n = 1:size(A,1)
        _kron!(C[:,n],A[n,:],B[n,:])
    end
    return C
end

function kr(A,B)
    C = Matrix{Float64}(undef,size(A,2)*size(B,2),size(A,1))
    kr!(C,A,B)
    return C
end

function normal2!(C,Cx,x,A,B,temp,n1,n2)
    @views kr!(temp,A[n1:n2,:],B[n1:n2,:])
    BLAS.gemm!('N','T',1.,temp,temp,1.,C)
    @views BLAS.gemv!('N',1.,temp,x[n1:n2],1.,Cx)
end
function normal2!(C,Cx,x,A,B,blocksize)
    n1 = 1
    if blocksize < size(A,1)
        temp = Matrix{Float64}(undef,size(A,2)*size(B,2),blocksize)
        @views while n1 <= size(A,1)-blocksize
            n2 = n1+blocksize-1
            normal2!(C,Cx,x,A,B,temp,n1,n2)
            n1 += blocksize
        end
    end
    n2 = size(A,1)
    temp = Matrix{Float64}(undef,size(A,2)*size(B,2),n2-n1+1)
    normal2!(C,Cx,x,A,B,temp,n1,n2)
    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,100)
@assert C ≈ CNormal
@assert Cx ≈ CxNormal

@btime normal1!($C,$Cx,$x,$A,$B)

# for blocksize in [10, 100, 1000, 10000, 100000]
for blocksize in [10, 100, 1000]
    @btime normal2!($C,$Cx,$x,$A,$B,$blocksize)
end

```

to use `mul` instead

```julia
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,:])
        # mul!(C,temp,temp',1.,1.)
        BLAS.ger!(1.,temp,temp,C)
        mul!(Cx,temp,x[n],1.,1.)
    end
    return Cx
end

function kr!(C,A,B)
    @views for n = 1:size(A,1)
        _kron!(C[:,n],A[n,:],B[n,:])
    end
    return C
end

function kr(A,B)
    C = Matrix{Float64}(undef,size(A,2)*size(B,2),size(A,1))
    kr!(C,A,B)
    return C
end

function normal2!(C,Cx,x,A,B,temp,n1,n2)
    @views kr!(temp,A[n1:n2,:],B[n1:n2,:])
    mul!(C,temp,temp',1.,1.)
    @views mul!(Cx,temp,x[n1:n2],1.,1.)
end
function normal2!(C,Cx,x,A,B,blocksize)
    n1 = 1
    if blocksize < size(A,1)
        temp = Matrix{Float64}(undef,size(A,2)*size(B,2),blocksize)
        @views while n1 <= size(A,1)-blocksize
            n2 = n1+blocksize-1
            normal2!(C,Cx,x,A,B,temp,n1,n2)
            n1 += blocksize
        end
    end
    n2 = size(A,1)
    temp = Matrix{Float64}(undef,size(A,2)*size(B,2),n2-n1+1)
    normal2!(C,Cx,x,A,B,temp,n1,n2)
    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,100)
@assert C ≈ CNormal
@assert Cx ≈ CxNormal

@btime normal1!($C,$Cx,$x,$A,$B)

# for blocksize in [10, 100, 1000, 10000, 100000]
for blocksize in [10, 100, 1000]
    @btime normal2!($C,$Cx,$x,$A,$B,$blocksize)
end

```

I noticed the following oddities.

1. `mul!(C,temp,temp',1.,1.)` seems to compute the same as `BLAS.ger!(1.,temp,temp,C)` but is _slooow_.
2. `mul!(Cx,temp,x[n],1.,1.)` seems to work (which I like very much), although this is not documented?
3. Runtime performance using BLAS

```julia
  371.686 ms (1 allocation: 6.38 KiB)
  100.720 ms (4 allocations: 125.16 KiB)
  62.885 ms (4 allocations: 1.22 MiB)
  69.160 ms (4 allocations: 12.21 MiB)

```

looks significantly different to the Julia `mul` version

```julia
  336.500 ms (1 allocation: 6.38 KiB)
  922.396 ms (4 allocations: 125.16 KiB)
  128.485 ms (4 allocations: 1.22 MiB)
  56.164 ms (4 allocations: 12.21 MiB)

```

which really surprises me. This is on Julia 1.6.4, 1.7.0-rc3 behaves similar.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 20, 2021, 4:02pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/2 "2021-11-20T16:02:43Z")

</div>

Regarding 3:

Not sure, but could it be [https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/matmul.jl#L351](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/matmul.jl#L351), which redirects to `syrk` instead of `gemm`? There is even [https://github.com/JuliaLang/julia/issues/659](https://github.com/JuliaLang/julia/issues/659). Interesting, seems numerically sensible.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 20, 2021, 4:11pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/3 "2021-11-20T16:11:21Z")

</div>

Regarding 1:

It seems to me that support for `ger` and `syr` is missing in [https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/matmul.jl](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/matmul.jl).

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [November 20, 2021, 10:10pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/4 "2021-11-20T22:10:24Z")

</div>

This would be a really good issue or PR.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 20, 2021, 10:35pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/5 "2021-11-20T22:35:25Z")

</div>

Thanks for your feedback, filed an [issue](https://github.com/JuliaLang/julia/issues/43175).

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 21, 2021, 4:45pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/6 "2021-11-21T16:45:47Z")

</div>

Is there some documentation which mappings of high level Julia functions to BLAS are intended to work?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [November 21, 2021, 4:50pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/7 "2021-11-21T16:50:35Z")

</div>

I don’t think so. Generally we just promise that whatever we decide to do will be reasonably accurate and reasonably fast.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 21, 2021, 4:56pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/8 "2021-11-21T16:56:42Z")

</div>

Which it sometimes isn’t, so that we have [discussions](https://discourse.julialang.org/t/outperformed-by-matlab/71603). I’m at a loss, what to do.

Edit: it seems for high performance applications it then would be appropriate to depend on the low level BLAS wrappers.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 21, 2021, 5:30pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/9 "2021-11-21T17:30:09Z")

</div>

Which brings problems of it’s own, because

```julia
using BenchmarkTools, Random, LinearAlgebra
elt = Float64
n = 32
Random.seed!(1234)
# A = rand(elt, n, n)
A = rand(elt, n, n)
b = rand(elt, n)
println(BLAS.get_config())
t1 = BLAS.syr!('U', elt(1), b, copy(A))
t2 = BLAS.syr!('U', elt(1), b, Symmetric(copy(A)))

```

fails with

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

```

Shouldn’t at least something like this work then?

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [November 21, 2021, 5:57pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/10 "2021-11-21T17:57:18Z")

</div>

BLAS can’t handle Julia types like `Symmetric` and you’re calling rather directly into BLAS. So I wouldn’t expect this to work, no. But what do I know 🙂

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 21, 2021, 6:02pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/11 "2021-11-21T18:02:00Z")

</div>

I’m asking the question because `BLAS` [`syr`](http://www.netlib.org/lapack/explore-html/d6/d30/group __single__ blas__level2_ga7b8a99048765ed2bf7c1e770bff0b622.html) specifically expects a symmetric matrix memory layout AFAIU, and I’d have expected `Symmetric` to exactly provide this?

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 21, 2021, 7:19pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/12 "2021-11-21T19:19:07Z")

</div>

OK, so what about this then?

```julia
using BenchmarkTools, Random, LinearAlgebra

elt = Float64
n = 32
Random.seed!(1234)
# A = rand(elt, n, n)
A = rand(elt, n, n)
b = Symmetric(rand(elt, n, n))
println(BLAS.get_config())
t1 = BLAS.syrk!('U', 'N', elt(1), b, elt(1), copy(A))
t2 = BLAS.syrk!('U', 'N', elt(1), b, elt(1), Symmetric(copy(A)))

```

failing again with

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

```

And you really expect us to answer MATLAB user questions regarding performance?

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 21, 2021, 10:36pm UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/13 "2021-11-21T22:36:32Z")

</div>

Filed another issue [here](https://github.com/JuliaLang/julia/issues/43182).

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [November 22, 2021, 12:49am UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/14 "2021-11-22T00:49:18Z")

</div>

`Symmetric` is only a wrapper. If you want to use BLAS.syrk! on your matrix `A`, you should use `A.data`.  
Example :

```julia
elt = Float64
n = 32
Random.seed!(1234)
A = rand(elt, n, n)
B = zeros(elt, n, n)
A = triu(A) # A = tril(A)
S = Symmetric(A, :U) # S = Symmetric(A, :L)
 
println(BLAS.get_config())
t1 = BLAS.syrk!(S.uplo, 'N', elt(1), S.data, elt(1), B)

```

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 22, 2021, 1:07am UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/15 "2021-11-22T01:07:25Z")

</div>

Yes, this seems to work:

```julia
using BenchmarkTools, Random, LinearAlgebra
elt = Float64
n = 32
Random.seed!(1234)
# A = rand(elt, n, n)
A = rand(elt, n, n)
b = rand(elt, n)
t1 = BLAS.syr!('U', elt(1), b, copy(A))
println(BLAS.get_config())
t2 = BLAS.syr!('U', elt(1), b, Symmetric(copy(A)).data) 

```

Edit: after @amontoison `s correction this

```julia
using BenchmarkTools, Random, LinearAlgebra

elt = Float64
n = 32
Random.seed!(1234)
# A = rand(elt, n, n)
A = rand(elt, n, n)
b = rand(elt, n, n)
println(BLAS.get_config())
t1 = BLAS.syrk!('U', 'N', elt(1), b, elt(1), copy(A))
t2 = BLAS.syrk!('U', 'N', elt(1), b, elt(1), Symmetric(copy(A)).data)

```

works too.

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [November 22, 2021, 1:13am UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/16 "2021-11-22T01:13:18Z")

</div>

`b` is a Symmetric matrix, you forgot `.data`.

```julia
elt = Float64
n = 32
Random.seed!(1234)
# A = rand(elt, n, n)
A = rand(elt, n, n)
b = Symmetric(rand(elt, n, n))
println(BLAS.get_config())
t1 = BLAS.syrk!('U', 'N', elt(1), b.data, elt(1), copy(A))
t2 = BLAS.syrk!('U', 'N', elt(1), b.data, elt(1), Symmetric(copy(A)).data)

```

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 22, 2021, 1:14am UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/17 "2021-11-22T01:14:40Z")

</div>

Edit: Nonsense on my side…

Hm. But the [spec](http://www.netlib.org/lapack/explore-html/db/dc9/group __single__ blas__level3_gae953a93420ca237670f5c67bbde9d9ff.html) allows for general matrices there?

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [November 22, 2021, 1:29am UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/18 "2021-11-22T01:29:48Z")

</div>

BLAS routines require general arrays (Vector{T} or Matrix{T}).  
I don’t know why you use a symmetric matrix for `b`, I just fixed the error of your previous message 🙂

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 22, 2021, 1:34am UTC](https://discourse.julialang.org/t/mul-dispatch-to-blas-incomplete/71831/19 "2021-11-22T01:34:06Z")

</div>

My mistake, excellent: this did the job, you have earned another solution point!