15 ms·

Composable multi-threaded parallelism in Julia

289 points by KenoFischer 7y ago

sjkelly 7y ago
Live stream here: https://www.youtube.com/watch?v=HfiRnfKxI64 https://www.youtube.com/watch?v=HfiRnfKxI64
eigenspace 7y ago
This is huge! I was testing out the master branch a few days ago and the parallelism improvements were amazing.
- pkofod 7y ago
  I think it's great that they went ahead and released a binary for people to try it out without having to compile from source.
  - staticfloat 7y ago
    Thanks! We do our best to make sure that people who want to get ahold of the latest Julia binaries are always able to. Our downloads page [0] always offers the latest stable builds, whereas the nightlies page [1] allows you to download binaries off of the tip of the `master` branch, which usually lag behind the github repository by only a couple of hours. [0] https://julialang.org/downloads/ https://julialang.org/downloads/ [1] https://julialang.org/downloads/nightlies.html https://julialang.org/downloads/nightlies.html
mindB 7y ago
It's going to be really neat to watch what the Julia package ecosystem comes up with now that they have access to truly composable parallelism at the language level. Up until now, most parallelism in the language had to live at the application level rather than being in library code since parallelism either relied on spawning new processes or the Threads.@threads macro which couldn't be nested.
StefanKarpinski 7y ago
It's hard to overstate what a big deal this is (and how much effort went into it). The only mainstream languages with this kind of composable parallelism are Go and Erlang, both of which target writing servers, rather than high-performance computing, and which therefore optimize for I/O throughput rather than CPU throughput. The only systems which have supported compute-oriented composable parallelism like this are Cilk and TBB, which are both really excellent and successful—and provided major inspiration for this work. However, they both seem to have fallen short of widespread adoption partly due to being non-free (until recently) and also, I suspect, because they are non-standard C/C++ language/compiler extensions, which seem to never catch on that broadly. This is a brave new world for high-performance technical computing, not just for Julia.
- eigenspace 7y ago
  As a physicist, it's been really cool seeing the awareness of Julia in the community gradually increase. Large swathes of the physics community have been migrating their code to Python over the past few years but it really seems now that Julia is getting so many amazing features and packages that it can't be ignored.
  - benjaminva 7y ago
    Same here, I worked in a computational biophysics department with lots of python/bash/R and I was the only one who wrote lots of high-performance code in Julia. People were curious about the language but it was still very much unknown. Hope we will see a broader adoption of Julia in the future - it's just that much better for the stuff we do on a daily basis.
- Arnavion 7y ago
  >The only mainstream languages with this kind of composable parallelism are Go and Erlang, both of which target writing servers, rather than high-performance computing, and which therefore optimize for I/O throughput rather than CPU throughput. Rust, while perhaps not as mainstream and not necessarily targeting scientific computing users, has data-bound parallelization like this via the rayon library.
3JPLW 7y ago
I'm super excited about this! For me the most exciting part is these four sentences: > The model is portable and free from low-level details. You don’t need to explicitly start and stop threads, and you don’t even need to know how many processors or threads there are (though you can find out if you want). > The model is nestable and composable: you can start parallel tasks that call library functions that themselves start parallel tasks, and everything works. Your CPUs will not be over-subscribed with threads. It means that you can put threaded parallelism inside libraries without worrying. It means that Julia itself will start making its builtins threaded. It means that Julia itself is much more threadsafe by default. It means that packages can implement threaded algorithms — and you can use multiple threaded packages together at the same time without worrying. And it means that you can put threaded parallelism in your application code and all these things will work together beautifully.
- thekhatribharat 7y ago
  Well, only if your code and your libraries do not share state.
- jstrong 7y ago
  still have to be careful about data races and thread synchronization, though: julia> xs = [1] 1-element Array{Int64,1}: 1 julia> for j in 1:1024 d = @spawn xs[1] += 1 end julia> xs 1-element Array{Int64,1}: 686
  - deleted 7y ago
    [deleted]
polskibus 7y ago
Is this something similar to PLINQ in C# ?
- ska 7y ago
  IIUC it is more general, and operates at a lower level.
