The Different Flavors of Parallelism

Chris Rackauckas
September 25th, 2020

Youtube Video Link

Now that you are aware of the basics of parallel computing, let's give a high level overview of the differences between different modes of parallelism.

Lowest Level: SIMD

Recall SIMD, the idea that processors can run multiple commands simultaneously on specially structured data. "Single Instruction Multiple Data". SIMD is parallelism within a single core.

High Level Idea of SIMD

Calculations can occur in parallel in the processor if there is sufficient structure in the computation.

How to do SIMD

The simplest way to do SIMD is simply to make sure that your values are aligned. If they are, then great, LLVM's autovectorizer pass has a good chance of automatic vectorization (in the world of computing, "SIMD" is synonymous with vectorization since it is taking specific values and instead computing on small vectors. That is not to be confused with "vectorization" in the sense of Python/R/MATLAB, which is a programming style which prefers using C-defined primitive functions, like broadcast or matrix multiplication).

You can check for auto-vectorization inside of the LLVM IR by looking for statements like:

%wide.load24 = load <4 x double>, <4 x double> addrspac(13)* %46, align 8
; └
; ┌ @ float.jl:395 within `+'
%47 = fadd <4 x double> %wide.load, %wide.load24

which means that 4 additions are happening simultaneously. The amount of vectorization is heavily dependent on your architecture. The ancient form of SIMD, the SSE(2) instructions, required that your data was aligned. Now there's a bit more leeway, but generally it holds that making your the data you're trying to SIMD over is aligned. Thus there can be major differences in computing using a struct of array format instead of an arrays of structs format. For example:

struct MyComplex
  real::Float64
  imag::Float64
end
arr = [MyComplex(rand(),rand()) for i in 1:100]
100-element Vector{MyComplex}:
 MyComplex(0.7302952470177198, 0.5865425609737499)
 MyComplex(0.22866594502824422, 0.7378546901164249)
 MyComplex(0.8545732103550867, 0.7431326934659985)
 MyComplex(0.2871537808843082, 0.7580834228260477)
 MyComplex(0.6044578113288847, 0.007172996524091979)
 MyComplex(0.27572962333918005, 0.03893669055526128)
 MyComplex(0.8133831421086706, 0.7110612365510052)
 MyComplex(0.7617371057392228, 0.30813709978118264)
 MyComplex(0.6587680508222387, 0.030382505581488917)
 MyComplex(0.500480663139846, 0.06820104396111737)
 ⋮
 MyComplex(0.040466707682546565, 0.9592861297724803)
 MyComplex(0.8387114285911098, 0.7463524676310868)
 MyComplex(0.1639290039029072, 0.7141651059352342)
 MyComplex(0.9322836528887831, 0.5648163106996179)
 MyComplex(0.9461734646601937, 0.5509146910909132)
 MyComplex(0.7572780881492088, 0.27782447742514793)
 MyComplex(0.479465567775742, 0.5814379735637724)
 MyComplex(0.8863621829174242, 0.9056593573420595)
 MyComplex(0.931213703369945, 0.09972313071599537)

is represented in memory as

[real1,imag1,real2,imag2,...]

while the struct of array formats are

struct MyComplexes
  real::Vector{Float64}
  imag::Vector{Float64}
end
arr2 = MyComplexes(rand(100),rand(100))
MyComplexes([0.3378185834559957, 0.10072193483180802, 0.03052528185770631, 
0.7874008387061221, 0.8427793089584803, 0.21800852162762197, 0.614477852179
1766, 0.7445108139534009, 0.9946095702710185, 0.3040193278673452  …  0.6301
836561621345, 0.6793943354709119, 0.12497074639200012, 0.41131900853849546,
 0.5608653969521878, 0.4034015621989776, 0.8246826300742593, 0.542089098162
4119, 0.20232275946334655, 0.0980837773276041], [0.7468032896731819, 0.9146
244130621338, 0.8254522471312051, 0.601679984082818, 0.8762626623289597, 0.
310851095350141, 0.8032829249344134, 0.08964167966088754, 0.914082483507435
, 0.6854698246317493  …  0.5909511213082141, 0.585719488768233, 0.045049907
49236726, 0.9512096468283845, 0.7437302301075415, 0.143128801870137, 0.7883
607026920928, 0.38952521013749264, 0.6537181918694862, 0.34958699676248106]
)

Now let's check what happens when we perform a reduction:

using InteractiveUtils
Base.:+(x::MyComplex,y::MyComplex) = MyComplex(x.real+y.real,x.imag+y.imag)
Base.:/(x::MyComplex,y::Int) = MyComplex(x.real/y,x.imag/y)
average(x::Vector{MyComplex}) = sum(x)/length(x)
@code_llvm average(arr)
;  @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/styles_of_paral
lelism.jmd:5 within `average`
define void @julia_average_16495([2 x double]* noalias nocapture sret([2 x 
double]) %0, {}* nonnull align 16 dereferenceable(40) %1) #0 {
top:
  %2 = alloca [4 x {}*], align 8
  %3 = alloca <2 x double>, align 16
  %tmpcast = bitcast <2 x double>* %3 to [2 x double]*
  %.sub = getelementptr inbounds [4 x {}*], [4 x {}*]* %2, i64 0, i64 0
; ┌ @ reducedim.jl:889 within `sum`
; │┌ @ reducedim.jl:889 within `#sum#738`
; ││┌ @ reducedim.jl:893 within `_sum`
; │││┌ @ reducedim.jl:893 within `#_sum#740`
; ││││┌ @ reducedim.jl:894 within `_sum`
; │││││┌ @ reducedim.jl:894 within `#_sum#741`
; ││││││┌ @ reducedim.jl:322 within `mapreduce`
; │││││││┌ @ reducedim.jl:322 within `#mapreduce#731`
; ││││││││┌ @ reducedim.jl:330 within `_mapreduce_dim`
; │││││││││┌ @ reduce.jl:399 within `_mapreduce`
; ││││││││││┌ @ indices.jl:457 within `LinearIndices`
; │││││││││││┌ @ abstractarray.jl:95 within `axes`
; ││││││││││││┌ @ array.jl:151 within `size`
               %4 = bitcast {}* %1 to {}**
               %5 = getelementptr inbounds {}*, {}** %4, i64 3
               %6 = bitcast {}** %5 to i64*
               %7 = load i64, i64* %6, align 8
