# Bad performances when using Multithreading and Distributed with heavy LinearAlgebra calculations

**URL:** <https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258>\
**Category:** Performance\
**Tags:** performance, multithreading, distributed\
**Created:** [July 19, 2024, 11:47pm UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258 "2024-07-19T23:47:29Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [July 19, 2024, 11:47pm UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/1 "2024-07-19T23:47:30Z")

</div>

Hello,

I have noticed that Julia multithreading and Distributed scale very badly when calculating relatively heavy computations involving BLAS and LAPACK. As an example, here I show you the calculation of the `expv` function, which calculates `exp(A) * b` (with `A` a square matrix and `b` a vector`) without directly calculating `exp(A)` but rather using the Arnoldi algorithm.

Here I show you the main code of the algorithm, but I don’t think it is necessary to study it in depth. It should be well written (correct me if it is not), involving only in-place matrix-vector multiplications `mul!(y, A, x)`, where `x` and `y` can also be a `view`, and then it involves `exp(B)` where `B` is a 30x30 matrix. I paste it here just for completeness

> **Arnoldi expv**
>
> ```julia
> struct ArnoldiSpace{VT<:AbstractMatrix{<:BlasFloat},HT<:AbstractMatrix{<:BlasFloat},mT<:Integer}
> V::VT
> H::HT
> Hcopy::HT
> m::mT
> end
> 
> function Base.copy(AS::ArnoldiSpace{<:AbstractMatrix{T1},<:AbstractMatrix{T1}}) where {T1<:BlasFloat}
> return ArnoldiSpace(copy(AS.V), copy(AS.H), copy(AS.Hcopy), AS.m)
> end
> 
> function Base.deepcopy(AS::ArnoldiSpace{<:AbstractMatrix{T1},<:AbstractMatrix{T1}}) where {T1<:BlasFloat}
> return ArnoldiSpace(deepcopy(AS.V), deepcopy(AS.H), deepcopy(AS.Hcopy), AS.m)
> end
> 
> function arnoldi_init!(A, b::AbstractVector{T}, V::AbstractMatrix{T}, H::AbstractMatrix{T}) where {T<:BlasFloat}
> v₁ = view(V, :, 1)
> v₂ = view(V, :, 2)
> v₁ .= b
> normalize!(v₁)
> 
> mul!(v₂, A, v₁)
> H[1, 1] = dot(v₁, v₂)
> axpy!(-H[1, 1], v₁, v₂)
> H[2, 1] = norm(v₂)
> return v₂ ./= H[2, 1]
> end
> 
> function arnoldi_step!(A, V::AbstractMatrix{T}, H::AbstractMatrix{T}, i::TI) where {T<:BlasFloat,TI<:Integer}
> vᵢ = view(V, :, i)
> vᵢ₊₁ = view(V, :, i + 1)
> mul!(vᵢ₊₁, A, vᵢ)
> for j in 1:i
> vⱼ = view(V, :, j)
> H[j, i] = dot(vⱼ, vᵢ₊₁)
> axpy!(-H[j, i], vⱼ, vᵢ₊₁)
> end
> β = H[i+1, i] = norm(vᵢ₊₁)
> vᵢ₊₁ ./= H[i+1, i]
> return β
> end
> 
> function arnoldi!(
> AS::ArnoldiSpace{<:AbstractMatrix{T1},<:AbstractMatrix{T1}},
> A,
> b::AbstractVector{T2},
> ) where {T1<:BlasFloat,T2<:BlasFloat}
> n = size(A, 2)
> V = AS.V
> H = AS.H
> m = AS.m
> 
> n == size(V, 1) || throw(DimensionMismatch())
> n == length(b) || throw(DimensionMismatch())
> 
> arnoldi_init!(A, b, V, H)
> for i in 2:m
> arnoldi_step!(A, V, H, i)
> end
> return AS
> end
> 
> function arnoldi(A, b::AbstractVector{T}, m::Integer) where {T<:BlasFloat}
> n = size(A, 2)
> V = similar(b, n, m + 1)
> H = zeros(T, m + 1, m)
> AS = ArnoldiSpace(V, H, copy(H), m)
> return arnoldi!(AS, A, b)
> end
> 
> ### EXPV TOOLS ###
> 
> function expv!(
> x::AbstractVector{T1},
> AS::ArnoldiSpace{<:AbstractMatrix{T1},<:AbstractMatrix{T1}},
> t::T2,
> b::AbstractVector{T1},
> ) where {T1<:BlasFloat,T2<:Union{BlasFloat,BlasInt}}
> H = AS.H
> Hcopy = AS.Hcopy
> V = AS.V
> m = AS.m
> 
> Hm = view(H, 1:m, 1:m)
> Vm = view(V, :, 1:m)
> lmul!(t, Hm)
> 
> Hcopym = view(Hcopy, 1:m, 1:m)
> copyto!(Hcopym, Hm)
> expH = LinearAlgebra.exp!(Hcopym)
> # expH = exp(Hm)
> 
> β = norm(b)
> expHe = expH[:, 1] # the view doesn't work when using copyto! with CUDA
> cache = similar(x, m)
> copyto!(cache, expHe) # just in case expHe is different from cache (e.g., with CUDA)
> mul!(x, Vm, cache)
> lmul!(β, x)
> 
> return x
> end
> 
> function expv!(
> x::AbstractVector{T1},
> A,
> t::T2,
> b::AbstractVector{T1};
> m::Int = min(30, cld(2 * length(b), 3)),
> ) where {T1<:BlasFloat,T2<:Union{BlasFloat,BlasInt}}
> AS = arnoldi(A, b, m)
> return expv!(x, AS, t, b)
> end
> 
> function expv(
> A,
> t::T1,
> b::AbstractVector{T2};
> m::Int = min(30, cld(2 * length(b), 3)),
> ) where {T1<:BlasFloat,T2<:BlasFloat}
> x = similar(b)
> return expv!(x, A, t, b, m = m)
> end
> 
> ```

And then I create a file called `expv.jl` to make some benchmarks

```julia
# expv.jl #

using Distributed
addprocs(Threads.nthreads())

##

@everywhere begin
    using LinearAlgebra
    using LinearAlgebra: BlasFloat, BlasInt
    using SparseArrays
    BLAS.set_num_threads(1)
end

using BenchmarkTools
using JLD2

##

N = parse(Int, ARGS[1])
T = ComplexF64

H = sprand(T, N, N, 0.1/N)
H = H + H' # Hermitianize 

V = sprand(T, N, N, 0.01/N)
V = V + V' # Hermitianize

ψ0 = rand(T, N)
normalize!(ψ0)

t = 0.4

expv(-1im * H, t, ψ0)

function calc_serial(H, V, ψ0, t, g_list)
    res = similar(g_list)
    for i in eachindex(g_list)
        g = g_list[i]
        A = -1im * (H + g * V)
        res[i] = abs2( dot(ψ0, expv(A, t, ψ0)) )
    end
    return res
end

function calc_multithread(H, V, ψ0, t, g_list)
    res = similar(g_list)
    Threads.@threads for i in eachindex(g_list)
        g = g_list[i]
        A = -1im * (H + g * V)
        res[i] = abs2( dot(ψ0, expv(A, t, ψ0)) )
    end
    return res    
end

function calc_distributed(H, V, ψ0, t, g_list)
    res = pmap(g_list) do g
        A = -1im * (H + g * V)
        abs2( dot(ψ0, expv(A, t, ψ0)) )
    end
    return res
end

##

g_list = range(0.0, 0.1, 32)

# calc_serial(H, V, ψ0, t, g_list)

# calc_multithread(H, V, ψ0, t, g_list)

# calc_distributed(H, V, ψ0, t, g_list)

##

if Threads.nthreads() == 1
    bench_serial = @benchmark calc_serial($H, $V, $ψ0, $t, $g_list)
    mean_time_serial = 1e-9 * sum(bench_serial.times) / length(bench_serial.times)
    jldopen("benchmarks.jld2", "a+") do file
        file["mean_time_serial_$(Threads.nthreads())_$N"] = mean_time_serial
    end
end

bench_multithread = @benchmark calc_multithread($H, $V, $ψ0, $t, $g_list)
bench_distributed = @benchmark calc_distributed($H, $V, $ψ0, $t, $g_list)

mean_time_multithread = 1e-9 * sum(bench_multithread.times) / length(bench_multithread.times)
mean_time_distributed = 1e-9 * sum(bench_distributed.times) / length(bench_distributed.times)

jldopen("benchmarks.jld2", "a+") do file
    file["mean_time_multithread_$(Threads.nthreads())_$N"] = mean_time_multithread
    file["mean_time_distributed_$(Threads.nthreads())_$N"] = mean_time_distributed
end

```

And then running the following bash script

```bash
#!/bin/bash

MAX_THREADS=32

for ((threads=1; threads<=MAX_THREADS; threads++)); do
    # for ((N=1000; N<=100000; N*=10)); do
    # julia -t $threads --project expv.jl $N
    # echo "$threads,$N"
    # done
    N=100000
    julia -t $threads --project expv.jl $N
    echo "$threads,$N"
done

```

The benchmarks are very disappointing. Especially for the Distributed case, were I would expect better performances, due to the heavy computations. The Serial case is a clean horizontal line because I did only one benchmark.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/e/9e59745b1657c15d16ce7b934f2ef2e57f1cc9ca.png)

The general system infos are

```julia-repl
julia> versioninfo()
Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × 13th Gen Intel(R) Core(TM) i9-13900KF
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, goldmont)
Threads: 16 default, 0 interactive, 8 GC (on 32 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 16

```

Just to check if multithreading is working properly, I tested a simple matrix-vector multiplication

```julia
function matvec_mul!(y, A, x)
    for i in 1:size(A, 1)
        y[i] = 0
        for j in 1:size(A, 2)
            y[i] += A[i, j] * x[j]
        end
    end
    y
end

function matvec_mul_threaded!(y, A, x)
    Threads.@threads for i in 1:size(A, 1)
        y[i] = 0
        for j in 1:size(A, 2)
            y[i] += A[i, j] * x[j]
        end
    end
    y
end

```

And I get **88ms** (Serial) vs **16ms** (Multithreaded) on a `ComplexF64` 5000x5000 matrix, which seems more reasonable (but still just a factor 5.5) when using 16 threads.

---

<div class="post-metadata">

**Author:** ![dqeeq](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@dqeeq](https://discourse.julialang.org/u/dqeeq)\
**Post date:** [July 20, 2024, 12:00am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/2 "2024-07-20T00:00:38Z")

</div>

Do you think [this](https://discourse.julialang.org/t/most-efficient-implementation-of-covariance-matrix/114797/10) is related?

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [July 20, 2024, 12:06am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/3 "2024-07-20T00:06:24Z")

</div>

Not sure if this applies to here, but BLAS is already [multi-threaded](https://docs.julialang.org/en/v1/manual/performance-tips/#man-multithreading-linear-algebra) which can complicate measurements like this. Probably worth setting the environment variable `OPENBLAS_NUM_THREADS=1` to start with.

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [July 20, 2024, 12:12am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/4 "2024-07-20T00:12:04Z")

</div>

I have just set `BLAS.set_num_threads(1)`, it should be the same, right?

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [July 20, 2024, 3:21am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/5 "2024-07-20T03:21:57Z")

</div>

Yes, I think that is the same.

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [July 20, 2024, 8:14am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/6 "2024-07-20T08:14:42Z")

</div>

What I find surprisingly strange is the Distributed case, which shouldn’t be correlated to the BLAS multithreading.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [July 20, 2024, 9:49am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/7 "2024-07-20T09:49:05Z")

</div>

You could try to use ThreadPinning.jl to see whether this improves things.

> [@albertomercurio](#):
>
> Just to check if multithreading is working properly, I tested a simple matrix-vector multiplication

Note that this naive implementation is memory bound. So I am not surprised that throwing more threads at it only helps so much after a point.

About your `expv` code: It could be that you hit a memory bound as well. In that case improvements are difficult. But I am not an expert on these things.

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [July 20, 2024, 11:30am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/8 "2024-07-20T11:30:02Z")

</div>

I suspect there might be some memory contention here. You have mutiple threads trying to access the same memory all ar once.

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [July 20, 2024, 11:38am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/9 "2024-07-20T11:38:46Z")

</div>

I tried ThreadPinning.jl.

I have just added `pinthreads(:cores)` before calling the multithreaded function, the results are pretty much the same. Here I show the benchmarks with a matrix smaller by a factor 10

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/9/59d2ccd0030df077df2f77ece2bbdea3e7d4556d.png)

Moreover, what do you mean by “hitting a memory bound”?

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [July 20, 2024, 11:42am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/10 "2024-07-20T11:42:24Z")

</div>

> [@mkitti](#):
>
> I suspect there might be some memory contention here. You have mutiple threads trying to access the same memory all ar once.

But this shouldn’t hold for the Distributed case, right? In this case each process should have its own copy of the memory.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 20, 2024, 11:47am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/11 "2024-07-20T11:47:14Z")

</div>

Well, but the bandwidth between the CPU and the memory (RAM) is limited, no matter how many cores you have. Which hardware are you using?

Or are you using Distributed with multiple machines?

- Ryzen 7950X: Max. Memory Bandwidth… 73.4 GB/s
- Apple M3 Max: Memory Bandwidth with 96GB RAM 300.0 GB/s

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [July 20, 2024, 11:49am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/12 "2024-07-20T11:49:07Z")

</div>

> [@albertomercurio](#):
>
> The benchmarks are very disappointing. Especially for the Distributed case, were I would expect better performances, due to the heavy computations

Do you have a gold standard implementation to compare to that you can share?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 20, 2024, 11:50am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/13 "2024-07-20T11:50:45Z")

</div>

Without sharing your `expv` it’s hard to know. There’s a lot of details in how it can and should be done. ExponentialUtilities.jl did a bunch of optimizations here and found that the `balance!` that’s required in the Arnoldi tends to be a limiting factor in some cases, and IIUC that does not multithread well. That’s where GenericSchur.jl and stuff then starts to come in, though that does not fix the multithreading as that library needs an update.

Also, OpenBLAS gets much worse thread performance than MKL, so switching the backends helped there.

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [July 20, 2024, 12:45pm UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/14 "2024-07-20T12:45:39Z")

</div>

> [@albertomercurio](#):
>
> `i9-13900KF`

Your cpu has 8 performance cores, you should perhaps not expect much more speedup with such a highly optimized library as BLAS. Also you have a maximum memory bandwidth of 90 GB/s. With many BLAS operations you might exceed it. You can check such metrics with packages like LIKWID.jl or LinuxPerf.jl

[https://www.intel.com/content/www/us/en/products/sku/230497/intel-core-i913900kf-processor-36m-cache-up-to-5-80-ghz/specifications.html](https://www.intel.com/content/www/us/en/products/sku/230497/intel-core-i913900kf-processor-36m-cache-up-to-5-80-ghz/specifications.html)

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [July 20, 2024, 1:04pm UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/15 "2024-07-20T13:04:29Z")

</div>

> [@ChrisRackauckas](#):
>
> Without sharing your `expv` it’s hard to know

The `expv` code is in the beginning of the thread. Actually, the code I have is much more complex, stil involving Arnoldi but also using LinearSolve.jl and diagonalization. I have chosen the `expv` code just to have a relative clean code to check with a heavy computation.

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [July 20, 2024, 1:08pm UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/16 "2024-07-20T13:08:18Z")

</div>

> [@sgaure](#):
>
> Your cpu has 8 performance cores, you should perhaps not expect much more speedup with such a highly optimized library as BLAS. Also you have a maximum memory bandwidth of 90 GB/s.

Ok, memory bandwidth seems a reasonable bottleneck. Studying the `expv` function, I think that is comes from the `mul!(y, A, x)` functions, since the other “expensive” function should be `exp!`, which is applied to a small 30x30 matrix. Do you think this is correct? Is there a way to avoid such memory efforts?

BTW, 90 GB/s is very large in my opinion. I mean, in the last plot I showed I considered a system of size 10000x10000 of type ComplexF64. This means that a single operation `mul!(y, A, x)` should move at least 10000x16 bytes = 156 KB, which is far below 90 GB/s, even when parallelizing with 16 cores.

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [July 23, 2024, 10:50pm UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/17 "2024-07-23T22:50:30Z")

</div>

I was trying to use LIKWID.jl to see the memory bandwidth usage of certain Julia code. I started using a single core `mul!(y, A, x)` case, where A is a dense matrix and x and y are vectors. If the dimension is `N = 10000`, the benchmarks give

```julia-repl

julia> @benchmark mul!($y, $A, $x)
BenchmarkTools.Trial: 299 samples with 1 evaluation.
 Range (min … max): 16.439 ms … 17.241 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 16.799 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 16.710 ms ± 241.237 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▁▆▆▇█                                                         
  █████▇▅▆▃▂▁▁▁▁▁▁▂▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▂▃▃█▅▄▇▆▃▃▃▅▅█▄▄▄▃▃▄▄█▇▄▄▃▁▃ ▃
  16.4 ms Histogram: frequency by time 17.1 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

While LIKWID gives

```julia-repl
julia> metrics, events = @perfmon "FLOPS_SP" mul!(y, A, x);
ERROR: The selected register PMC0 is in use.
Please run likwid with force option (-f, --force) to overwrite settings

Group: FLOPS_SP
┌───────────────────────────┬───────────┬──────────┬──────────┬──────────┬──────────┬──────────┬──────────┬──────────┬──────────┬────
│ Event │ Thread 1 │ Thread 2 │ Thread 3 │ Thread 4 │ Thread 5 │ Thread 6 │ Thread 7 │ Thread 8 │ Thread 9 │ T ⋯
├───────────────────────────┼───────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼────
│ ACTUAL_CPU_CLOCK │ 1.02241e8 │ 133354.0 │ 323192.0 │ 132258.0 │ 131953.0 │ 133884.0 │ 137479.0 │ 117064.0 │ 128708.0 │ ⋯
│ MAX_CPU_CLOCK │ 8.23769e7 │ 130235.0 │ 313180.0 │ 128100.0 │ 127610.0 │ 130270.0 │ 133035.0 │ 126910.0 │ 129395.0 │ ⋯
│ RETIRED_INSTRUCTIONS │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ ⋯
│ CPU_CLOCKS_UNHALTED │ 7.19432e7 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
│ RETIRED_SSE_AVX_FLOPS_ALL │ 2.25e8 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
│ MERGE │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
└───────────────────────────┴───────────┴──────────┴──────────┴──────────┴──────────┴──────────┴──────────┴──────────┴──────────┴────
                                                                                                                   23 columns omitted
┌──────────────────────┬───────────┬────────────┬────────────┬────────────┬────────────┬────────────┬────────────┬────────────┬──────
│ Metric │ Thread 1 │ Thread 2 │ Thread 3 │ Thread 4 │ Thread 5 │ Thread 6 │ Thread 7 │ Thread 8 │ T ⋯
├──────────────────────┼───────────┼────────────┼────────────┼────────────┼────────────┼────────────┼────────────┼────────────┼──────
│ Runtime (RDTSC) [s] │ 0.0170919 │ 0.0170919 │ 0.0170919 │ 0.0170919 │ 0.0170919 │ 0.0170919 │ 0.0170919 │ 0.0170919 │ 0. ⋯
│ Runtime unhalted [s] │ 0.0292671 │ 3.81736e-5 │ 9.25162e-5 │ 3.78599e-5 │ 3.77725e-5 │ 3.83253e-5 │ 3.93544e-5 │ 3.35105e-5 │ 3.6 ⋯
│ Clock [MHz] │ 4335.71 │ 3577.02 │ 3605.04 │ 3606.75 │ 3612.25 │ 3590.27 │ 3610.05 │ 3222.33 │ ⋯
│ CPI │ Inf │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ ⋯
│ SP [MFLOP/s] │ 13164.1 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
└──────────────────────┴───────────┴────────────┴────────────┴────────────┴────────────┴────────────┴────────────┴────────────┴──────
                                                                                                                   24 columns omitted

```

I don’t understand why I have 2.25 \* 10^8 FLOPS, since the correct one should be 2 \* 10^8. Aside from that, the calculated MFLOPS/s seems reasonable, since I divided the FLOPS by the 16 ms of the benchmark.

Now I do the memory case

```julia-repl
julia> metrics, events = @perfmon "MEM_SP" mul!(y, A, x);
ERROR: The selected register PMC0 is in use.
Please run likwid with force option (-f, --force) to overwrite settings

Group: MEM_SP
┌───────────────────────────┬───────────┬──────────┬──────────┬──────────┬──────────┬──────────┬──────────┬──────────┬──────────┬────
│ Event │ Thread 1 │ Thread 2 │ Thread 3 │ Thread 4 │ Thread 5 │ Thread 6 │ Thread 7 │ Thread 8 │ Thread 9 │ T ⋯
├───────────────────────────┼───────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼──────────┼────
│ ACTUAL_CPU_CLOCK │ 1.0262e8 │ 134993.0 │ 140043.0 │ 136181.0 │ 136959.0 │ 139807.0 │ 135305.0 │ 121054.0 │ 132740.0 │ ⋯
│ MAX_CPU_CLOCK │ 8.27059e7 │ 130760.0 │ 135450.0 │ 131635.0 │ 132475.0 │ 135905.0 │ 132020.0 │ 129045.0 │ 133315.0 │ ⋯
│ RETIRED_INSTRUCTIONS │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ ⋯
│ CPU_CLOCKS_UNHALTED │ 7.19892e7 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
│ RETIRED_SSE_AVX_FLOPS_ALL │ 2.25e8 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
│ MERGE │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
│ DRAM_CHANNEL_0 │ 888725.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
│ DRAM_CHANNEL_1 │ 888544.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
└───────────────────────────┴───────────┴──────────┴──────────┴──────────┴──────────┴──────────┴──────────┴──────────┴──────────┴────
                                                                                                                   23 columns omitted
┌───────────────────────────────────┬───────────┬────────────┬────────────┬────────────┬────────────┬────────────┬────────────┬──────
│ Metric │ Thread 1 │ Thread 2 │ Thread 3 │ Thread 4 │ Thread 5 │ Thread 6 │ Thread 7 │ T ⋯
├───────────────────────────────────┼───────────┼────────────┼────────────┼────────────┼────────────┼────────────┼────────────┼──────
│ Runtime (RDTSC) [s] │ 0.016978 │ 0.016978 │ 0.016978 │ 0.016978 │ 0.016978 │ 0.016978 │ 0.016978 │ 0 ⋯
│ Runtime unhalted [s] │ 0.0293757 │ 3.86428e-5 │ 4.00884e-5 │ 3.89828e-5 │ 3.92056e-5 │ 4.00208e-5 │ 3.87321e-5 │ 3.4 ⋯
│ Clock [MHz] │ 4334.49 │ 3606.45 │ 3611.81 │ 3614.0 │ 3611.6 │ 3593.66 │ 3580.28 │ ⋯
│ CPI │ Inf │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ ⋯
│ SP [MFLOP/s] │ 13252.4 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
│ Memory bandwidth [MBytes/s] │ 26798.2 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
│ Memory data volume [GBytes] │ 0.454981 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ 0.0 │ ⋯
│ Operational intensity [FLOP/Byte] │ 0.494526 │ NaN │ NaN │ NaN │ NaN │ NaN │ NaN │ ⋯
└───────────────────────────────────┴───────────┴────────────┴────────────┴────────────┴────────────┴────────────┴────────────┴──────
                                                                                                                   25 columns omitted

```

I see 26 GB/s of memory bandwidth. So I guess that it is obtained by doing

```julia
sizeof(A) / 1024^3 / 0.0167
22.307127535700086

```

Assuming that every is correct, this means that I could parallelize this operations just 3 times, before reaching the memory bandwidth limits of my hardware (assuming 75 GB/s). Is that right?

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [July 24, 2024, 8:20am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/18 "2024-07-24T08:20:35Z")

</div>

> [@albertomercurio](#):
>
> Assuming that every is correct, this means that I could parallelize this operations just 3 times, before reaching the memory bandwidth limits of my hardware (assuming 75 GB/s). Is that right?

More or less, yes. It’s not necessarily that simple because more threads may cause a different pattern of memory access, but more or less, yes.

---

<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:** [July 24, 2024, 10:07am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/19 "2024-07-24T10:07:05Z")

</div>

Some (hopefully helpful) comments.

> [@albertomercurio](#):
>
> I was trying to use [LIKWID.jl](https://juliahub.com/ui/Packages/General/LIKWID) to see the memory bandwidth usage of certain Julia code.

Great. I’m slightly confused, however. Is your CPU really [this one](https://www.intel.com/content/www/us/en/products/sku/230497/intel-core-i913900kf-processor-36m-cache-up-to-5-80-ghz/specifications.html)? According to your `versioninfo()` output above it is. That would be strange because LIKWID doesn’t support Raptor Lake yet. Did you run the experiment on a different system/CPU?

> [@albertomercurio](#):
>
> I see 26 GB/s of memory bandwidth. So I guess that it is obtained by doing
> 
> ```julia
> sizeof(A) / 1024^3 / 0.0167
> 22.307127535700086
> 
> ```

Your estimation seems reasonable. Every element of `A` must must be read from memory at some point (because `sizeof(Float32) * 10000^2 / 1000^2` = 400 MB, `A` shoudn’t fit into cache). The values of `y` and `x` also need to be read (and written). But those vectors are much smaller and likely reside in cache. (You made a minor mistake though: MB = Mega = 1000^3 not Mebi = 1024^3.)

On a more general note, whether one is bounded by memory access or compute can depend on the specific implementation. For example, a naive implementation of a matrix-matrix product is memory bounded, a highly optimized implementation (say in BLAS) OTOH is bound by compute. Hence, estimating the data transfer isn’t always simple and requires care (and knowledge about the implementation).

> [@albertomercurio](#):
>
> Assuming that every is correct, this means that I could parallelize this operations just 3 times, before reaching the memory bandwidth limits of my hardware (assuming 75 GB/s). Is that right?

Depending on the system you are working on, it can be non-trivial to achieve the maximum memory bandwidth. For example, there might be multiple memory channels that can only be accessed from certain cores. If your 3 parallel computations run on the wrong cores, you are not bounded by the maximal memory bandwidth (of the entire system) but of the same of a single memory channel. In such a case, thread pinning plays a crucial role.

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [July 24, 2024, 10:35am UTC](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258/20 "2024-07-24T10:35:51Z")

</div>

> [@carstenbauer](#):
>
> Is your CPU really [this one](https://www.intel.com/content/www/us/en/products/sku/230497/intel-core-i913900kf-processor-36m-cache-up-to-5-80-ghz/specifications.html)? According to your `versioninfo()` output above it is. That would be strange because LIKWID doesn’t support Raptor Lake yet. Did you run the experiment on a different system/CPU?

Yes, LIKWID.jl doesn’t support my current CPU, so I switched to another machine to do the benchmarks.

> [@carstenbauer](#):
>
> For example, a naive implementation of a matrix-matrix product is memory bounded, a highly optimized implementation (say in BLAS) OTOH is bound by compute

Could you explain better this? What do you mean by “memory bound” and “compute bounds?

[Next page](https://discourse.julialang.org/t/bad-performances-when-using-multithreading-and-distributed-with-heavy-linearalgebra-calculations/117258.md?page=2)