- mr_crankypants 7y ago
  It seems like it's more like C#'s PLINQ and TPL merged into a single abstraction.
- manigandham 7y ago
  It's basically the Task-based async model.
- davidacoder 7y ago
  You could use this to build something like PLINQ, but this is at a lower level in the stack than PLINQ.
cutler 7y ago
Why isn't some big company putting its weight behind Julia considering what it has to offer in the field of data science? I just don't get why it's still playing second fiddle to Python. Tech is full of inertia, paradoxically.
- 6thaccount2 7y ago
  I think some finance groups are probably putting in some donations if I had to guess.
  - cutler 7y ago
    That's not the same as Google's support for Python, Go and Kotlin or Microsoft's support for Typescript and maybe soon Rust. I would have thought the case for supporting Julia was very strong considering Python's performance and concurrency problems.
    - 6thaccount2 7y ago
      Python's performance is less of an issue when you use libraries like Numpy that are basically running optimized C code. Depending on what you're doing, Python is usually great at performance where the time it saves the user writing code is more valuable than the performance itself (a lot of exploratory scientific research).
      - DNF2 7y ago
        Numpy is great if you can vectorize your code in a straightforward manner. But when you start doing anything non-standard it becomes difficult. I can tell you that programmer productivity can plummet drastically when you need to get everything into vectorized form. Also, custom classes and types do not vectorize well at all.
        6thaccount2 7y ago
        Agreed.
argd678 7y ago
Does anyone have any recommendations for using Julia on real world projects at this point? Is it mature enough to use and debug on a day to day basis for example? Mainly I’ve been using R for context.
- npr11 7y ago
  I work at company using Julia in production. What sort of recommendations are you looking for? :)
  - afiori 7y ago
    Not the parent, but I am interested in a good list of libraries that might be non obvious to find.
    - ddragon 7y ago
      You can use the search tools to find libraries for the categories you're interested at: https://juliaobserver.com/packages https://juliaobserver.com/packages https://pkg.julialang.org/docs/ https://pkg.julialang.org/docs/ As a personal answer, I really like unitful for keeping track of physical units and libraries like MLStyle.jl or Match.jl which add stuff like pattern matching to the language (though they are more general helper libraries compared to stuff like Turing.jl and DifferentialEquations.jl): https://painterqubits.github.io/Unitful.jl/stable/highlights/ https://painterqubits.github.io/Unitful.jl/stable/highlights... https://github.com/thautwarm/MLStyle.jl https://github.com/thautwarm/MLStyle.jl
  - argd678 7y ago
    Mainly, am I going to spend hours fussing with the basics of the language. For example in R, it’s too slow, so I had to write Rcpp, only to discover there’s no real way to debug it from RStudio, down into the rabbit hole, lldb, Rinside, everything that can crash does, terse errors messages with out line numbers, opaque SEXP types XCode can’t decode etc. So basically are the normal everyday use cases fraught with quirks you need to know more about the internals of Julia or can you just be a user? If that makes sense.