; │││││││││└└└└
; │││││││││┌ @ reduce.jl:401 within `_mapreduce`
            switch i64 %7, label %L14 [
    i64 0, label %L8
    i64 1, label %L12
  ]

L8:                                               ; preds = %top
; │││││││││└
; │││││││││┌ @ reduce.jl:402 within `_mapreduce`
            store {}* inttoptr (i64 140610866599904 to {}*), {}** %.sub, al
ign 8
            %8 = getelementptr inbounds [4 x {}*], [4 x {}*]* %2, i64 0, i6
4 1
            store {}* inttoptr (i64 140610888299072 to {}*), {}** %8, align
 8
            %9 = getelementptr inbounds [4 x {}*], [4 x {}*]* %2, i64 0, i6
4 2
            store {}* %1, {}** %9, align 8
            %10 = getelementptr inbounds [4 x {}*], [4 x {}*]* %2, i64 0, i
64 3
            store {}* inttoptr (i64 140610855818704 to {}*), {}** %10, alig
n 8
            %11 = call nonnull {}* @jl_invoke({}* inttoptr (i64 14061086838
2144 to {}*), {}** nonnull %.sub, i32 4, {}* inttoptr (i64 140607983488624 
to {}*))
            call void @llvm.trap()
            unreachable

L12:                                              ; preds = %top
; │││││││││└
; │││││││││┌ @ reduce.jl:404 within `_mapreduce`
; ││││││││││┌ @ array.jl:861 within `getindex`
             %12 = bitcast {}* %1 to <2 x double>**
             %13 = load <2 x double>*, <2 x double>** %12, align 8
             %14 = load <2 x double>, <2 x double>* %13, align 8
; │││││││││└└
; │││││││││┌ @ reduce.jl:405 within `_mapreduce`
            br label %L46

L14:                                              ; preds = %top
; │││││││││└
; │││││││││┌ @ reduce.jl:406 within `_mapreduce`
; ││││││││││┌ @ int.jl:83 within `<`
             %15 = icmp ugt i64 %7, 15
; ││││││││││└
            br i1 %15, label %L42, label %L16

L16:                                              ; preds = %L14
; │││││││││└
; │││││││││┌ @ reduce.jl:408 within `_mapreduce`
; ││││││││││┌ @ array.jl:861 within `getindex`
             %16 = bitcast {}* %1 to [2 x double]**
             %17 = load [2 x double]*, [2 x double]** %16, align 8
             %18 = bitcast [2 x double]* %17 to <2 x double>*
             %19 = load <2 x double>, <2 x double>* %18, align 8
; │││││││││└└
; │││││││││┌ @ reduce.jl:409 within `_mapreduce`
; ││││││││││┌ @ array.jl:861 within `getindex`
             %.elt65 = getelementptr inbounds [2 x double], [2 x double]* %
17, i64 1, i64 0
             %20 = bitcast double* %.elt65 to <2 x double>*
             %21 = load <2 x double>, <2 x double>* %20, align 8
; │││││││││└└
; │││││││││┌ @ reduce.jl:410 within `_mapreduce`
; ││││││││││┌ @ reduce.jl:24 within `add_sum`
; │││││││││││┌ @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/sty
les_of_parallelism.jmd:3 within `+` @ float.jl:399
              %22 = fadd <2 x double> %19, %21
; │││││││││└└└
; │││││││││┌ @ reduce.jl:411 within `_mapreduce`
; ││││││││││┌ @ int.jl:83 within `<`
             %.not6974 = icmp ugt i64 %7, 2
; ││││││││││└
            br i1 %.not6974, label %L33, label %L46

L33:                                              ; preds = %L33, %L16
            %value_phi277 = phi i64 [ %24, %L33 ], [ 2, %L16 ]
            %23 = phi <2 x double> [ %27, %L33 ], [ %22, %L16 ]
; │││││││││└
; │││││││││┌ @ reduce.jl:412 within `_mapreduce`
; ││││││││││┌ @ int.jl:87 within `+`
             %24 = add nuw nsw i64 %value_phi277, 1
; ││││││││││└
; ││││││││││┌ @ array.jl:861 within `getindex`
             %.elt70 = getelementptr inbounds [2 x double], [2 x double]* %
17, i64 %value_phi277, i64 0
             %25 = bitcast double* %.elt70 to <2 x double>*
             %26 = load <2 x double>, <2 x double>* %25, align 8
