# Are "optimized" single thread algorithms typically different than "optimized" multithread implementations?

**URL:** <https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412>\
**Category:** Performance\
**Tags:** multithreading\
**Created:** [November 17, 2022, 6:09pm UTC](https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412 "2022-11-17T18:09:10Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![JohnS](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johns/32/38331_2.png) [@JohnS](https://discourse.julialang.org/u/JohnS)\
**Post date:** [November 17, 2022, 6:09pm UTC](https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412/1 "2022-11-17T18:09:10Z")

</div>

Thanks for your help everyone!

I’m relatively new to multithreading and I’m probably wrong, but it seems like the fastest function for the desired output using a single thread doesn’t parallelize well. However, a for-loop implementation for the output is faster using multithreading. I’m using a Macbook Pro 16 inch M1 pro.

1. Does this generalize, or is it unique to the functions I’m using?
2. If that’s true, should you build a package optimized for single-thread or multithread performance?

Minimum working example - a Matrix-Scalar multiplication

#Single Thread Optimized Version

```julia
x = collect(1:1:10)
y = rand(length(x),length(x), length(x))

@benchmark dropdims(sum(y .* x', dims=2), dims=2)

BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max): 1.792 μs … 498.375 μs ┊ GC (min … max): 0.00% … 98.50%
 Time (median): 1.879 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 2.200 μs ± 10.861 μs ┊ GC (mean ± σ): 11.95% ± 2.41%
 Memory estimate: 9.12 KiB, allocs estimate: 14.

```

Multi Thread Optimized Version

```julia
x = collect(1:1:10)
y = rand(length(x),length(x), length(x))

function Method2(y::Array{Float64,3}, x::Vector{Int64})
    #Preallocate Output
    output = Matrix{Float64}(undef, length(x), size(y,3))
    
    @fastmath @inbounds @spawn for i=1:size(y,3)
        output[:,i] = @view(y[:,:,i]) * x
    end
    return output
end

@benchmark Method2(y, x)

BenchmarkTools.Trial: 10000 samples with 320 evaluations.
 Range (min … max): 348.697 ns … 41.262 μs ┊ GC (min … max): 0.00% … 95.21%
 Time (median): 626.106 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 800.846 ns ± 2.444 μs ┊ GC (mean ± σ): 18.31% ± 5.89%
  349 ns Histogram: frequency by time 1.05 μs <
Memory estimate: 2.59 KiB, allocs estimate: 15.

```

Thanks for much better-educated perspectives than my own!

---

<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:** [November 17, 2022, 6:30pm UTC](https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412/2 "2022-11-17T18:30:51Z")

</div>

> [@JohnS](#):
>
> it seems like the fastest function for the desired output using a single thread doesn’t parallelize well.

I’m not very experienced with parallelism, but I have normally made the opposite observation: if you want good parallel performance, you first optimize single-threaded performance, and then add parallelism. This is at least what I’ve seen for shared-memory parallelism (multi-threading).

The exception is when the algorithm isn’t inherently easily parallelizable, and you need to change the approach completely to achieve it.

In your examples there are, I think, some confounding factors. Your first example isn’t really optimized, and needlessly allocates an intermediate array. The benchmarking approach is also not optimal. You should wrap that code in a function, like you do with the second.

The other function uses `@spawn`, which I’m unsure whether is correct, I think it should be `@threads`. Currently, I think you are just spawning tasks and not waiting for them to finish.

Also, the second function calls BLAS, which is inherently multi-threaded, unless you disabled it explicitly. Also, by using `*` you are creating unnecessary temporary arrays. Try

```julia
@views mul!(output[:,i], y[:,:,i], x) 

```

instead, to work in-place.

---

<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:** [November 17, 2022, 6:49pm UTC](https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412/3 "2022-11-17T18:49:57Z")

</div>

In a very, very general sense: the runtime of a parallelized program is dominated by the portion of the work that can’t be parallelized. Why is this?

If you have a total runtime T, which can be split up into two parts, the parallelizable part T\_{parallel} (abbreviated T\_p) and the non-parallelizable (or serial) part T\_{serial} (abbreviated T\_s) such that T = T\_s + T\_p, then as you increase the number of threads/processes to work on your program, T\_p goes to 0 (as that time is divided between all processes/threads), leaving T = T\_s. This is because we can use p = T\_p / T, i.e. the proportion of the total runtime that we can parallelize, to write the above like:

T = T\_s + T\_p = (1-p)T + pT

If we now throw n parallel processors at the problem, we get

(1-\frac{p}{n})T + \frac{p}{n}T

We can see that as n \to \infty, the first half tends to T, while the second half tends to 0, resulting in the fact that the non-parallelizable portion dominates the runtime (in the limit).

This is called [Amdahl’s Law](https://en.wikipedia.org/wiki/Amdahl's_law).

The trouble is of course, that p depends heavily on the algorithm you want to calculate, some of which are easily parallelizeable (i.e. have a very large p close to 1), while others are very serial in nature (i.e. have a small p close to 0). Regardless of which it is though, it’s pretty much always beneficial to start optimizing the serial portion of the parallelizable work itself, since that’s the part that will benefit from parallelization.

---

<div class="post-metadata">

**Author:** ![jar1](https://avatars.discourse-cdn.com/v4/letter/j/c0e974/32.png) [@jar1](https://discourse.julialang.org/u/jar1)\
**Post date:** [November 17, 2022, 7:41pm UTC](https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412/4 "2022-11-17T19:41:37Z")

</div>

Parallel and sequential performance do pull in different directions.

- In sequential code, delay is strongly correlated with the number of operations, so that’s what you minimize. In parallel code, we perform redundant operations to reduce communication/synchronization.
- In sequential code, minimize space usage. In parallel code, use extra space to reduce synchronization. For example, fast sequential code reuses buffers; parallel code doesn’t.
- In sequential code, it’s common to have a single accumulator modified from a single thread. In parallel code, there is no single accumulator - values are accumulated in trees.

For more info [Organizing Functional Code for Parallel Execution; or, foldl and foldr Considered Slightly Harmful - YouTube](https://www.youtube.com/watch?v=JbvnhuiHYxA) and [JuliaCon 2016 | A functional approach to High Performance Computing | Erik Schnetter - YouTube](https://www.youtube.com/watch?v=h4VpT0YcnH8).

---

<div class="post-metadata">

**Author:** ![JohnS](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johns/32/38331_2.png) [@JohnS](https://discourse.julialang.org/u/JohnS)\
**Post date:** [November 17, 2022, 10:10pm UTC](https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412/5 "2022-11-17T22:10:26Z")

</div>

Thanks for your feedback!  
Function wrapping was negligible, so I didn’t bother, but will in the future!  
I’d read that re @spawn, but it returns the correct answer consistently…  
Love your specific solution for the mwe regardless!

---

<div class="post-metadata">

**Author:** ![jmair](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmair/32/35117_2.png) [@jmair](https://discourse.julialang.org/u/jmair)\
**Post date:** [November 18, 2022, 9:01am UTC](https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412/6 "2022-11-18T09:01:43Z")

</div>

I would prefer `Threads.@threads` for a `for` loop, it will always be correct unless you have race conditions and will likely be cheaper than using `@spawn` in most cases.

---

<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:** [November 18, 2022, 9:14am UTC](https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412/7 "2022-11-18T09:14:46Z")

</div>

Since `@threads` is using the same functionality as `@spawn` internally, I’m not sure I follow how `@threads` can be faster. Both are using the julia native task machinery.

Also, correctness depends on the code being wrapped - if you share data between iterations, `@threads` will not save you from race conditions. Both `@threads` and `@spawn` have the same core limitations & difficulties, one is just more convenient for throwing on a loop than manually dividing the work.

---

<div class="post-metadata">

**Author:** ![jmair](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmair/32/35117_2.png) [@jmair](https://discourse.julialang.org/u/jmair)\
**Post date:** [November 18, 2022, 9:21am UTC](https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412/8 "2022-11-18T09:21:56Z")

</div>

I may be wrong but `@threads` should only spawn (at most) as many threads as you have cores, and reuse them for later iterations? Is this right?

EDIT:

I’m beginning to think the example was wrong anyway:

> [@JohnS](#):
>
> ```julia
> @fastmath @inbounds @spawn for i=1:size(y,3)
> output[:,i] = @view(y[:,:,i]) * x
> end
> 
> ```

Will just spawn a thread to do the for loop sequentially, right? And worse, it won’t wait for the result so benchmarking will be quick. I don’t `@threads` would be faster, but it would be correct. I was under the impression that there would be a thread spawned for each part of the loop, which I thought would be slower for a large number of iterations than `@threads`, but with `@spawn` we also need to manually wait for each thread as well.

---

<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:** [November 18, 2022, 9:33am UTC](https://discourse.julialang.org/t/are-optimized-single-thread-algorithms-typically-different-than-optimized-multithread-implementations/90412/9 "2022-11-18T09:33:41Z")

</div>

> [@jmair](#):
>
> I may be wrong but `@threads` should only spawn (at most) as many threads as you have cores

Minor nitpick: neither `@spawn` nor `@threads` is spawning new threads (in the OS sense of the word). Both are creating julia native `Task` objects that are scheduled to run on OS threads by the julia runtime. It has been said that calling it `@threads` instead of `@tasks` has been a mistake, but there’s nothing we can do about it now since it’s supported API and changing that would be breaking.

> [@jmair](#):
>
> and reuse them for later iterations? Is this right?

Generally yes, but that’s something you can do with `@spawn` as well if you want, e.g. by querying `Threads.nthreads()` as the upper limit for the loop. Like I said, both use the same mechanism under the hood. One is just more manual, but may be more beneficial for nontrivial parallelization tasks that don’t neatly fit into a loop.

In the case of OP, it looks like they just spawned a single task in which the `for` loop runs, without actually waiting on the result, which is why the function “just” takes nanoseconds to run in the first place. It should have been written like

```julia
function Method2(y::Array{Float64,3}, x::Vector{Int64})
    #Preallocate Output
    output = Matrix{Float64}(undef, length(x), size(y,3))
    
    task = @spawn @fastmath @inbounds for i=1:size(y,3)
        output[:,i] = @view(y[:,:,i]) * x
    end
    wait(task)
    return output
end

```

which should have about the same runtime as the first version due to not being parallelized at all (the `@fastmath` probably doesn’t do a whole lot either - its usefulness is greatly overstated). You can of course wrap the `for i=1:size(y,3)` in a task spawning loop for dividing the work up and parallelize that way, but that’s more or less exactly what `@threads` is doing internally. One is not going to be faster than the other and both are equal in terms of possible race conditions.