karmakaze 7y ago
Not only is this huge for Julia, but they've just thrown down the gauntlet. The status quo has been upset. I expect Julia to start eating everyone's lunch starting with Python. Every language can use good concurrency & parallelism support and this is the biggest news for all dynamic languages.
cshenton 7y ago
Julia is already insanely productive for CS and statistics research and this will only make it more so. Super excited to see how this will be used, especially in PPLs, as algorithms like PMCMC benefit hugely from this kind of cooperative, interruptible concurrency.
joehl 7y ago
Maybe I misunderstand, but: In R allocating and copying huge return vectors for small changes on input vectors ruins performance. My understanding of Julia was that one advantage over R is the possiblity to mutate input vectors in place. Referring to the psort! example for Julia's multithreading I wonder why "We simply need to allocate one initially-empty array per thread" which to me looks like allocating huge vectors for small changes down the call recursion (and wastes memory by the number of threads). If Julia can mutate psort!s first argument global vector 'v' from multiple threads, why can't Julia also mutate a single allocated 'temp' argument global vector? Bonus: then psort! could just merge from one argument to the other and save copying back on each recursion level (along the lines of Segewick's nocopy mergesort)?
- joehl 7y ago
  OK, I found out myself: Julia can mutate a single temp buffer across all threads, and implementing Sedgewicks nocopy sort (with ping-pong merge) is faster, particularly when passing down preallocated buffer memory to the single threaded sorting Routine once N <= BIG_THRESHOLD. The following are measurements run with JULIA_NUM_THREADS=4 on a Intel(R) Core(TM) i5-4300M CPU @ 2.60GHz, 2594 MHz, 2 physical cores, 4 logical cores, replication code follows. Note that not only runtime of nocopy mergesort ist better, it also causes less garbage collection time. Conclusion: it would make sense for Julia's sorting interface to allow passing in preallocated buffer memory and to use nocopy mergesort behind the curtain. Of course true parallel sorting also requires a parallel merge. Best julia> # 4 threads, mergesort, dynamic allocation julia> @time for i=1:100 c = copy(a); psort!(c); end; 166.270792 seconds (756.58 k allocations: 89.470 GiB, 15.34% gc time) julia> # 4 threads, mergesort, preallocation per thread, re-allocation in single threaded subroutine julia> @time for i=1:100 c = copy(a); tsort!(c); end; 143.973485 seconds (658.00 k allocations: 36.972 GiB, 5.89% gc time) julia> # 4 threads, nocopy mergesort, one preallocation, re-allocation in single threaded subroutine julia> @time for i=1:100 c = copy(a); ksort!(c); end; 129.715467 seconds (654.51 k allocations: 37.312 GiB, 6.43% gc time) julia> # 4 threads, nocopy mergesort, one preallocation reused in single threaded subroutine julia> @time for i=1:100 c = copy(a); jsort!(c); end; 119.139551 seconds (449.91 k allocations: 29.846 GiB, 5.14% gc time) import Base.Threads.@spawn const SMALL_THRESHOLD = 2^6 const BIG_THRESHOLD = 2^16 # insertionsort function isort!(v, lo::Int=1, hi::Int=length(v)) @inbounds for i = lo+1:hi j = i x = v[i] while j > lo if x < v[j-1] v[j] = v[j-1] j -= 1 continue end break end v[j] = x end return v end # single threaded merge (from t to v) function nmerge!(v, lo::Int, mid::Int, hi::Int, t) i, k, j = lo, lo, mid+1 @inbounds while i <= mid && j <= hi if t[j] < t[i] v[k] = t[j] j += 1 else v[k] = t[i] i += 1 end k += 1 end @inbounds while i <= mid v[k] = t[i] k += 1 i += 1 end @inbounds while j <= hi v[k] = t[j] k += 1 j += 1 end return v end # single threaded nocopy mergesort, one preallocation function nsort!(v, lo::Int, hi::Int, t) @inbounds if hi-lo <= SMALL_THRESHOLD return isort!(v, lo, hi) end mid = (lo+hi)>>>1 nsort!(t, lo, mid, v) nsort!(t, mid+1, hi, v) nmerge!(v, lo, mid, hi, t) return v end function nsort!(v, lo::Int=1, hi::Int=length(v)) t = copy(v) nsort!(v,lo,hi,t) return v end # multithreaded mergesort, dynamic allocation function psort!(v, lo::Int=1, hi::Int=length(v)) @inbounds if hi-lo <= BIG_THRESHOLD sort!(view(v, lo:hi), alg = MergeSort) return v end mid = (lo+hi)>>>1 # find the midpoint half = @spawn psort!(v, lo, mid) # task to sort the lower half; will run psort!(v, mid+1, hi) # in parallel with the current call sorting # the upper half wait(half) # wait for the lower half to finish temp = v[lo:mid] # workspace for merging i, k, j = 1, lo, mid+1 # merge the two sorted sub-arrays @inbounds while k < j <= hi if v[j] < temp[i] v[k] = v[j] j += 1 else v[k] = temp[i] i += 1 end k += 1 end @inbounds while k < j v[k] = temp[i] k += 1 i += 1 end return v end # mergesort, preallocation per thread, re-allocation in single threaded subroutine function tsort!(v, lo::Int=1, hi::Int=length(v), temps=[similar(v, 0) for i = 1:Threads.nthreads()]) if hi - lo < BIG_THRESHOLD # below some cutoff, run in serial sort!(view(v, lo:hi), alg = MergeSort) return v end mid = (lo+hi)>>>1 # find the midpoint half = @spawn tsort!(v, lo, mid, temps) # task to sort the lower half; will run tsort!(v, mid+1, hi, temps) # in parallel with the current call sorting # the upper half wait(half) # wait for the lower half to finish temp = temps[Threads.threadid()] # workspace for merging length(temp) < mid-lo+1 && resize!(temp, mid-lo+1) copyto!(temp, 1, v, lo, mid-lo+1) i, k, j = 1, lo, mid+1 # merge the two sorted sub-arrays @inbounds while k < j <= hi if v[j] < temp[i] v[k] = v[j] j += 1 else v[k] = temp[i] i += 1 end k += 1 end @inbounds while k < j v[k] = temp[i] k += 1 i += 1 end return v end # multithreaded nocopy mergesort, one preallocation, re-allocation in single threaded subroutine function ksort!(v, lo::Int, hi::Int, t) if hi - lo < BIG_THRESHOLD # below some cutoff, run in serial sort!(view(v, lo:hi), alg = MergeSort) return v end mid = (lo+hi)>>>1 # find the midpoint half = @spawn ksort!(t, lo, mid, v) # task to sort the lower half; will run ksort!(t, mid+1, hi, v) # in parallel with the current call sorting wait(half) # wait for the lower half to finish nmerge!(v, lo, mid, hi, t) return v end function ksort!(v, lo::Int=1, hi::Int=length(v)) t = copy(v) ksort!(v,lo,hi,t) return v end # multithreaded nocopy mergesort, one preallocation reused in single threaded subroutine function jsort!(v, lo::Int, hi::Int, t) if hi - lo < BIG_THRESHOLD # below some cutoff, run in serial return nsort!(v, lo, hi, t) end mid = (lo+hi)>>>1 # find the midpoint half = @spawn jsort!(t, lo, mid, v) # task to sort the lower half; will run jsort!(t, mid+1, hi, v) # in parallel with the current call sorting wait(half) # wait for the lower half to finish nmerge!(v, lo, mid, hi, t) return v end function jsort!(v, lo::Int=1, hi::Int=length(v)) t = copy(v) jsort!(v,lo,hi,t) return v end a = rand(2000); b = copy(a); @time sort!(b, alg = MergeSort); # reference c = copy(a); @time isort!(c); # insertionsort all(b==c) # regression test a = rand(20000000); b = copy(a); @time sort!(b, alg = MergeSort); # single-threaded system mergesort used by psort and tsort c = copy(a); @time nsort!(c); # single threaded nocopy mergesort, one preallocation all(b==c) # regression test c = copy(a); @time psort!(c); # 4 threads, mergesort, dynamic allocation all(b==c) # regression test c = copy(a); @time tsort!(c); # 4 threads, mergesort, preallocation per thread, re-allocation in single threaded subroutine all(b==c) # regression test c = copy(a); @time ksort!(c); # 4 threads, nocopy mergesort, one preallocation, re-allocation in single threaded subroutine all(b==c) # regression test c = copy(a); @time jsort!(c); # 4 threads, nocopy mergesort, one preallocation reused in single threaded subroutine all(b==c) # regression test # 4 threads, mergesort, dynamic allocation @time for i=1:100 c = copy(a); psort!(c); end; # 4 threads, mergesort, preallocation per thread, re-allocation in single threaded subroutine @time for i=1:100 c = copy(a); tsort!(c); end; # 4 threads, nocopy mergesort, one preallocation, re-allocation in single threaded subroutine @time for i=1:100 c = copy(a); ksort!(c); end; # 4 threads, nocopy mergesort, one preallocation reused in single threaded subroutine @time for i=1:100 c = copy(a); jsort!(c); end;

289 points by KenoFischer 7y ago

sjkelly 7y ago

Live stream here: https://www.youtube.com/watch?v=HfiRnfKxI64 https://www.youtube.com/watch?v=HfiRnfKxI64

eigenspace 7y ago

This is huge! I was testing out the master branch a few days ago and the parallelism improvements were amazing.

pkofod 7y ago
I think it's great that they went ahead and released a binary for people to try it out without having to compile from source.
- staticfloat 7y ago
  Thanks! We do our best to make sure that people who want to get ahold of the latest Julia binaries are always able to. Our downloads page [0] always offers the latest stable builds, whereas the nightlies page [1] allows you to download binaries off of the tip of the `master` branch, which usually lag behind the github repository by only a couple of hours. [0] https://julialang.org/downloads/ https://julialang.org/downloads/ [1] https://julialang.org/downloads/nightlies.html https://julialang.org/downloads/nightlies.html

mindB 7y ago

It's going to be really neat to watch what the Julia package ecosystem comes up with now that they have access to truly composable parallelism at the language level. Up until now, most parallelism in the language had to live at the application level rather than being in library code since parallelism either relied on spawning new processes or the Threads.@threads macro which couldn't be nested.

StefanKarpinski 7y ago

It's hard to overstate what a big deal this is (and how much effort went into it). The only mainstream languages with this kind of composable parallelism are Go and Erlang, both of which target writing servers, rather than high-performance computing, and which therefore optimize for I/O throughput rather than CPU throughput. The only systems which have supported compute-oriented composable parallelism like this are Cilk and TBB, which are both really excellent and successful—and provided major inspiration for this work. However, they both seem to have fallen short of widespread adoption partly due to being non-free (until recently) and also, I suspect, because they are non-standard C/C++ language/compiler extensions, which seem to never catch on that broadly. This is a brave new world for high-performance technical computing, not just for Julia.

eigenspace 7y ago
As a physicist, it's been really cool seeing the awareness of Julia in the community gradually increase. Large swathes of the physics community have been migrating their code to Python over the past few years but it really seems now that Julia is getting so many amazing features and packages that it can't be ignored.
- benjaminva 7y ago
  Same here, I worked in a computational biophysics department with lots of python/bash/R and I was the only one who wrote lots of high-performance code in Julia. People were curious about the language but it was still very much unknown. Hope we will see a broader adoption of Julia in the future - it's just that much better for the stuff we do on a daily basis.
Arnavion 7y ago
>The only mainstream languages with this kind of composable parallelism are Go and Erlang, both of which target writing servers, rather than high-performance computing, and which therefore optimize for I/O throughput rather than CPU throughput. Rust, while perhaps not as mainstream and not necessarily targeting scientific computing users, has data-bound parallelization like this via the rayon library.

3JPLW 7y ago

I'm super excited about this! For me the most exciting part is these four sentences: > The model is portable and free from low-level details. You don’t need to explicitly start and stop threads, and you don’t even need to know how many processors or threads there are (though you can find out if you want). > The model is nestable and composable: you can start parallel tasks that call library functions that themselves start parallel tasks, and everything works. Your CPUs will not be over-subscribed with threads. It means that you can put threaded parallelism inside libraries without worrying. It means that Julia itself will start making its builtins threaded. It means that Julia itself is much more threadsafe by default. It means that packages can implement threaded algorithms — and you can use multiple threaded packages together at the same time without worrying. And it means that you can put threaded parallelism in your application code and all these things will work together beautifully.

thekhatribharat 7y ago
Well, only if your code and your libraries do not share state.
jstrong 7y ago
still have to be careful about data races and thread synchronization, though: julia> xs = [1] 1-element Array{Int64,1}: 1 julia> for j in 1:1024 d = @spawn xs[1] += 1 end julia> xs 1-element Array{Int64,1}: 686
- deleted 7y ago
  [deleted]

polskibus 7y ago

Is this something similar to PLINQ in C# ?

ska 7y ago
IIUC it is more general, and operates at a lower level.
mr_crankypants 7y ago
It seems like it's more like C#'s PLINQ and TPL merged into a single abstraction.
manigandham 7y ago
It's basically the Task-based async model.
davidacoder 7y ago
You could use this to build something like PLINQ, but this is at a lower level in the stack than PLINQ.

cutler 7y ago

Why isn't some big company putting its weight behind Julia considering what it has to offer in the field of data science? I just don't get why it's still playing second fiddle to Python. Tech is full of inertia, paradoxically.

6thaccount2 7y ago
I think some finance groups are probably putting in some donations if I had to guess.
- cutler 7y ago
  That's not the same as Google's support for Python, Go and Kotlin or Microsoft's support for Typescript and maybe soon Rust. I would have thought the case for supporting Julia was very strong considering Python's performance and concurrency problems.
  - 6thaccount2 7y ago
    Python's performance is less of an issue when you use libraries like Numpy that are basically running optimized C code. Depending on what you're doing, Python is usually great at performance where the time it saves the user writing code is more valuable than the performance itself (a lot of exploratory scientific research).
    - DNF2 7y ago
      Numpy is great if you can vectorize your code in a straightforward manner. But when you start doing anything non-standard it becomes difficult. I can tell you that programmer productivity can plummet drastically when you need to get everything into vectorized form. Also, custom classes and types do not vectorize well at all.
      - 6thaccount2 7y ago
        Agreed.

argd678 7y ago

Does anyone have any recommendations for using Julia on real world projects at this point? Is it mature enough to use and debug on a day to day basis for example? Mainly I’ve been using R for context.

npr11 7y ago
I work at company using Julia in production. What sort of recommendations are you looking for? :)
- afiori 7y ago
  Not the parent, but I am interested in a good list of libraries that might be non obvious to find.
  - ddragon 7y ago
    You can use the search tools to find libraries for the categories you're interested at: https://juliaobserver.com/packages https://juliaobserver.com/packages https://pkg.julialang.org/docs/ https://pkg.julialang.org/docs/ As a personal answer, I really like unitful for keeping track of physical units and libraries like MLStyle.jl or Match.jl which add stuff like pattern matching to the language (though they are more general helper libraries compared to stuff like Turing.jl and DifferentialEquations.jl): https://painterqubits.github.io/Unitful.jl/stable/highlights/ https://painterqubits.github.io/Unitful.jl/stable/highlights... https://github.com/thautwarm/MLStyle.jl https://github.com/thautwarm/MLStyle.jl