; │││││││││└└
; │││││││││┌ @ reduce.jl:413 within `_mapreduce`
; ││││││││││┌ @ reduce.jl:24 within `add_sum`
; │││││││││││┌ @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/sty
les_of_parallelism.jmd:3 within `+` @ float.jl:399
              %27 = fadd <2 x double> %23, %26
; │││││││││└└└
; │││││││││┌ @ reduce.jl:411 within `_mapreduce`
; ││││││││││┌ @ int.jl:83 within `<`
             %exitcond.not = icmp eq i64 %24, %7
; ││││││││││└
            br i1 %exitcond.not, label %L46, label %L33

L42:                                              ; preds = %L14
; │││││││││└
; │││││││││┌ @ reduce.jl:417 within `_mapreduce`
; ││││││││││┌ @ reduce.jl:259 within `mapreduce_impl`
             call void @j_mapreduce_impl_16497([2 x double]* noalias nocapt
ure nonnull sret([2 x double]) %tmpcast, {}* nonnull %1, i64 signext 1, i64
 signext %7, i64 signext 1024) #0
; └└└└└└└└└└└
  %28 = load <2 x double>, <2 x double>* %3, align 16
; ┌ @ reducedim.jl:889 within `sum`
; │┌ @ reducedim.jl:889 within `#sum#738`
; ││┌ @ reducedim.jl:893 within `_sum`
; │││┌ @ reducedim.jl:893 within `#_sum#740`
; ││││┌ @ reducedim.jl:894 within `_sum`
; │││││┌ @ reducedim.jl:894 within `#_sum#741`
; ││││││┌ @ reducedim.jl:322 within `mapreduce`
; │││││││┌ @ reducedim.jl:322 within `#mapreduce#731`
; ││││││││┌ @ reducedim.jl:330 within `_mapreduce_dim`
; │││││││││┌ @ reduce.jl:417 within `_mapreduce`
            br label %L46

L46:                                              ; preds = %L42, %L33, %L1
6, %L12
; └└└└└└└└└└
  %29 = phi <2 x double> [ %28, %L42 ], [ %14, %L12 ], [ %22, %L16 ], [ %27
, %L33 ]
; ┌ @ array.jl:215 within `length`
   %30 = bitcast {}* %1 to { i8*, i64, i16, i16, i32 }*
   %31 = getelementptr inbounds { i8*, i64, i16, i16, i32 }, { i8*, i64, i1
6, i16, i32 }* %30, i64 0, i32 1
   %32 = load i64, i64* %31, align 8
; └
; ┌ @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/styles_of_para
llelism.jmd:4 within `/` @ promotion.jl:382
; │┌ @ promotion.jl:350 within `promote`
; ││┌ @ promotion.jl:327 within `_promote`
; │││┌ @ number.jl:7 within `convert`
; ││││┌ @ float.jl:146 within `Float64`
       %33 = sitofp i64 %32 to double
; │└└└└
; │ @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/styles_of_para
llelism.jmd:4 within `/` @ promotion.jl:382 @ float.jl:408
   %34 = insertelement <2 x double> poison, double %33, i32 0
   %35 = shufflevector <2 x double> %34, <2 x double> undef, <2 x i32> zero
initializer
   %36 = fdiv <2 x double> %29, %35
; └
  %37 = bitcast [2 x double]* %0 to <2 x double>*
  store <2 x double> %36, <2 x double>* %37, align 8
  ret void
}

What this is doing is creating small little vectors and then parallelizing the operations of those vectors by calling specific vector-parallel instructions. Keep this in mind.

Explicit SIMD

The following was all a form of loop-level parallelism known as loop vectorization. It's simply easier for compilers to reason at the array level, prove iterates are independent, and automatically generate SIMD code from that. This is not necessary, and compilers can produce SIMD code from non-looping code through a process known as SLP supervectorization, but the results are far from optimal and the compiler requires a lot of time to do this calculation, meaning that it's usually not a pass used by default.

If you want to pack the vectors yourself, then primitives for doing so from within Julia are available in SIMD.jl. This is for "real" performance warriors. This looks like for example:

using SIMD
v = Vec{4,Float64}((1,2,3,4))
@show v+v # basic arithmetic is supported
@show sum(v) # basic reductions are supported
v + v = <4 x Float64>[2.0, 4.0, 6.0, 8.0]
sum(v) = 10.0
10.0

Using this you can pull apart code and force the usage of SIMD vectors. One library which makes great use of this is LoopVectorization.jl. However, one word of "caution":

Most performance optimization is not trying to do something really good for performance. Most performance optimization is trying to not do something that is actively bad for performance.

Summary of SIMD

  • Communication in SIMD is due to locality: if things are local the processor can automatically setup the operations.

  • There's no real worry about "getting it wrong": you cannot overwrite pieces from different parts of the arithmetic unit, and if SIMD is unsafe then it just won't auto-vectorize.

  • Suitable for operations measured in ns.

Next Level Up: Multithreading

Last time we briefly went over multithreading and described how every process has multiple threads which share a single heap, and when multiple threads are executed simultaneously we have multithreaded parallelism. Note that you can have multiple threads which aren't executed simultaneously, like in the case of I/O operations, and this is an example of concurrency without parallelism and is commonly referred to as green threads.

