Scientific Computing, TU Berlin, WS 2019/2020, Lecture 26
Jürgen Fuhrmann, WIAS Berlin
Multithreading in Julia¶
- Start a function as a Task on a available thread
Threads.@spawn. wait(task): wait for the completion of the taskfetch(task)wait for the completion of the taks and retrieve result
Packages, setup¶
In [1]:
using BenchmarkTools
using LinearAlgebra
using Printf
using PyPlot
injupyter()=isdefined(Main, :IJulia) && Main.IJulia.inited
Out[1]:
- In order to allow multithreading, one has to start Julia with
the environment variable
JULIA_NUM_THREADSset to the desired number - Figure out the number of threads
In [2]:
Threads.nthreads()
Out[2]:
Starting threads using spawn()¶
A task which just returns the number of the thread executing it
In [3]:
function mytask()
return Threads.threadid()
end
Out[3]:
In [4]:
function threads_hello(;ntasks=10)
println("number of tasks: $(ntasks)")
tasks=[Threads.@spawn mytask() for i=1:ntasks]
for i=1:length(tasks)
ithd=fetch(tasks[i])
println("task #$(i) was executed thread #$(ithd)")
end
end
Out[4]:
In [5]:
threads_hello(ntasks=10)
Multithreaded dot product calculation¶
- Interesting example as one has to collect the result into one variable
- Start with creating a subdivision of the loop length into equal parts
In [6]:
function partition(N,ntasks)
loop_begin=zeros(Int64,ntasks)
loop_end=zeros(Int64,ntasks)
for itask=1:ntasks
ltask=Int(floor(N/ntasks))
loop_begin[itask]=(itask-1)*ltask+1
if itask==ntasks # adjust last task length
ltask=N-(ltask*(ntasks-1))
end
loop_end[itask]=loop_begin[itask]+ltask-1
end
return (loop_begin,loop_end)
end
Out[6]:
Check this
In [7]:
partition(100000,3)
Out[7]:
Calculate part of scalar product from n0 to n1
In [8]:
function mydot(a,b,n0,n1)
result=0.0
for i=n0:n1
result+=a[i]*b[i]
end
result
end
Out[8]:
Calculate scalar product in parallell
In [9]:
function threaded_mydot(a,b,N,ntasks)
loop_begin,loop_end=partition(N,ntasks)
tasks=[Threads.@spawn mydot(a,b,loop_begin[i],loop_end[i]) for i=1:ntasks]
return mapreduce(task->fetch(task),+,tasks)
end
Out[9]:
Compare times, check accuracy
In [10]:
function test_threaded_mydot(;N=400000, ntasks=4)
a=rand(N)
b=rand(N)
res_s=@btime mydot($a,$b,1,$N)
res_p=@btime threaded_mydot($a,$b,$N,$ntasks)
res_s ≈ res_p
end
Out[10]:
In [11]:
test_threaded_mydot(N=400000, ntasks=4)
Out[11]:
Multithreaded dot product calculation based on fork-join model¶
The fork-join model promises that we can skip the scheduling "by hand"
In [12]:
function forkjoin_mydot_primitive(a,b)
result=0.0
N=length(a)
Threads.@threads for i=1:N
result+=a[i]*b[i]
end
result
end
Out[12]:
In [13]:
function test_forkjoin_mydot_primitive(;N=400000)
a=rand(N)
b=rand(N)
res_s=@btime mydot($a,$b,1,$N)
res_p=@btime forkjoin_mydot_primitive($a,$b)
res_s ≈ res_p
end
Out[13]:
In [14]:
test_forkjoin_mydot_primitive(N=400000)
Out[14]:
What was wrong ?
- In the parallel loop, there is no control over the access to
result: - A possible scenario:
- Thread 1 wants to perform
result+=x - Thread 2 wants to perform
result+=y
- Thread 1 reads
resultand gets its actual valuer0 - Thread 1 reads
resultand gets its actual valuer0 - Thread 1 performs its operation and writes back
r0+x - Thread 2 performs its operation and writes back
r0+y
- In the result, we have
result=r0+yinstead of the correct valuer0+x+y
- Thread 1 wants to perform
Atomic variables¶
- Atomic variables are protected from unscheduled update:
- Thread 2 would have to wait until Thread 1 is done with its operation
- This is expensive due to the necessary communication infrastructure
In [15]:
function forkjoin_mydot_atomic(a,b)
result = Threads.Atomic{Float64}(0.0)
N=length(a)
Threads.@threads for i=1:N
Threads.atomic_add!(result,a[i]*b[i])
end
return result[]
end
Out[15]:
In [16]:
function test_forkjoin_mydot_atomic(;N=400000)
a=rand(N)
b=rand(N)
res_s=@btime mydot($a,$b,1,$N)
res_p=@btime forkjoin_mydot_atomic($a,$b)
res_s ≈ res_p
end
Out[16]:
In [17]:
test_forkjoin_mydot_atomic(N=400000)
Out[17]:
... at least it is correct
Reduction variables¶
- Our threaded result was correct, because each thread had its own temporary variable
- Transfer this concept to the fork-join model
- Introduce a reduction variable
In [18]:
function forkjoin_mydot_reduction(a,b)
N=length(a)
result=zeros(Threads.nthreads())
Threads.@threads for i=1:N
ithd=Threads.threadid()
result[ithd]+=a[i]*b[i]
end
return sum(result)
end
function test_forkjoin_mydot_reduction(;N=400000)
a=rand(N)
b=rand(N)
res_s=@btime mydot($a,$b,1,$N)
res_p=@btime forkjoin_mydot_reduction($a,$b)
res_s ≈ res_p
end
Out[18]:
In [19]:
test_forkjoin_mydot_reduction(N=400000)
Out[19]:
- Results hint on the experimental character of the implementation
- Handling of the implicit barriers can be a problem
Schönauer vector triad:¶
d[i]=a[i]+b[i]*c[i]- Vary length of the array
- -> Memory performance issues
In [20]:
function vtriad(N,nrepeat)
a = Array{Float64,1}(undef,N)
b = Array{Float64,1}(undef,N)
c = Array{Float64,1}(undef,N)
d = Array{Float64,1}(undef,N)
Threads.@threads for i=1:N
a[i]=i
b[i]=N-i
c[i]=i
d[i]=-i
end
t_scalar=@elapsed begin
@inbounds @fastmath for j=1:nrepeat
for i=1:N
d[i]=a[i]+b[i]*c[i]
end
end
end
t_parallel=@elapsed begin
for j=1:nrepeat
Threads.@threads for i=1:N
@inbounds @fastmath d[i]=a[i]+b[i]*c[i]
end
end
end
GFlops=N*nrepeat*2.0/1.0e9
@printf("% 10d % 10.3f % 10.3f % 10.3f\n", N, GFlops/t_scalar,GFlops/t_parallel,t_scalar/t_parallel)
return [N, GFlops/t_scalar,GFlops/t_parallel,t_scalar/t_parallel]
GC.gc()
end
Out[20]:
In [21]:
function run_triad()
# Approximate number of FLOPs per measurement
flopcount=5.0e8
# Smallest array size
N0=1000
# Data points per decade (of array size)
ppdec=8
# Number of array size increases
nrun=40
# File header
@printf("# nthreads=%d\n",Threads.nthreads())
@printf("# N S_GFlops/s P_GFlops/s speedup\n")
# Loop over measurements
N=N0
result=zeros(4,nrun)
for irun=0:nrun-1
# Have exact powers of 10 in the measurement
if (irun%ppdec==0)
N=N0
N0*=10
end
# Calculate number of repeats so that overall effort stays constant
nrepeat=flopcount/N
result[:,irun+1].=vtriad(Int(ceil(N)),Int(ceil(nrepeat)))
N=N*10^(1.0/ppdec)
end
PyPlot.figure(1)
PyPlot.semilogx(result[1,:], result[2,:])
PyPlot.grid()
PyPlot.ylim(1,20)
PyPlot.xlabel("Array Size")
PyPlot.ylabel("GFlops/s")
end
Out[21]:
In [22]:
run_triad()
Out[22]:
This notebook was generated using Literate.jl.