# Why is BLAS dot product so much faster than Julia loop?

**URL:** <https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994>\
**Category:** Performance\
**Created:** [August 15, 2020, 2:53pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994 "2020-08-15T14:53:01Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![JeffFessler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jefffessler/32/6650_2.png) [@JeffFessler](https://discourse.julialang.org/u/JeffFessler)\
**Post date:** [August 15, 2020, 2:53pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/1 "2020-08-15T14:53:01Z")

</div>

I wanted to illustrate the beauty of Julia for my class by showing the speed of a simple vector dot product, but to my surprise the Julia loop version was about 5x slower than the built-in `dot()` that (I think) calls BLAS.cdotc\_  
[http://www.netlib.org/lapack/explore-html/dd/db2/cdotc\_8f\_source.html](http://www.netlib.org/lapack/explore-html/dd/db2/cdotc_8f_source.html)  
That FORTAN code uses a loop that looks essentially identical to the Julia code.  
Why is Julia so much slower? Neither `@inbounds` nor `@simd` helped BTW.  
Is Julia calling a fancier version of BLAS with assembly code or such?  
Here is my test code:

```julia
using BenchmarkTools: @btime
using LinearAlgebra: dot
N = 2^16; x = rand(ComplexF32, N); y = rand(ComplexF32, N)
f2(x,y) = dot(y,x)
function f5(x,y) # basic loop
	accum = zero(eltype(x))
	for i in eachindex(x)
		accum += x[i] * conj(y[i])
	end
	return accum
end
@assert f2(x,y) ≈ f5(x,y) # check
# times below are Julia 1.5 on 2017 iMac Pro with 8 threads
@btime f2($x,$y); # dot(y,x) 13.2 us
@btime f5($x,$y); # loop 60.2 us (@inbounds and @simd did not help)

```

This question is kind of the opposite of [How to make Julia slow?](https://discourse.julialang.org/t/how-to-make-julia-slow/28742)

Edit: answered below (BLAS uses cpu-specific assembly code).

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [August 15, 2020, 3:01pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/2 "2020-08-15T15:01:13Z")

</div>

I find a significant improvement with `@simd` + `@inbounds`

```julia
julia> function f5(x,y) # basic loop
               accum = zero(eltype(x))
               for i in eachindex(x)
                       accum += x[i] * conj(y[i])
               end
               return accum
       end
f5 (generic function with 1 method)

julia> function f6(x,y) # basic loop
               accum = zero(eltype(x))
               @simd for i in eachindex(x)
                       @inbounds accum += x[i] * conj(y[i])
               end
               return accum
       end
f6 (generic function with 1 method)

julia> @btime f5($x,$y);
  138.478 μs (0 allocations: 0 bytes)

julia> @btime f6($x,$y);
  56.938 μs (0 allocations: 0 bytes)

julia> @assert f5(x,y) ≈ f6(x,y)

```

---

<div class="post-metadata">

**Author:** ![JeffFessler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jefffessler/32/6650_2.png) [@JeffFessler](https://discourse.julialang.org/u/JeffFessler)\
**Post date:** [August 15, 2020, 3:02pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/3 "2020-08-15T15:02:27Z")

</div>

That is exactly what I had tried.  
Could you please show the time for f2 (dot) for comparison?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 15, 2020, 3:03pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/4 "2020-08-15T15:03:34Z")

</div>

I had a very similar question, with a lot of nice answers

[Simple Mat-Vec multiply (understanding performance, without the bugs)](https://discourse.julialang.org/t/simple-mat-vec-multiply-understanding-performance-without-the-bugs/44762)

my favorite by far was to use @tullio to avoid coding loops at all, just use Einstein tensor notation

```julia
return @tullio x[i]*y[i]

```

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [August 15, 2020, 3:04pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/5 "2020-08-15T15:04:26Z")

</div>

> [@JeffFessler](#):
>
> `@btime f2($x,$y);`

Ah sorry had missed that, here it is

```julia
julia> @btime f2($x,$y);
  44.532 μs (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [August 15, 2020, 3:04pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/6 "2020-08-15T15:04:56Z")

</div>

LoopVectorization.jl also has an example for the scalar product here: [https://github.com/chriselrod/LoopVectorization.jl#dot-product](https://github.com/chriselrod/LoopVectorization.jl#dot-product). Unfortunately, you’ll have to use StructArrays.jl for this to work with complex numbers.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [August 15, 2020, 3:07pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/7 "2020-08-15T15:07:42Z")

</div>

Is Blas running single-threaded?

---

<div class="post-metadata">

**Author:** ![JeffFessler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jefffessler/32/6650_2.png) [@JeffFessler](https://discourse.julialang.org/u/JeffFessler)\
**Post date:** [August 15, 2020, 3:10pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/8 "2020-08-15T15:10:30Z")

</div>

Thanks for the tip. I will look into Tullio.jl, but I am still curious why BLAS is so much faster.

---

<div class="post-metadata">

**Author:** ![JeffFessler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jefffessler/32/6650_2.png) [@JeffFessler](https://discourse.julialang.org/u/JeffFessler)\
**Post date:** [August 15, 2020, 3:11pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/9 "2020-08-15T15:11:14Z")

</div>

I did “top” while running my timing test and it showed 100% cpu, not 400% cpu so I am guessing it is single threaded, but not sure.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [August 15, 2020, 3:15pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/10 "2020-08-15T15:15:22Z")

</div>

You can check with `LinearAlgebra.BLAS.numthreads()` (taking function name from memory.)

---

<div class="post-metadata">

**Author:** ![JeffFessler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jefffessler/32/6650_2.png) [@JeffFessler](https://discourse.julialang.org/u/JeffFessler)\
**Post date:** [August 15, 2020, 3:23pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/11 "2020-08-15T15:23:49Z")

</div>

> [@DNF](#):
>
> LinearAlgebra.BLAS.numthreads()

The closest I found to that was:

```julia
julia-1.5.0κ8> LinearAlgebra.BLAS.vendor()
:openblas64
julia-1.5.0κ8> LinearAlgebra.BLAS.openblas_get_config()
"OpenBLAS 0.3.9 USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=32"

```

I also found `LinearAlgebra.BLAS.set_num_threads(1)` and that led to exactly the same timing results for `dot` and same for setting it to 4 threads, so I think it must be single threaded.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 15, 2020, 3:23pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/12 "2020-08-15T15:23:56Z")

</div>

> [@JeffFessler](#):
>
> curious why BLAS is so much faster.

Because it does in essence all the stuff Tullio is doing under the hood… @avx and soforth.

---

<div class="post-metadata">

**Author:** ![JeffFessler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jefffessler/32/6650_2.png) [@JeffFessler](https://discourse.julialang.org/u/JeffFessler)\
**Post date:** [August 15, 2020, 3:28pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/13 "2020-08-15T15:28:48Z")

</div>

> [@dlakelan](#):
>
> Because it does in essence all the stuff Tullio is doing under the hood… @avx and soforth.

If so, I’d like to see the BLAS source code for that. The code I am seeing doesn’t have any such features:

> <https://github.com/OpenMathLib/OpenBLAS/blob/ce3651516f12079f3ca2418aa85b9ad571c3a391/lapack-netlib/BLAS/SRC/cdotc.f>

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 15, 2020, 3:32pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/14 "2020-08-15T15:32:52Z")

</div>

It’s probably the Fortran compiler that knows people use Fortran to write matrix ops and optimizes the heck out of them

---

<div class="post-metadata">

**Author:** ![onButtonUp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/onbuttonup/32/3098_2.png) [@onButtonUp](https://discourse.julialang.org/u/onButtonUp)\
**Post date:** [August 15, 2020, 3:38pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/15 "2020-08-15T15:38:38Z")

</div>

> [@How to get number of current BLAS threads?](https://discourse.julialang.org/t/how-to-get-number-of-current-blas-threads/32090):
>
> LinearAlgebra.BLAS.set\_num\_threads() sets the number of threads. But, how to get the number of threads (so that I could save it, set another one, then restore the value)? Thanks.

cut and paste worked here . . .

```julia
const get_num_threads = function() # anonymous so it will be serialized when called
    blas = LinearAlgebra.BLAS.vendor()
    # Wrap in a try to catch unsupported blas versions
    try
        if blas == :openblas
            return ccall((:openblas_get_num_threads, Base.libblas_name), Cint, ())
        elseif blas == :openblas64
            return ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
        elseif blas == :mkl
            return ccall((:MKL_Get_Max_Num_Threads, Base.libblas_name), Cint, ())
        end

        # OSX BLAS looks at an environment variable
        if Sys.isapple()
            return tryparse(Cint, get(ENV, "VECLIB_MAXIMUM_THREADS", "1"))
        end
    catch
    end

    return nothing
end

julϊ̇a> using LinearAlgebra

julϊ̇a> get_num_threads()
4

```

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [August 15, 2020, 3:39pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/16 "2020-08-15T15:39:22Z")

</div>

[https://github.com/xianyi/OpenBLAS/blob/d5e6940253b2ee638509de283b8b1d7695fefbbf/kernel/x86\_64/cdot\_microk\_haswell-2.c](https://github.com/xianyi/OpenBLAS/blob/d5e6940253b2ee638509de283b8b1d7695fefbbf/kernel/x86_64/cdot_microk_haswell-2.c)

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [August 15, 2020, 3:41pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/17 "2020-08-15T15:41:27Z")

</div>

Faulty memory there, sorry.

Odd, though, that you can set the number of threads, but not query it. Exactly the opposite of native Julia threads, which you can query but not set. I wonder why.

---

<div class="post-metadata">

**Author:** ![JeffFessler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jefffessler/32/6650_2.png) [@JeffFessler](https://discourse.julialang.org/u/JeffFessler)\
**Post date:** [August 15, 2020, 4:11pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/18 "2020-08-15T16:11:53Z")

</div>

Thank you @yuyichao for the answer: BLAS is using architecture-specific assembly code for the dot product.

I retried my timing results on a linux server, and there I found that `@simd` with `@inbounds` was nearly as fast as as the `dot()` that calls BLAS. Not sure why my iMac “Pro” didn’t benefit from those macros.

Anyway, hopefully someday the Julia/LLVM compiler will do it for us 🙂

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [August 15, 2020, 4:48pm UTC](https://discourse.julialang.org/t/why-is-blas-dot-product-so-much-faster-than-julia-loop/44994/19 "2020-08-15T16:48:19Z")

</div>

I see similar performance to BLAS with `@fastmath` and `@inbounds`.

```julia
N = 2^16
x = rand(ComplexF32, N)
y = rand(ComplexF32, N)
f2(x,y) = dot(y,x)
@fastmath function f5(x,y) # basic loop
    accum = zero(eltype(x))
    @inbounds for i in eachindex(x)
        accum += x[i] * conj(y[i])
    end
    return accum
end
@assert f2(x,y) ≈ f5(x,y) 
@btime f2($x,$y) # dot(y,x) 19.699 us (BLAS)
@btime f5($x,$y) # loop 20.799 us (@fastmath and @inbounds)

```