Last time we described a simple multithreaded program and noticed that multithreading has an overhead cost of around 50ns-100ns. This is due to the construction of the new stack (among other things) each time a new computational thread is spun up. This means that, unlike SIMD, some thought needs to be put in as to when to perform multithreading: it's not always a good idea. It needs to be high enough on the cost for this to be counter-balanced.

One abstraction that was glossed over was the memory access style. Before, we were considering a single heap, or an UMA style:

However, this is the case for all shared memory devices. For example, compute nodes on the HPC tend to be "dual Xeon" or "quad Xeon", where each Xeon processor is itself a multicore processor. But each processor on its own accesses its own local caches, and thus one has to be aware that this is setup in a NUMA (non-uniform memory access) manner:

where there is a cache that is closer to the processor and a cache that is further away. Care should be taken in this to localize the computation per thread, otherwise a cost associated with the memory sharing will be hit (but all sharing will still be automatic).

In this sense, interthread communication is naturally done through the heap: if you want other threads to be able to touch a value, then you can simply place it on the heap and then it'll be available. We saw this last time by how overlapping computations can re-use the same heap-based caches, meaning that care needs to be taken with how one writes into a dynamically-allocated array.

A simple example that demonstrates this is. First, let's make sure we have multithreading enabled:

using Base.Threads
Threads.nthreads() # should not be 1
1
using BenchmarkTools
acc = 0
@threads for i in 1:10_000
    global acc
    acc += 1
end
acc
10000

The reason for this behavior is that there is a difference between the reading and the writing step to an array. Here, values are being read while other threads are writing, meaning that they see a lower value than when they are attempting to write into it. The result is that the total summation is lower than the true value because of this clashing. We can prevent this by only allowing one thread to utilize the heap-allocated variable at a time. One abstraction for doing this is atomics:

acc = Atomic{Int64}(0)
@threads for i in 1:10_000
    atomic_add!(acc, 1)
end
acc
Atomic{Int64}(10000)

When an atomic add is being done, all other threads wishing to do the same computation are blocked. This of course can have a massive effect on performance since atomic computations are not parallel.

Julia also exposes a lower level of heap control in threading using locks

const acc_lock = Ref{Int64}(0)
const splock = SpinLock()
function f1()
    @threads for i in 1:10_000
        lock(splock)
        acc_lock[] += 1
        unlock(splock)
    end
end
const rsplock = ReentrantLock()
function f2()
    @threads for i in 1:10_000
        lock(rsplock)
        acc_lock[] += 1
        unlock(rsplock)
    end
end
acc2 = Atomic{Int64}(0)
function g()
  @threads for i in 1:10_000
      atomic_add!(acc2, 1)
  end
end
const acc_s = Ref{Int64}(0)
function h()
  global acc_s
  for i in 1:10_000
      acc_s[] += 1
  end
end
@btime f1()
162.000 μs (6 allocations: 544 bytes)

SpinLock is non-reentrant, i.e. it will block itself if a thread that calls a lock does another lock. Therefore it has to be used with caution (every lock goes with one unlock), but it's fast. ReentrantLock alleviates those concerns, but trades off a bit of performance:

@btime f2()
482.299 μs (6 allocations: 544 bytes)

But if you can use atomics, they will be faster:

@btime g()
64.900 μs (6 allocations: 544 bytes)

and if your computation is actually serial, then use serial code:

@btime h()
4.457 μs (0 allocations: 0 bytes)

Why is this so fast? Check the code:

@code_llvm h()
;  @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/styles_of_paral
lelism.jmd:26 within `h`
define void @julia_h_16679() #0 {
top:
;  @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/styles_of_paral
lelism.jmd:29 within `h`
; ┌ @ refvalue.jl:56 within `getindex`
; │┌ @ Base.jl:42 within `getproperty`
    %.pre = load i64, i64* inttoptr (i64 140609009595344 to i64*), align 16
; └└
;  @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/styles_of_paral
lelism.jmd:28 within `h`
  br label %L2

L2:                                               ; preds = %L2, %top
;  @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/styles_of_paral
lelism.jmd:29 within `h`
; ┌ @ refvalue.jl:56 within `getindex`
; │┌ @ Base.jl:42 within `getproperty`
    %0 = phi i64 [ %.pre, %top ], [ %1, %L2 ]
    %value_phi = phi i64 [ 1, %top ], [ %2, %L2 ]
; └└
; ┌ @ int.jl:87 within `+`
   %1 = add i64 %0, 1
; └
; ┌ @ refvalue.jl:57 within `setindex!`
; │┌ @ Base.jl:43 within `setproperty!`
    store i64 %1, i64* inttoptr (i64 140609009595344 to i64*), align 16
; └└
; ┌ @ range.jl:837 within `iterate`
; │┌ @ promotion.jl:468 within `==`
    %.not = icmp eq i64 %value_phi, 10000
; │└
   %2 = add nuw nsw i64 %value_phi, 1
; └
  br i1 %.not, label %L19, label %L2

L19:                                              ; preds = %L2
  ret void
}

It just knows to add 10000. So to get a proper timing let's make the size mutable:

const len = Ref{Int}(10_000)
function h2()
  global acc_s
  global len
  for i in 1:len[]
      acc_s[] += 1
  end
end
@btime h2()
3.237 μs (0 allocations: 0 bytes)
@code_llvm h2()
;  @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/styles_of_paral
lelism.jmd:3 within `h2`
define void @julia_h2_16687() #0 {
top:
;  @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/styles_of_paral
lelism.jmd:6 within `h2`
; ┌ @ refvalue.jl:56 within `getindex`
; │┌ @ Base.jl:42 within `getproperty`
    %0 = load i64, i64* inttoptr (i64 140608011976960 to i64*), align 256
