Merge branch 'master' of github.com:JuliaTeachingCTU/Scientific-Programming-in-Julia

Author: smidl

docs/make.jl

@@ -93,6 +93,9 @@ lecture_09 = [
 ]
 
 lecture_10 = [
+    "Lecture" => "./lecture_10/lecture.md"
+    "Lab" => "./lecture_10/lab.md"
+    "Homework" => "./lecture_10/hw.md"
 ]
 
 lecture_11 = [

docs/src/lecture_09/lecture.md

@@ -491,13 +491,11 @@ Let's now use the IRTools to insert the timing statements into the code for `foo
 using IRTools: xcall, insert!, insertafter!
 
 ir = @code_ir foo(1.0, 1.0)
-for b in IRTools.blocks(ir)
-    for (v, ex) in b
-        if timable(ex.expr)
-            fname = exportname(ex.expr)
-            insert!(b, v, xcall(LoggingProfiler, :record_start, fname))
-            insertafter!(b, v, xcall(LoggingProfiler, :record_end, fname))
-        end
+for (v, ex) in ir
+    if timable(ex.expr)
+        fname = exportname(ex.expr)
+        insert!(ir, v, xcall(LoggingProfiler, :record_start, fname))
+        insertafter!(ir, v, xcall(LoggingProfiler, :record_end, fname))
     end
 end
 
@@ -531,13 +529,11 @@ profile_fun(f::Core.Builtin, args...) = f(args...)
 
 @dynamo function profile_fun(f, args...)
     ir = IRTools.Inner.IR(f, args...)
-    for b in IRTools.blocks(ir)
-        for (v, ex) in b
-            if timable(ex.expr)
-                fname = exportname(ex.expr)
-                insert!(b, v, xcall(Main, :record_start, fname))
-                insertafter!(b, v, xcall(Main, :record_end, fname))
-            end
+    for (v, ex) in ir
+        if timable(ex.expr)
+            fname = exportname(ex.expr)
+            insert!(ir, v, xcall(LoggingProfiler, :record_start, fname))
+            insertafter!(ir, v, xcall(LoggingProfiler, :record_end, fname))
         end
     end
     for (x, st) in ir

docs/src/lecture_10/42-years-processor-trend.png

@@ -0,0 +1,2 @@
+# [Homework 9: Something parallel](@id hw09)
+Will be published somewhere during the weekend.

docs/src/lecture_10/hw.md

@@ -0,0 +1,75 @@
+using Pkg
+Pkg.activate(@__DIR__)
+using Plots
+using BenchmarkTools
+using Distributed
+using SharedArrays
+
+function juliaset_pixel(z₀, c)
+    z = z₀
+    for i in 1:255
+        abs2(z)> 4.0 && return (i - 1)%UInt8
+        z = z*z + c
+    end
+    return UInt8(255)
+end
+
+function juliaset_column!(img, c, n, colj, j)
+    x = -2.0 + (j-1)*4.0/(n-1)
+    for i in 1:n
+        y = -2.0 + (i-1)*4.0/(n-1)
+        @inbounds img[i, colj] = juliaset_pixel(x+im*y, c)
+    end
+    nothing
+end
+
+function juliaset_columns(c, n, columns)
+    img = Array{UInt8,2}(undef, n, length(columns))
+    for (colj, j) in enumerate(columns)
+        juliaset_column!(img, c, n, colj, j)
+    end
+    img
+end
+
+function juliaset_distributed(x, y, partitions = nworkers(), n = 1000)
+    c = x + y*im
+    columns = Iterators.partition(1:n, div(n, partitions))
+    slices = pmap(cols -> juliaset_columns(c, n, cols), columns)
+    reduce(hcat, slices)
+end
+
+# @btime juliaset_distributed(-0.79, 0.15)
+
+# frac = juliaset_distributed(-0.79, 0.15)
+# plot(heatmap(1:size(frac,1), 1:size(frac,2), frac, color=:Spectral))
+
+
+####
+#   Let's work out the shared array approach
+####
+function juliaset_column!(img, c, n, j)
+    x = -2.0 + (j-1)*4.0/(n-1)
+    for i in 1:n
+        y = -2.0 + (i-1)*4.0/(n-1)
+        @inbounds img[i, j] = juliaset_pixel(x+im*y, c)
+    end
+    nothing
+end
+
+function juliaset_range!(img, c, n, columns)
+    for j in columns
+        juliaset_column!(img, c, n, j)
+    end
+    nothing
+end
+
+function juliaset_shared(x, y, partitions = nworkers(), n = 1000)
+    c = x + y*im
+    columns = Iterators.partition(1:n, div(n, partitions))
+    img = SharedArray{UInt8,2}((n, n))
+    slices = pmap(cols -> juliaset_range!(img, c, n, cols), columns)
+    img
+end
+
+# juliaset_shared(-0.79, 0.15)
+# juliaset_shared(-0.79, 0.15, 16)

docs/src/lecture_10/juliaset_p.jl

@@ -0,0 +1,77 @@
+using Plots, BenchmarkTools
+Threads.nthreads()
+function juliaset_pixel(z₀, c)
+    z = z₀
+    for i in 1:255
+        abs2(z)> 4.0 && return (i - 1)%UInt8
+        z = z*z + c
+    end
+    return UInt8(255)
+end
+
+function juliaset_column!(img, c, n, j)
+    x = -2.0 + (j-1)*4.0/(n-1)
+    for i in 1:n
+        y = -2.0 + (i-1)*4.0/(n-1)
+        @inbounds img[i,j] = juliaset_pixel(x+im*y, c)
+    end
+    nothing
+end
+
+function juliaset_single!(img, c, n)
+    for j in 1:n
+        juliaset_column!(img, c, n, j)
+    end
+    nothing
+end
+
+function juliaset(x, y, n=1000, method = juliaset_single!, extra...)
+    c = x + y*im
+    img = Array{UInt8,2}(undef,n,n)
+    method(img, c, n, extra...)
+    return img
+end
+
+# frac = juliaset(-0.79, 0.15)
+# plot(heatmap(1:size(frac,1),1:size(frac,2), frac, color=:Spectral))
+
+
+@btime juliaset(-0.79, 0.15);
+
+function juliaset_static!(img, c, n)
+    Threads.@threads for j in 1:n
+        juliaset_column!(img, c, n, j)
+    end
+    nothing
+end
+
+@btime juliaset(-0.79, 0.15, 1000, juliaset_static!);
+
+
+using Folds
+function juliaset_folds!(img, c, n)
+    Folds.foreach(j -> juliaset_column!(img, c, n, j), 1:n)
+    nothing
+end
+julia> @btime juliaset(-0.79, 0.15, 1000, juliaset_folds!);
+  16.267 ms (25 allocations: 978.20 KiB)
+
+function juliaset_folds!(img, c, n, nt)
+    parts = collect(Iterators.partition(1:n, cld(n, nt)))
+    Folds.foreach(parts) do ii
+        foreach(j ->juliaset_column!(img, c, n, j), ii)
+    end
+    nothing
+end
+julia> @btime juliaset(-0.79, 0.15, 1000, (args...) -> juliaset_folds!(args..., 16));
+  16.716 ms (25 allocations: 978.61 KiB)
+
+
+using FLoops, FoldsThreads
+function juliaset_folds!(img, c, n)
+    @floop ThreadedEx(basesize = 2) for j in 1:n
+        juliaset_column!(img, c, n, j)
+    end
+    nothing
+end
+@btime juliaset(-0.79, 0.15, 1000, juliaset_folds!);