- argd678 7y ago
  Mainly, am I going to spend hours fussing with the basics of the language. For example in R, it’s too slow, so I had to write Rcpp, only to discover there’s no real way to debug it from RStudio, down into the rabbit hole, lldb, Rinside, everything that can crash does, terse errors messages with out line numbers, opaque SEXP types XCode can’t decode etc. So basically are the normal everyday use cases fraught with quirks you need to know more about the internals of Julia or can you just be a user? If that makes sense.

karmakaze 7y ago

Not only is this huge for Julia, but they've just thrown down the gauntlet. The status quo has been upset. I expect Julia to start eating everyone's lunch starting with Python. Every language can use good concurrency & parallelism support and this is the biggest news for all dynamic languages.

cshenton 7y ago

Julia is already insanely productive for CS and statistics research and this will only make it more so. Super excited to see how this will be used, especially in PPLs, as algorithms like PMCMC benefit hugely from this kind of cooperative, interruptible concurrency.

joehl 7y ago

Maybe I misunderstand, but: In R allocating and copying huge return vectors for small changes on input vectors ruins performance. My understanding of Julia was that one advantage over R is the possiblity to mutate input vectors in place. Referring to the psort! example for Julia's multithreading I wonder why "We simply need to allocate one initially-empty array per thread" which to me looks like allocating huge vectors for small changes down the call recursion (and wastes memory by the number of threads). If Julia can mutate psort!s first argument global vector 'v' from multiple threads, why can't Julia also mutate a single allocated 'temp' argument global vector? Bonus: then psort! could just merge from one argument to the other and save copying back on each recursion level (along the lines of Segewick's nocopy mergesort)?

joehl 7y ago
OK, I found out myself: Julia can mutate a single temp buffer across all threads, and implementing Sedgewicks nocopy sort (with ping-pong merge) is faster, particularly when passing down preallocated buffer memory to the single threaded sorting Routine once N <= BIG_THRESHOLD. The following are measurements run with JULIA_NUM_THREADS=4 on a Intel(R) Core(TM) i5-4300M CPU @ 2.60GHz, 2594 MHz, 2 physical cores, 4 logical cores, replication code follows. Note that not only runtime of nocopy mergesort ist better, it also causes less garbage collection time. Conclusion: it would make sense for Julia's sorting interface to allow passing in preallocated buffer memory and to use nocopy mergesort behind the curtain. Of course true parallel sorting also requires a parallel merge. Best julia> # 4 threads, mergesort, dynamic allocation julia> @time for i=1:100 c = copy(a); psort!(c); end; 166.270792 seconds (756.58 k allocations: 89.470 GiB, 15.34% gc time) julia> # 4 threads, mergesort, preallocation per thread, re-allocation in single threaded subroutine julia> @time for i=1:100 c = copy(a); tsort!(c); end; 143.973485 seconds (658.00 k allocations: 36.972 GiB, 5.89% gc time) julia> # 4 threads, nocopy mergesort, one preallocation, re-allocation in single threaded subroutine julia> @time for i=1:100 c = copy(a); ksort!(c); end; 129.715467 seconds (654.51 k allocations: 37.312 GiB, 6.43% gc time) julia> # 4 threads, nocopy mergesort, one preallocation reused in single threaded subroutine julia> @time for i=1:100 c = copy(a); jsort!(c); end; 119.139551 seconds (449.91 k allocations: 29.846 GiB, 5.14% gc time) import Base.Threads.@spawn const SMALL_THRESHOLD = 2^6 const BIG_THRESHOLD = 2^16 # insertionsort function isort!(v, lo::Int=1, hi::Int=length(v)) @inbounds for i = lo+1:hi j = i x = v[i] while j > lo if x < v[j-1] v[j] = v[j-1] j -= 1 continue end break end v[j] = x end return v end # single threaded merge (from t to v) function nmerge!(v, lo::Int, mid::Int, hi::Int, t) i, k, j = lo, lo, mid+1 @inbounds while i <= mid && j <= hi if t[j] < t[i] v[k] = t[j] j += 1 else v[k] = t[i] i += 1 end k += 1 end @inbounds while i <= mid v[k] = t[i] k += 1 i += 1 end @inbounds while j <= hi v[k] = t[j] k += 1 j += 1 end return v end # single threaded nocopy mergesort, one preallocation function nsort!(v, lo::Int, hi::Int, t) @inbounds if hi-lo <= SMALL_THRESHOLD return isort!(v, lo, hi) end mid = (lo+hi)>>>1 nsort!(t, lo, mid, v) nsort!(t, mid+1, hi, v) nmerge!(v, lo, mid, hi, t) return v end function nsort!(v, lo::Int=1, hi::Int=length(v)) t = copy(v) nsort!(v,lo,hi,t) return v end # multithreaded mergesort, dynamic allocation function psort!(v, lo::Int=1, hi::Int=length(v)) @inbounds if hi-lo <= BIG_THRESHOLD sort!(view(v, lo:hi), alg = MergeSort) return v end mid = (lo+hi)>>>1 # find the midpoint half = @spawn psort!(v, lo, mid) # task to sort the lower half; will run psort!(v, mid+1, hi) # in parallel with the current call sorting # the upper half wait(half) # wait for the lower half to finish temp = v[lo:mid] # workspace for merging i, k, j = 1, lo, mid+1 # merge the two sorted sub-arrays @inbounds while k < j <= hi if v[j] < temp[i] v[k] = v[j] j += 1 else v[k] = temp[i] i += 1 end k += 1 end @inbounds while k < j v[k] = temp[i] k += 1 i += 1 end return v end # mergesort, preallocation per thread, re-allocation in single threaded subroutine function tsort!(v, lo::Int=1, hi::Int=length(v), temps=[similar(v, 0) for i = 1:Threads.nthreads()]) if hi - lo < BIG_THRESHOLD # below some cutoff, run in serial sort!(view(v, lo:hi), alg = MergeSort) return v end mid = (lo+hi)>>>1 # find the midpoint half = @spawn tsort!(v, lo, mid, temps) # task to sort the lower half; will run tsort!(v, mid+1, hi, temps) # in parallel with the current call sorting # the upper half wait(half) # wait for the lower half to finish temp = temps[Threads.threadid()] # workspace for merging length(temp) < mid-lo+1 && resize!(temp, mid-lo+1) copyto!(temp, 1, v, lo, mid-lo+1) i, k, j = 1, lo, mid+1 # merge the two sorted sub-arrays @inbounds while k < j <= hi if v[j] < temp[i] v[k] = v[j] j += 1 else v[k] = temp[i] i += 1 end k += 1 end @inbounds while k < j v[k] = temp[i] k += 1 i += 1 end return v end # multithreaded nocopy mergesort, one preallocation, re-allocation in single threaded subroutine function ksort!(v, lo::Int, hi::Int, t) if hi - lo < BIG_THRESHOLD # below some cutoff, run in serial sort!(view(v, lo:hi), alg = MergeSort) return v end mid = (lo+hi)>>>1 # find the midpoint half = @spawn ksort!(t, lo, mid, v) # task to sort the lower half; will run ksort!(t, mid+1, hi, v) # in parallel with the current call sorting wait(half) # wait for the lower half to finish nmerge!(v, lo, mid, hi, t) return v end function ksort!(v, lo::Int=1, hi::Int=length(v)) t = copy(v) ksort!(v,lo,hi,t) return v end # multithreaded nocopy mergesort, one preallocation reused in single threaded subroutine function jsort!(v, lo::Int, hi::Int, t) if hi - lo < BIG_THRESHOLD # below some cutoff, run in serial return nsort!(v, lo, hi, t) end mid = (lo+hi)>>>1 # find the midpoint half = @spawn jsort!(t, lo, mid, v) # task to sort the lower half; will run jsort!(t, mid+1, hi, v) # in parallel with the current call sorting wait(half) # wait for the lower half to finish nmerge!(v, lo, mid, hi, t) return v end function jsort!(v, lo::Int=1, hi::Int=length(v)) t = copy(v) jsort!(v,lo,hi,t) return v end a = rand(2000); b = copy(a); @time sort!(b, alg = MergeSort); # reference c = copy(a); @time isort!(c); # insertionsort all(b==c) # regression test a = rand(20000000); b = copy(a); @time sort!(b, alg = MergeSort); # single-threaded system mergesort used by psort and tsort c = copy(a); @time nsort!(c); # single threaded nocopy mergesort, one preallocation all(b==c) # regression test c = copy(a); @time psort!(c); # 4 threads, mergesort, dynamic allocation all(b==c) # regression test c = copy(a); @time tsort!(c); # 4 threads, mergesort, preallocation per thread, re-allocation in single threaded subroutine all(b==c) # regression test c = copy(a); @time ksort!(c); # 4 threads, nocopy mergesort, one preallocation, re-allocation in single threaded subroutine all(b==c) # regression test c = copy(a); @time jsort!(c); # 4 threads, nocopy mergesort, one preallocation reused in single threaded subroutine all(b==c) # regression test # 4 threads, mergesort, dynamic allocation @time for i=1:100 c = copy(a); psort!(c); end; # 4 threads, mergesort, preallocation per thread, re-allocation in single threaded subroutine @time for i=1:100 c = copy(a); tsort!(c); end; # 4 threads, nocopy mergesort, one preallocation, re-allocation in single threaded subroutine @time for i=1:100 c = copy(a); ksort!(c); end; # 4 threads, nocopy mergesort, one preallocation reused in single threaded subroutine @time for i=1:100 c = copy(a); jsort!(c); end;