; └└
; ┌ @ range.jl:5 within `Colon`
; │┌ @ range.jl:354 within `UnitRange`
; ││┌ @ range.jl:359 within `unitrange_last`
     %.inv = icmp sgt i64 %0, 0
     %1 = select i1 %.inv, i64 %0, i64 0
; └└└
  br i1 %.inv, label %L14.preheader, label %L31

L14.preheader:                                    ; preds = %top
;  @ /home/runner/work/SciMLBook/SciMLBook/_weave/lecture06/styles_of_paral
lelism.jmd:7 within `h2`
; ┌ @ refvalue.jl:56 within `getindex`
; │┌ @ Base.jl:42 within `getproperty`
    %.pre = load i64, i64* inttoptr (i64 140609009595344 to i64*), align 16
; └└
  br label %L14

L14:                                              ; preds = %L14, %L14.preh
eader
; ┌ @ refvalue.jl:56 within `getindex`
; │┌ @ Base.jl:42 within `getproperty`
    %2 = phi i64 [ %3, %L14 ], [ %.pre, %L14.preheader ]
    %value_phi2 = phi i64 [ %4, %L14 ], [ 1, %L14.preheader ]
; └└
; ┌ @ int.jl:87 within `+`
   %3 = add i64 %2, 1
; └
; ┌ @ refvalue.jl:57 within `setindex!`
; │┌ @ Base.jl:43 within `setproperty!`
    store i64 %3, i64* inttoptr (i64 140609009595344 to i64*), align 16
; └└
; ┌ @ range.jl:837 within `iterate`
; │┌ @ promotion.jl:468 within `==`
    %.not = icmp eq i64 %value_phi2, %1
; │└
   %4 = add nuw i64 %value_phi2, 1
; └
  br i1 %.not, label %L31, label %L14

L31:                                              ; preds = %L14, %top
  ret void
}

It's still optimizing it!

non_const_len = 10000
function h3()
  global acc_s
  global non_const_len
  len2::Int = non_const_len
  for i in 1:len2
      acc_s[] += 1
  end
end
@btime h3()
3.375 μs (0 allocations: 0 bytes)

Note that what is shown here is a type-declaration. a::T = ... forces a to be of type T throughout the whole function. By giving the compiler this information, I am able to use the non-constant global in a type-stable manner.

One last thing to note about multithreaded computations, and parallel computations, is that one cannot assume that the parallelized computation is computed in any given order. For example, the following will has a quasi-random ordering:

const a2 = zeros(nthreads()*10)
const acc_lock2 = Ref{Int64}(0)
const splock2 = SpinLock()
function f_order()
    @threads for i in 1:length(a2)
        lock(splock2)
        acc_lock2[] += 1
        a2[i] = acc_lock2[]
        unlock(splock2)
    end
end
f_order()
a2
10-element Vector{Float64}:
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  8.0
  9.0
 10.0

Note that here we can see that Julia 1.5 is dividing up the work into groups of 10 for each thread, and then one thread dominates the computation at a time, but which thread dominates is random.

The Dining Philosophers Problem

A classic tale in parallel computing is the dining philosophers problem. In this case, there are N philosophers at a table who all want to eat at the same time, following all of the same rules. Each philosopher must alternatively think and then eat. They need both their left and right fork to start eating, but cannot start eating until they have both forks. The problem is how to setup a concurrent algorithm that will not cause any philosophers to starve.

The difficulty is a situation known as deadlock. For example, if each philosopher was told to grab the right fork when it's available, and then the left fork, and put down the fork after eating, then they will all grab the right fork and none will ever eat because they will all be waiting on the left fork. This is analogous to two blocked computations which are waiting on the other to finish. Thus, when using blocking structures, one needs to be careful about deadlock!

Two Programming Models: Loop-Level Parallelism and Task-Based Parallelism

As described in the previous lecture, one can also use Threads.@spawn to do multithreading in Julia v1.3+. The same factors all applay: how to do locks and Mutex etc. This is a case of a parallelism construct having two alternative programming models. Threads.@spawn represents task-based parallelism, while Threads.@threads represents Loop-Level Parallelism or a parallel iterator model. Loop-based parallelization models are very high level and, assuming every iteration is independent, almost requires no code change. Task-based parallelism is a more expressive parallelism model, but usually requires modifying the code to be explicitly written as a set of parallelizable tasks. Note that in the case of Julia, Threads.@threads is implemented using Threads.@spawn's model.

Summary of Multithreading

  • Communication in multithreading is done on the heap. Locks and atomics allow for a form of safe message passing.

  • 50ns-100ns of overhead. Suitable for 1μs calculations.

  • Be careful of ordering and heap-allocated values.

GPU Computing

GPUs are not fast. In fact, the problem with GPUs is that each processor is slow. However, GPUs have a lot of cores... like thousands.

An RTX2080, a standard "gaming" GPU (not even the ones in the cluster), has 2944 cores. However, not only are GPUs slow, but they also need to be programmed in a style that is SPMD, which standard for Single Program Multiple Data. This means that every single thread must be running the same program but on different pieces of data. Exactly the same program. If you have

if a > 1
  # Do something
else
  # Do something else
end

where some of the data goes on one branch and other data goes on the other branch, every single thread will run both branches (performing "fake" computations while on the other branch). This means that GPU tasks should be "very parallel" with as few conditionals as possible.

GPU Memory

GPUs themselves are shared memory devices, meaning they have a heap that is shared amongst all threads. However, GPUs are heavily in the NUMA camp, where different blocks of the GPU have much faster access to certain parts of the memory. Additionally, this heap is disconnected from the standard processor, so data must be passed to the GPU and data must be returned.

GPU memory size is relatively small compared to CPUs. Example: the RTX2080Ti has 8GB of RAM. Thus one needs to be doing computations that are memory compact (such as matrix multiplications, which are O(n^3) making the computation time scale quicker than the memory cost).

Note on GPU Hardware

Standard GPU hardware "for gaming", like RTX2070, is just as fast as higher end GPU hardware for Float32. Higher end hardware, like the Tesla, add more memory, memory safety, and Float64 support. However, these require being in a server since they have alternative cooling strategies, making them a higher end product.

SPMD Kernel Generation GPU Computing Models

The core programming models for GPU computing are SPMD kernel compilers, of which the most well-known is CUDA. CUDA is a C++-like programming language which compiles to .ptx kernels, and GPU execution on NVIDIA GPUs is done by "all steams" of a GPU doing concurrent execution of the kernel (generally, without going into more details, you can of "all streams" as just meaning "all cores". More detailed views of GPU execution will come later).

.ptx CUDA kernels can be compiled from LLVM IR, and thus since Julia is a programming language which emits LLVM IR for all of its operations, native Julia programs are compatible with compilation to CUDA. The helper functions to enable this separate compilation path is CUDA.jl. Let's take a look at a basic CUDA.jl kernel generating example:

using CUDA

N = 2^20
x_d = CUDA.fill(1.0f0, N)  # a vector stored on the GPU filled with 1.0 (Float32)
y_d = CUDA.fill(2.0f0, N)  # a vector stored on the GPU filled with 2.0

function gpu_add2!(y, x)
    index = threadIdx().x    # this example only requires linear indexing, so just use `x`
    stride = blockDim().x
    for i = index:stride:length(y)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y_d, 2)
@cuda threads=256 gpu_add2!(y_d, x_d)
all(Array(y_d) .== 3.0f0)
ERROR: Could not find the CUDA driver library. Please make sure you have installed the NVIDIA driver for your GPU.
If you're sure it's installed, look for `libcuda.so` in your system and make sure it's discoverable by the linker.
Typically, that involves an entry in '/etc/ld.so.conf', or setting LD_LIBRARY_PATH.

The key to understanding the SPMD kernel approach is the index = threadIdx().x and stride = blockDim().x portions.

The way kernels are expected to run in parallel is that they are given a specific block of the computation and are expected to write a kernel which only on that small block of the input. This kernel is then called on every separate thread on the GPU, making each CUDA core simultaneously compute each block. Thus as a user in such a SPMD programming model, you never specify the computation globally but instead simply specify how chunks should behave, giving the compiler the leeway to determine the optimal global execution.

Array-Based GPU Computing Models

The simplest version of GPU computing is the array-based programming model.

A = rand(100,100); B = rand(100,100)
using CUDA
# Pass to the GPU
cuA = cu(A); cuB = cu(B)
cuC = cuA*cuB
# Pass to the CPU
C = Array(cuC)
ERROR: Could not find the CUDA driver library. Please make sure you have installed the NVIDIA driver for your GPU.
If you're sure it's installed, look for `libcuda.so` in your system and make sure it's discoverable by the linker.
Typically, that involves an entry in '/etc/ld.so.conf', or setting LD_LIBRARY_PATH.

Let's see the transfer times:

@btime cu(A)
ERROR: Could not find the CUDA driver library. Please make sure you have installed the NVIDIA driver for your GPU.
If you're sure it's installed, look for `libcuda.so` in your system and make sure it's discoverable by the linker.
Typically, that involves an entry in '/etc/ld.so.conf', or setting LD_LIBRARY_PATH.
@btime Array(cuC)
ERROR: UndefVarError: cuC not defined

The cost transferring is about 20μs-50μs in each direction, meaning that one needs to be doing operations that cost at least 200μs for GPUs to break even. A good rule of thumb is that GPU computations should take at least a millisecond, or GPU memory should be re-used.

Summary of GPUs

  • GPUs cores are slow

  • GPUs are SPMD

  • GPUs are generally used for linear algebra

  • Suitable for SPMD 1ms computations

Xeon Phi Accelerators and OpenCL

Other architectures exist to keep in mind. Xeon Phis are a now-defunct accelerator that used X86 (standard processors) as the base, using hundreds of them. For example, the Knights Landing series had 256 core accelerator cards. These were all clocked down, meaning they were still slower than a standard CPU, but there were less restrictions on SPMD (though SPMD-like computations were still preferred in order to heavily make use of SIMD). However, because machine learning essentially only needs linear algebra, and linear algebra is faster when restricting to SPMD-architectures, this failed. These devices can still be found on many high end clusters.

One alternative to CUDA is OpenCL which supports alternative architectures such as the Xeon Phi at the same time that it supports GPUs. However, one of the issues with OpenCL is that its BLAS implementation currently does not match the speed of CuBLAS, which makes NVIDIA-specific libraries still the king of machine learning and most scientific computing.

TPU Computing

TPUs are tensor processing units, which is Google's newest accelerator technology. They are essentially just "tensor operation compilers", which in computer science speak is simply higher dimensional linear algebra. To do this, they internally utilize a BFloat16 type, which is a 16-bit floating point number with the same exponent size as a Float32 with an 8-bit significant. This means that computations are highly prone to catastrophic cancellation. This computational device only works because BFloat16 has primitive operations for FMA which allows 32-bit-like accuracy of multiply-add operations, and thus computations which are only dot products (linear algebra) end up okay. Thus this is simply a GPU-like device which has gone further to completely specialize in linear algebra.

Multiprocessing (Distributed Computing)

While multithreading computes with multiple threads, multiprocessing computes with multiple independent processes. Note that processes do not share any memory, not heap or data, and thus this mode of computing also allows for distributed computations, which is the case where processes may be on separate computing hardware. However, even if they are on the same hardware, the lack of a shared address space means that multiprocessing has to do message passing, i.e. send data from one process to the other.

Distributed Tasks with Explicit Memory Handling: The Master-Worker Model

Given the amount of control over data handling, there are many different models for distributed computing. The simplest, the one that Julia's Distributed Standard Library defaults to, is the master-worker model. The master-worker model has one process, deemed the master, which controls the worker processes.

Here we can start by adding some new worker processes:

using Distributed
addprocs(4)

This adds 4 worker processes for the master to control. The simplest computations are those where the master process gives the worker process a job which returns the value afterwards. For example, a pmap operation or @distributed loop gives the worker a function to execute, along with the data, and the worker then computes and returns the result.

At a lower level, this is done by Distributed.@spawning jobs, or using a remotecall and fetching the result. ParallelDataTransfer.jl gives an extended set of primitive message passing operations. For example, we can explicitly tell it to compute a function f on the remote process like:

@everywhere f(x) = x.^2 # Define this function on all processes
t = remotecall(f,2,randn(10))

remotecall is a non-blocking operation that returns a Future. To access the data, one should use the blocking operation fetch to receive the data:

xsq = fetch(t)

Distributed Tasks with Implicit Memory Handling: Distributed Task-Based Parallelism

Another popular programming model for distributed computation is task-based parallelism but where all of the memory handling is implicit. Since, unlike the shared memory parallelism case, data transfers are required for given processes to share heap allocated values, distributed task-based parallelism libraries tend to want a global view of the whole computation in order to build a sophisticated schedule that includes where certain data lives and when transfers will occur. Because of this, distributed task-based parallelism libraries tend to want the entire computational graph of the computation, to be able to restructure the graph as necessary with their own data transfer portions spliced into the compute. Examples of this kind of framework are:

  • Tensorflow

  • dask ("distributed tasks")

  • Dagger.jl

Using these kinds of libraries requires building a directed acyclic graph (DAG). For example, the following showcases how to use Dagger.jl to represent a bunch of summations:

using Dagger

add1(value) = value + 1
add2(value) = value + 2
combine(a...) = sum(a)

p = delayed(add1)(4)
q = delayed(add2)(p)
r = delayed(add1)(3)
s = delayed(combine)(p, q, r)

@assert collect(s) == 16

Once the global computation is specified, commands like collect are used to instantiate the graph on given input data, which then run the computation in a (potentially) distributed manner, depending on internal scheduler heuristics.

Distributed Array-Based Parallelism: SharedArrays, Elemental, and DArrays

Because array operations are a standard way to compute in scientific computing, there are higher level primitives to help with message passing. A SharedArray is an array which acts like a shared memory device. This means that every change to a SharedArray causes message passing to keep them in sync, and thus this should be used with a performance caution. DistributedArrays.jl is a parallel array type which has local blocks and can be used for writing higher level abstractions with explicit message passing. Because it is currently missing high-level parallel linear algebra, currently the recommended tool for distributed linear algebra is Elemental.jl.

MapReduce, Hadoop, and Spark: The Map-Reduce Model

Many data-parallel operations work by mapping a function f onto each piece of data and then reducing it. For example, the sum of squares maps the function x -> x^2 onto each value, and then these values are reduced by performing a summation. MapReduce was a Google framework in the 2000's built around this as the parallel computing concept, and current data-handling frameworks, like Hadoop and Spark, continue this as the core distributed programming model.

In Julia, there exists the mapreduce function for performing serial mapreduce operations. It also work on GPUs. However, it does not auto-distribute. For distributed map-reduce programming, the @distributed for-loop macro can be used. For example, sum of squares of random numbers is:

@distributed (+) for i in 1:1000
  rand()^2
end

One can see that computing summary statistics is easily done in this framework which is why it was majorly adopted among "big data" communities.

@distributed uses a static scheduler. The dynamic scheduling equivalent is pmap:

pmap(i->rand()^2,1:100)

which will dynamically allocate jobs to processes as they declare they have finished jobs. This thus has the same performance difference behavior as Threads.@threads vs Threads.@spawn.

MPI: The Distributed SPMD Model

The main way to do high-performance multiprocessing is MPI, which is an old distributed computing interface from the C/Fortran days. Julia has access to the MPI programming model through MPI.jl. The programming model for MPI is that every computer is running the same program, and synchronization is performed by blocking communication. For example, let's look at the following:

using MPI
MPI.Init()

comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
size = MPI.Comm_size(comm)

dst = mod(rank+1, size)
src = mod(rank-1, size)

N = 4

send_mesg = Array{Float64}(undef, N)
recv_mesg = Array{Float64}(undef, N)

fill!(send_mesg, Float64(rank))

rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)

print("$rank: Sending   $rank -> $dst = $send_mesg\n")
sreq = MPI.Isend(send_mesg, dst, rank+32, comm)

stats = MPI.Waitall!([rreq, sreq])

print("$rank: Received $src -> $rank = $recv_mesg\n")

MPI.Barrier(comm)
> mpiexecjl -n 3 julia examples/04-sendrecv.jl
1: Sending   1 -> 2 = [1.0, 1.0, 1.0, 1.0]
0: Sending   0 -> 1 = [0.0, 0.0, 0.0, 0.0]
1: Received 0 -> 1 = [0.0, 0.0, 0.0, 0.0]
2: Sending   2 -> 0 = [2.0, 2.0, 2.0, 2.0]
0: Received 2 -> 0 = [2.0, 2.0, 2.0, 2.0]
2: Received 1 -> 2 = [1.0, 1.0, 1.0, 1.0]

Let's investigate this a little bit. Think about having two computers run this line-by-line side by side. They will both locally build arrays, and then call MPI.Irecv!, which is an asynchronous non-blocking call to listen for a message from a given rank (a rank is the ID for a given process). Then they call their sreq = MPI.Isend function, which is an asynchronous non-blocking call to send a message send_mesg to the chosen rank. When the expected message is found, MPI.Irecv! will then run on its green thread and finish, updating the recv_mesg with the information from the message. However, in order to make sure all of the messages are received, we have added in a blocking operation MPI.Waitall!([rreq, sreq]), which will block all further execution on the given rank until both its rreq and sreq tasks are completed. After that is done, each given rank will have its updated data, and the script will continue on all ranks.

This model is thus very asynchronous and allows for many different computers to run one highly parallelized program, managing the data transmissions in a sparse way without a single computer in charge of managing the whole computation. However, it can be prone to deadlock, since errors in the program may for example require rank 1 to receive a message from rank 2 before continuing the program, but rank 2 won't continue to program until it receives a message from rank 1. For this reason, while MPI has been the most successful large-scale distributed computing model and almost all major high-performance computing (HPC) cluster competitions have been won by codes utilizing the MPI model, the MPI model is nowadays considered a last resort due to these safety issues.

Summary of Multiprocessing

  • Cost is hardware dependent: only suitable for 1ms or higher depending on the connections through which the messages are being passed and the topology of the network.

  • The Master-worker programming model is Julia's Distributed model

  • The Map-reduce programming model is a common data-handling model

  • Array-based distributed computations are another abstraction, used in all forms of parallelism.

  • MPI is a SPMD model of distributed computing, where each process is completely independent and one just controls the memory handling.

The Bait-and-switch: Parallelism is about Programming Models

While this looked like a lecture about parallel programming at the different levels and types of hardware, this wide overview showcases that the real underlying commonality within parallel program is in the parallel programming models, of which there are not too many. There are:

  • Map-reduce parallelism models. pmap, MapReduce (Hadoop/Spark)

    • Pros: Easy to use

    • Cons: Requires that your program is specifically only mapping functions f and reducing them. That said, many data science operations like mean, variance, maximum, etc. can be represented as map-reduce calls, which lead to the popularity of these approaches for "big data" operations.

  • Array-based parallelism models. SIMD (at the compiler level), CuArray, DistributedArray, PyTorch.torch, ...

    • Pros: Easy to use, can have very fast library implementations for specific functions

    • Cons: Less control and restricted to specific functions implemented by the library. Parallelism matches the data structure, so it requires the user to be careful and know the best way to split the data.

  • Loop-based parallelism models. Threads.@threads, @distributed, OpenMP, MATLAB's parfor, Chapel's iterator parallelism

    • Pros: Easy to use, almost no code change can make existing loops parallelized

    • Cons: Refined operations, like locking and sharing data, can be awkward to write. Less control over fine details like scheduling, meaning less opportunities to optimize.

  • Task-based parallelism models with implicit distributed data handling. Threads.@spawn, Dagger.jl, TensorFlow, dask

    • Pros: Relatively high level, low risk of errors since parallelism is mostly handled for the user. User simply describes which functions to call in what order.

    • Cons: When used on distributed systems, implicit data handling is hard, meaning it's generally not as efficient if you don't optimize the code yourself or help the optimizer, and these require specific programming constructs for building the computational graph. Note this is only a downside for distributed data parallelism, whereas when applied to shared memory systems these aspects no longer require handling by the task scheduler.

  • Task-based parallelism models with explicit data handling. Distributed.@spawn

    • Pros: Allows for control over what compute hardware will have specific pieces of data and allows for transferring data manually.

    • Cons: Requires transferring data manually. All computations are managed by a single process/computer/node and thus it can have some issues scaling to extreme (1000+ node) computing situations.

  • SPMD kernel parallelism models. CUDA, MPI, KernelAbstractions.jl

    • Pros: Reduces the problem for the user to only specify what happens in small chunks of the problem. Works on accelerator hardware like GPUs, TPUs, and beyond.

    • Cons: Only works for computations that be represented block-wise, and relies on the compiler to generate good code.

In this sense, the different parallel programming "languages" and features are much more similar than they are all different, falling into similar categories.