Jdi na navigaci předmětu

10: Paralelní výpočty a profilování

Tento notebook je výukovým materiálem v předmětu BI-JUL.21 vyučovaném v zimním semestru akademického roku 2021/2022 Tomášem Kalvodou. Tvorba těchto materiálů byla podpořena NVS FIT.

Hlavní stránkou předmětu, kde jsou i další notebooky a zajímavé informace, je jeho Course Pages stránka.

versioninfo()

1. Úvod

V Julia prostředí se můžeme setkat s čtyřmi druhy paralelismu:

asynchroní "tasky" (coroutines),
více vláknové výpočty (multi-threading),
distribuované výpočty (distributed computing),
výpočty na GPU.

V tomto notebooku se zaměříme na druhý bod a třetí bod. Prezentace bude více experimentální. Zájemci mohou prozkoumat následující balíčky a zdroje

2. Paralelní výpočty

2.1: Více vláknové výpočty

Julia nám umožňuje mezi vlákna rozdělovat tasky, vlákna sdílí paměťové prostředky.

Z příkazové řádky můžeme Julia spustit s více vlákny pomocí přepínače -t/--threads:

$ julia -t4
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.6.3 (2021-09-23)
 _/ |\__'_|_|_|\__'_|  |  
|__/                   |

julia> Threads.nthreads()
4

Pokud chceme pracovat s více vlákny i v Jupyter notebooku/labu tak musíme nastavit proměnnou prostředí JULIA_NUM_THREADS na požadovaný počet vláken. Pokud toto provedete (například pomocí export JULIA_NUM_THREADS=4 na Linuxu), pak by vám výstup následující buňky měl dát 4.

Threads.nthreads()

Pokud jste spustili Jupyter notebook/lab standardním způsobem, pravděpodobně je výstup předešlé buňky 1.

`Threads.@threads` makro

Toto makro z modulu Threads slouží k jednoduchému paralelizování for cyklu. Pokud ho aplikujeme na for cyklus, tak Julia rozdělí práci rovnoměrně (co do počtu prvků v iterátoru) mezi jednotlivá vlákna a počká až všechny vykonají svou práci.

Aktuálně (tj. zima 2021) je k dispozici pouze toto "statické" rozdělení práce, v budoucnou přibudou i další možné strategie, což se už stalo (viz dokumentaci zmíněného makra, k dispozici máme přepínače :dynamic a :static). Toto je samozřejmě vhodné pouze pokud každý z úkonů v cyklu trvá přibližně stejně dlouho.

Jako jednoduchý názorný příklad uvažme pole nul a do každého zapišme id vlákna. Nejprve znovu zkontrolujme, že nám běží jádro s více vlákny:

Threads.nthreads()

Julia vlákna interně čísluje od 1. Demonstrační příklad je tedy následující:

zrs = zeros(Int64, 10)

for j in 1:10
    zrs[j] = Threads.threadid()
end

zrs

# pole o deseti prvcích
L   = 10
zrs = zeros(Int64, L)

Threads.@threads :dynamic for j in 1:L
    # do každého prvku zapíšeme id vlákna
    zrs[j] = Threads.threadid()
end

zrs

# pole o deseti prvcích
L   = 10
zrs = zeros(Int64, L)

Threads.@threads :static for j in 1:L
    # do každého prvku zapíšeme id vlákna
    zrs[j] = Threads.threadid()
end

zrs

Co když každý krok trvá jinou dobu?

# pole o deseti prvcích
L   = 20
zrs = zeros(Int64, L)

Threads.@threads :dynamic for j in 1:L
    sleep(rand()/2)
    # do každého prvku zapíšeme id vlákna
    zrs[j] = Threads.threadid()
end

zrs

# pole o deseti prvcích
L   = 20
zrs = zeros(Int64, L)

Threads.@threads :static for j in 1:L
    sleep(rand()/2)
    # do každého prvku zapíšeme id vlákna
    zrs[j] = Threads.threadid()
end

zrs

`Threads.@spawn` makro, `fetch` a `wait`

@spawn makro z modulu Threads dokáže zadat práci nějakému vláknu bez čekání na výsledek. Jeho návratová hodnota je typu Task, pomocí které můžeme později výsledek získat pomocí metody fetch nebo počkat až výpočet doběhne pomocí metody wait.

using BenchmarkTools
Threads.nthreads()

Jako demonstrační příklad si představme následující ospalou funkci:

"""
Chvilku spinkej a pak vrať 42!
"""
function sleepy_function(t=2)
    sleep(t)
    return 42
end

sleepy_function()

Představme si, že chceme zavolat tuto funkci tolikrát, kolik máme k dispozici vláken a poté výsledky sečíst. Nejprve sériově:

function f()
    values = zeros(Int64, Threads.nthreads())

    for j = 1:Threads.nthreads()
        values[j] = sleepy_function()
    end

    return sum(values)
end

Jaký je výsledek?

f()

Kolik času to zabralo?

@benchmark f()

Osm sekund je přesně to co jsme čekali. Pojďme tuto funkci nyní skutečně paralelizovat a to rovnou dvěma způsoby:

"""
@spawn-fetch verze.
"""
function g()
    values = Array{Task}(undef, Threads.nthreads())

    for j = 1:Threads.nthreads()
        values[j] = Threads.@spawn sleepy_function()
    end

    return sum(fetch.(values))
end

"""
@threads verze (:dynamic).
"""
function h()
    values = zeros(Int64, Threads.nthreads())

    Threads.@threads :dynamic for j = 1:Threads.nthreads()
        values[j] = sleepy_function()
    end

    return sum(values)
end

"""
@threads verze (:static).
"""
function ch()
    values = zeros(Int64, Threads.nthreads())

    Threads.@threads :static for j = 1:Threads.nthreads()
        values[j] = sleepy_function()
    end

    return sum(values)
end

g()

h()

ch()

@benchmark g()

@benchmark h()

@benchmark ch()

Opět, dvě sekundy na čtyřech vláknech jsou přesně to co jsme čekali.

Cvičení, resp. poznámka k 1. domácímu úkolu (z roku 2021)

Porovnejme rychlost několika způsobů jak k jednomu řádku matice přičíst násobek jiného řádku matice (úloha z prvního domácího úkolu). Rozebrme rozdíly mezi jednotlivými přístupy.

function add1!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
    m[target, :] += α * m[source, :]
end

function add2!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
    m[target, :] += α .* m[source, :]
end

function add3!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
    for j in axes(m, 2)
        m[target, j] += α * m[source, j]
    end
end

function add4!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
    @inbounds for j in axes(m, 2)
        m[target, j] += α * m[source, j]
    end
end

function add5!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
    @simd for j in axes(m, 2)
        @inbounds m[target, j] += α * m[source, j]
    end
end

function add6!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
    @simd for j in axes(m, 2)
        m[target, j] += α * m[source, j]
    end
end

function add7!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
    Threads.@threads :dynamic for j in axes(m, 2)
        @inbounds m[target, j] += α * m[source, j]
    end
end

function add8!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
    Threads.@threads :static for j in axes(m, 2)
        @inbounds m[target, j] += α * m[source, j]
    end
end

using BenchmarkTools

mat = 20 * (rand(4, 10^8) .- 1);

@benchmark add1!(a, -2.376, 1, 2) setup=(a=copy(mat))

@benchmark add2!(a, -2.376, 1, 2) setup=(a=copy(mat))

@benchmark add3!(a, -2.376, 1, 2) setup=(a=copy(mat))

@benchmark add4!(a, -2.376, 1, 2) setup=(a=copy(mat))

@benchmark add5!(a, -2.376, 1, 2) setup=(a=copy(mat))

@benchmark add6!(a, -2.376, 1, 2) setup=(a=copy(mat))

@benchmark add7!(a, -2.376, 1, 2) setup=(a=copy(mat))

@benchmark add8!(a, -2.376, 1, 2) setup=(a=copy(mat))

Jaký učiníme závěr z těchto experimentů?

Cvičení: Paralelizace výpočtu Julia fraktálu

Pojďme se vrátit k dřívějšímu příkladu na výpočet Julia fraktálu (notebook bi-jul-07.ipynb). Aby se nám kód dobře paralelizoval, tak ho lehce rozdělíme:

using PyPlot, BenchmarkTools

"""
Výpočet počtu iterací `f` v bodě `z` nutných k opuštění koule o poloměru `R` a středu
v počátku komplexní roviny.
"""
function julia_iterations(z::Complex{Float64}, f::Function, R::Float64; imax::Int64=1000)
    for n = 1:imax
        z = f(z)
        if abs(z) > R
            return n
        end
    end
    
    return imax
end

"""
Sériový výpočet Julia fraktálu.
"""
function julia_matrix!(c::Complex{Float64}, R::Float64, reals, imags, data; imax::Int64=1000)
    f(u) = u^2 + c   
    
    for j = axes(reals, 1), k = axes(imags, 1)
        data[j, k] = julia_iterations(reals[j] + 1im * imags[k], f, R, imax = imax)
    end
end

Následující výpočet zatím není paralelizovaný, ověřme jeho základní funkčnost.

res = LinRange(-2, 2, 1920)
ims = LinRange(-1, 1, 1080)
data = zeros(Int64, length(res), length(ims))
julia_matrix!(-0.8 + 0.256im, 10.0, res, ims, data)

fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])

A podívejme se na dobu běhu.

@benchmark julia_matrix!(-0.8 + 0.256im, 10.0, $res, $ims, $data)

Kolik máme k dispozici vláken? Případně nastavte proměnnou prostředí JULIA_NUM_THREADS před spuštěním Jupyter notebooku/labu.

Threads.nthreads()

Paralelizace se v tomto případě zdá poměrně přímočarou. Ovšem každý "job" (iterace $f$ ) nepoběží úplně stejně dlouho. Pozor, Threads.@threads zatím nepodporuje "složené cykly", musíme kód přepsat pomocí dvou vnořených for cyklů (tj. paralelizujeme plnění matice po řádcích):

function julia_matrix_p1!(c::Complex{Float64}, R::Float64, reals, imags, data; imax::Int64=1000)
    f(u) = u^2 + c   
    
    Threads.@threads for j = axes(reals, 1)
        for k = axes(imags, 1)
            data[j, k] = julia_iterations(reals[j] + 1im * imags[k], f, R, imax = imax)
        end
    end
end

res = LinRange(-2, 2, 1920)
ims = LinRange(-1, 1, 1080)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p1!(-0.8 + 0.256im, 10.0, res, ims, data)

fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])

@benchmark julia_matrix_p1!(-0.8 + 0.256im, 10.0, $res, $ims, $data)

Získali jsme přibližně trojnásobné zrychlení.

function julia_matrix_p2!(c::Complex{Float64}, R::Float64, reals, imags, data; imax::Int64=1000)
    f(u) = u^2 + c   
    
    Threads.@threads :static for j = axes(reals, 1)
        for k = axes(imags, 1)
            data[j, k] = julia_iterations(reals[j] + 1im * imags[k], f, R, imax = imax)
        end
    end
end

res = LinRange(-2, 2, 1920)
ims = LinRange(-1, 1, 1080)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p2!(-0.8 + 0.256im, 10.0, res, ims, data)

fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])

@benchmark julia_matrix_p2!(-0.8 + 0.256im, 10.0, $res, $ims, $data)

Zajímavé, zkusme změnit parametry.

@benchmark julia_matrix_p1!(-0.8 + 0.256im, 5.0, $res, $ims, $data, imax=5_000)

@benchmark julia_matrix_p2!(-0.8 + 0.256im, 5.0, $res, $ims, $data, imax=5_000)

Zkusme aleternativní přístup pomocí @sync a Threads.@spawn. V tomto případě opět ručně iterujeme přes řádky matice.

function julia_matrix_p3!(c::Complex{Float64}, R::Float64, reals, imags, data; imax::Int64=1000)
    f(u) = u^2 + c   
    nthreads = Threads.nthreads()
    
    @sync for t in 1:nthreads
        Threads.@spawn for j = t:nthreads:length(reals), k = axes(imags, 1)
            data[j, k] = julia_iterations(reals[j] + 1im * imags[k], f, R, imax = imax)
        end
    end
end

res = LinRange(-2, 2, 1920)
ims = LinRange(-1, 1, 1080)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p3!(-0.8 + 0.256im, 10.0, res, ims, data)

fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])

@benchmark julia_matrix_p3!(-0.8 + 0.256im, 5.0, $res, $ims, $data, imax=5_000)

To je prakticky stejný výsledek, o fous rychlejší!

res = LinRange(-0.4, 0.4, 2000)
ims = LinRange(-0.6, -0.1, 1600)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p3!(-0.8 + 0.156im, 300.0, res, ims, data, imax=1000)

fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])

res = LinRange(-0.4, 0.4, 2000)
ims = LinRange(-0.6, -0.1, 1600)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p1!(-0.8 + 0.156im, 300.0, res, ims, data, imax=1000)

fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])

2.2 Distribuované výpočty

Pokud chceme využít k výpočtům více strojů, například celé clustery, můžeme k tomu použít modul Distributed ze standardní knihovny. Julia procesy pak běží odděleně (potenciálně i na jiných strojích) a kontrola běhu programu se stává obtížnější/zábavnější.

Procesy jsou opět číslovány od jedné. Ten s id rovným jedné odpovídá procesu v němž běží REPL/Jupyter (interaktivní sezení).

Pokud chceme spustit Julia s více procesy (workery) na jednom stroji, použijeme k tomu přepínač -p/--procs. Například požadujeme-li čtyři procesy na jednom stroji:

$ julia -p 4
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.6.3 (2021-09-23)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> nworkers()
4

Pokud chceme spustit více procesů na více strojích, musíme tyto uvést v souboru a ten předat Julia při spuštění pomocí argumentu --machine-file FILE. Pro komunikaci se používá ssh protokol (bez hesla, klíč), řádek v daném souboru má formát [count*][user@]host[:port] [bind_addr[:port]].user (výchozí hodnota user je aktuální uživatel).

Alternativně lze procesy přidávat/odebírat programově pomocí metod addprocs/rmprocs z modulu Distributed. Viz jejich dokumentaci.

Pokud Julia spustíme tímto způsobem, tak automaticky dojde k importu modulu Distributed.

Logika kontroly výpočtu je podobná jako u více vláknových výpočtů. Máme k dispozici následující makra/metody:

nworkers(): počet procesů (workerů),
workers(): pole s id workerů,
@spawnat ID EXPRESSION: vyhodnotí výraz EXPRESSION na procesu s id ID (může být rovno :any, pokud nám nezáleží na kterém), vrací objekt typu Future.
fetch(future), wait(future): získá výsledek vzdáleného výpočtu future, resp. počká na doběhnutí výpočtu.
@everywhere EXPRESSION: vyhodnotí výraz EXPRESSION na všech procesech (například funkce, s kterými pracujeme, musí být dostupné všem procesům).
pmap(f, array): paralelní map.
@distributed [OP] for: paralelní for cyklus s případnou redukcí operátorem OP.

V tento okamžik doporučuji restartovat jupyter notebooku a spustit ho s jedním vláknem.

Threads.nthreads()

Demonstrační ukázka: náhodná procházka

Aktuálně nám v Jupyter notebooku pravděpodobně běží Julia s jedním procesem/workerem:

using Distributed

nworkers()

Threads.nthreads()

Zvyšme počet procesů na 4:

addprocs(4)
nworkers()

K demonstraci paralelizace zvolme modelování tzv. diskrétní náhodné procházky. Představme si opilého chodce, který vyráží na procházku z bodu $0$ na reálné ose a v každém kroku si s pravděpodobností 0.5 vybere jednotkový krok vpravo, nebo vlevo.

Do jaké vzdálenosti od počátku se dostane po $n=1\,000\,000$ krocích? Jak bude vypadat hustota pravděpodobnosti, že se po těchto $n$ krocích nachází v jistém bodě $k \in \mathbb{Z}$ ?

K experimentálnímu zodpovězení těchto otázek budeme potřebovat tyto procházky simulovat mnohokrát. To je ideální příležitost pro paralelizaci. Sestrojme nejprve sériové řešení:

function random_walk(steps)
    position = 0
    
    for _ in 1:steps
        if rand() >= 0.5
            position += 1
        else
            position -= 1
        end
    end
    
    return position
end

random_walk(0)

random_walk(10)

Vygenerujme sériově 10 takovýchto náhodných procházek o 1000 krocích.

function get_walks(n, steps)
    return [ random_walk(steps) for _ in 1:n ]
end

get_walks(10, 1_000)

Kolik času nám zabere $10\,000$ procházek o stejném počtu kroků?

@time get_walks(10_000, 1_000_000);

Nejprve musíme naší funkci vykonávající náhodnou procházku zadefinovat i v ostatních procesech. Teď by v našem procesu výpočet prošel, ale v jiném ne:

t = @spawnat 1 random_walk(100) # v tento okamžik OK

fetch(t)

t = @spawnat 2 random_walk(100) # v tento okamžik selže

fetch(t)

Stačí použít makro @everywhere (typicky můžeme udělat třeba @everywhere include("common.jl"), kde v daném souboru máme kód pořebný k provádění výpočtů).

@everywhere function random_walk(steps)
    position = 0
    
    for _ in 1:steps
        if rand() >= 0.5
            position += 1
        else
            position -= 1
        end
    end
    
    return position
end

# @everywhere include("source.jl")

Pro jistotu zkontrolujme, že vše funguje tak jak očekáváme.

t = @spawnat 3 random_walk(100)

fetch(t)

K samotné paralelizaci použijeme metodu pmap nebo paralelní for cyklus (makro @distributed). pmap rozděluje práci prvek po prvku a hodí se proto na úlohy, které jsou větší, nebo nerovnoměrné. Naopak @distributed for rozdělí práci mezi workery najednou.

test = zeros(10)
pmap(x -> (sleep(0.2); myid()), test)

workers()

myid()

function get_walks_p1(n, steps)
    output = fill(steps, n)
    output = pmap(random_walk, output)
    return output
end

get_walks_p1(10, 1_000)

Kolik času nám zabere $10\,000$ procházek o stejném počtu kroků?

@time get_walks_p1(10_000, 1_000_000);

Při použití paralelního for cyklu potřebujeme mít možnost zaznamenávat výsledky do sdílené paměti. Přesně k tomu slouží SharedArrays.

using SharedArrays

function get_walks_p2(n, steps)
    output = SharedArray{Int64}(n)
    
    @sync @distributed for j = 1:n
        output[j] = random_walk(steps)
    end
    
    return output
end

get_walks_p2(10, 1_000)

Kolik času nám zabere $10\,000$ procházek o stejném počtu kroků?

@time get_walks_p2(10_000, 1_000_000);

Proč byl bylo řešení s pmap nejpomalejší? Protože úloha, kterou řešíme byla příliš jednoduchá a výhoda paralelizace byla utopena v komunikaci mezi procesy (pmap distribuje každou z úloh zvlášť).

Ukažme si tento efekt na ještě jednodušším příkladě, uvažme jednoduchý výpočet kvadrátu a faktorizaci.

@everywhere begin
    using Primes
    f(n) = n^2       # "levná" operace
    g(n) = factor(n) # občas "drahá" operace
end

A tři různé způsoby výpočtu těchto funkcí v jistém rozsahu.

function serial(func, rng)
    map(func, rng)
end

function para2(func, rng)
    pmap(func, rng)
end

function para3(func, rng)
    @sync @distributed for j in rng
        func(j)
    end
end

integer_range = BigInt(10^18):BigInt(10^18 + 1_000)

Výpočty kvadrátů:

@time serial(f, integer_range);

@time para2(f, integer_range);

@time para3(f, integer_range);

Výpočty faktorizací.

@time serial(g, integer_range);

@time para2(g, integer_range);

@time para3(g, integer_range);

Vraťme se zpět k naší náhodné procházce, napočtětem 10000 procházek o milionu krocích (hodně bloudící opilec):

ws = get_walks_p2(10_000, 1_000_000)

using PyPlot

hist(ws, bins=1000);

Chudák opilec se tedy bude nejpravděpodobněji nacházet v bodě 0. Spočtěme si pár statistických ukazatelů.

using Statistics

mean(ws)

median(ws)

std(ws)

minimum(ws), maximum(ws)

Cvičení: Collatz conjecture

Collatzova hypotéza je pravděpodobně jedním z nejznámějších nevyřešených matematických problémů, které je možné formulovat pomocí naprosto elementárních prostředků. Uvažme zobrazení $f(n) = \begin{cases} n/2, & n \text{ je sudé}, \\ 3n+1, & n \text{ je liché}. \end{cases}$ Collatzova hypotéza pak zní: Pro každé přirozené $n$ existuje přirozené $k$ takové, že $f^{\circ k}(n) = 1$ .

nworkers()

Definujme metodu collatz, která pro dané $n$ provede Collatzovu iteraci a vrátí počet iterací nutný k dosažení $1$

@everywhere function collatz(n)
    i = 0
    
    while n != 1
        if iseven(n)
            n = div(n, 2)
        else
            n = 3n + 1
        end
        
        i += 1
    end
    
    return i
end

Prvních pár hodnot:

[ collatz(n) for n = 1:10 ]

Uvažme opět dvě možné paralelizace (pomocí pmap a pomocí @distributed for):

using SharedArrays

function collatz_pmap(nmax)
    pmap(collatz, 1:nmax)
end

function collatz_dfor(nmax)
    data = SharedArray{Int64}(nmax)
    
    @sync @distributed for n = 1:nmax
        data[n] = collatz(n)
    end
    
    return data
end

@time collatz_pmap(1_000_000)

@time collatz_dfor(1_000_000)

Evidentně jsou jednotlivé iterace stále dost levné, takže druhá varianta vítězí. Napočtěme si pořádnou dávku dat:

@time iterations = collatz_dfor(10^8)

using DataFrames, PyPlot

df = DataFrame(:n => 1:10^8, :iterations => iterations)

Jak vypadá distribuce počtu iterací na prvních pár číslech?

plt.xlabel(L"n")
plt.ylabel("iterations")
scatter(df.n[1:1_000_000], df.iterations[1:1_000_000], s=1)

Jaká je distribuce počtu iterací?

counts = combine(groupby(df, :iterations), nrow => :count)

plt.xlabel("iterations")
plt.ylabel("count")
scatter(counts.iterations, counts.count, s=1)

3. Profilování

Profilování pomáhá při snaze vylepšit výkon naší implementace. Oproti obyčejným benchmarkům, které jsme si ukazovali dříve, se snažíme profilováním zjistit co je podstatné pro dobu běhu programu.

V Julia lze použít modul Profile ze standardní knihovny, který sleduje běh programu a zaznamenává kde běh "tráví nejvíce času". Profile periodicky zaznamenává backtrace aktuálně prováděného řádku kódu a tyto údaje agreguje.

3.1 Ukázka

Názornější bude asi jednoduchá ukázka na nějaké elementární funkci. Uvažme následující funkci, která vyhodnocuje funkční hodnotu polynomu $p(x) = \sum_{k=0}^n a_k x^k$ zadaného pomocí vektoru koeficientů [a[0], a[1], ..., a[n]]. Naivní implementace by byla následující:

function poly_eval_v1(coeffs::Vector{T}, x::T) where { T <: Number }
    value = zero(T)
    for k in 1:length(coeffs)
        value += coeffs[k] * x^(k - 1)
    end
    
    return value
end

poly_eval_v1 (generic function with 1 method)

Před profilování je dobré funkci alespoň jednou vyhodnotit a tím vyvolat její kompilaci (JIT). Jinak budeme profilovat LLVM kompilátor.

poly_eval_v1([1.0, 2.0], 10.0)

21.0

Nejprve musíme importovat zmíněný modul.

using Profile

A poté k měření použít makro @profile. Také použijeme větší vstup, aby bylo co měřit. Zdůrazněme, že Profile nesleduje každý krok programu, jen po jistých časových krocích zkontroluje běh programu. To má jistě své výhody (malé zpomalení samotného výpočtu) i nevýhody (ne zcela přesné měření). Toto makro akumuluje výsledky i z několika běhů. Pokud chceme předchozí data odstranit, použijeme k tomu metodu Profile.clear().

coeffs = rand(10^8)
Profile.clear()
@profile poly_eval_v1(coeffs, -0.8945)

0.46208491495222437

K zobrazení dat o běhu programu můžeme použít hned několik nástrojů. Základní je zobrazení přímo pomocí metody print z modulu Profile:

Profile.print()

Overhead ╎ [+additional indent] Count File:Line; Function
=========================================================
   1╎1    @Base/math.jl:0; pow_body(x::Float64, n::Int64)
    ╎9047 @Base/task.jl:514; (::IJulia.var"#15#18")()
    ╎ 9047 ...lia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket)
    ╎  9047 @Base/essentials.jl:816; invokelatest
    ╎   9047 @Base/essentials.jl:819; #invokelatest#2
    ╎    9047 ...execute_request.jl:67; execute_request(socket::ZMQ.S...
    ╎     9047 ...SoftGlobalScope.jl:65; softscope_include_string(m::...
    ╎    ╎ 9047 @Base/loading.jl:1903; include_string(mapexpr::typ...
   1╎    ╎  9047 @Base/boot.jl:370; eval
  39╎    ╎   39   In[1]:0; poly_eval_v1(coeffs::Vect...
    ╎    ╎   8886 In[1]:4; poly_eval_v1(coeffs::Vect...
  48╎    ╎    48   ...e/essentials.jl:13; getindex
 121╎    ╎    121  @Base/float.jl:410; *
    ╎    ╎    8717 @Base/math.jl:1165; ^
 182╎    ╎     182  @Base/math.jl:0; pow_body(x::Float64, n::...
    ╎    ╎     23   @Base/math.jl:1171; pow_body(x::Float64, n::...
  23╎    ╎    ╎ 23   @Base/promotion.jl:499; ==
 277╎    ╎     588  @Base/math.jl:1179; pow_body(x::Float64, n::...
    ╎    ╎    ╎ 311  @Base/operators.jl:369; >
 311╎    ╎    ╎  311  @Base/int.jl:83; <
  55╎    ╎     168  @Base/math.jl:1180; pow_body(x::Float64, n::...
    ╎    ╎    ╎ 113  @Base/operators.jl:369; >
 113╎    ╎    ╎  113  @Base/int.jl:83; <
    ╎    ╎     1205 @Base/math.jl:1181; pow_body(x::Float64, n::...
 163╎    ╎    ╎ 163  @Base/float.jl:410; *
1041╎    ╎    ╎ 1042 @Base/float.jl:413; muladd
    ╎    ╎     1482 @Base/math.jl:1182; pow_body(x::Float64, n::...
    ╎    ╎    ╎ 1252 @Base/math.jl:47; two_mul
1252╎    ╎    ╎  1252 @Base/float.jl:410; *
    ╎    ╎    ╎ 230  @Base/math.jl:48; two_mul
    ╎    ╎    ╎  230  ...e/floatfuncs.jl:439; fma
 230╎    ╎    ╎   230  .../floatfuncs.jl:434; fma_llvm
    ╎    ╎     1367 @Base/math.jl:1183; pow_body(x::Float64, n::...
1367╎    ╎    ╎ 1367 @Base/float.jl:408; +
    ╎    ╎     1054 @Base/math.jl:1185; pow_body(x::Float64, n::...
    ╎    ╎    ╎ 1054 @Base/operators.jl:578; *
 860╎    ╎    ╎  860  @Base/float.jl:410; *
    ╎    ╎    ╎  194  @Base/promotion.jl:411; *
 194╎    ╎    ╎   194  @Base/float.jl:410; *
    ╎    ╎     991  @Base/math.jl:1186; pow_body(x::Float64, n::...
    ╎    ╎    ╎ 895  @Base/math.jl:47; two_mul
 895╎    ╎    ╎  895  @Base/float.jl:410; *
    ╎    ╎    ╎ 96   @Base/math.jl:48; two_mul
    ╎    ╎    ╎  96   ...e/floatfuncs.jl:439; fma
  96╎    ╎    ╎   96   .../floatfuncs.jl:434; fma_llvm
    ╎    ╎     104  @Base/math.jl:1187; pow_body(x::Float64, n::...
 104╎    ╎    ╎ 104  @Base/float.jl:408; +
    ╎    ╎     1111 @Base/math.jl:1188; pow_body(x::Float64, n::...
    ╎    ╎    ╎ 1111 @Base/int.jl:512; >>>
1111╎    ╎    ╎  1111 @Base/int.jl:504; >>>
    ╎    ╎     14   @Base/math.jl:1190; pow_body(x::Float64, n::...
   1╎    ╎    ╎ 1    @Base/float.jl:410; *
  13╎    ╎    ╎ 13   @Base/float.jl:413; muladd
    ╎    ╎     428  @Base/math.jl:1191; pow_body(x::Float64, n::...
 214╎    ╎    ╎ 214  ...e/essentials.jl:575; ifelse
    ╎    ╎    ╎ 99   @Base/float.jl:622; isfinite
  99╎    ╎    ╎  99   @Base/float.jl:409; -
 115╎    ╎    ╎ 115  @Base/float.jl:413; muladd
    ╎    ╎   104  In[1]:5; poly_eval_v1(coeffs::Vect...
    ╎    ╎    104  @Base/range.jl:891; iterate
 104╎    ╎     104  @Base/promotion.jl:499; ==
  17╎    ╎   17   @Base/math.jl:1168; pow_body(x::Float64, n::I...
   1╎1    ...ia/src/heartbeat.jl:22; heartbeat_thread(sock::Ptr{Noth...
Total snapshots: 9050. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task.

Alternativně:

Profile.print(format=:flat)

 Count  Overhead File                    Line Function
 =====  ======== ====                    ==== ========
    39        39 In[1]                      ? poly_eval_v1(coeffs::Vector{Flo...
  8886         0 In[1]                      4 poly_eval_v1(coeffs::Vector{Flo...
   104         0 In[1]                      5 poly_eval_v1(coeffs::Vector{Flo...
  9047         1 @Base/boot.jl            370 eval
  9047         0 @Base/essentials.jl      819 #invokelatest#2
    48        48 @Base/essentials.jl       13 getindex
   214       214 @Base/essentials.jl      575 ifelse
  9047         0 @Base/essentials.jl      816 invokelatest
  3486      3486 @Base/float.jl           410 *
  1471      1471 @Base/float.jl           408 +
    99        99 @Base/float.jl           409 -
    99         0 @Base/float.jl           622 isfinite
  1170      1169 @Base/float.jl           413 muladd
   326         0 @Base/floatfuncs.jl      439 fma
   326       326 @Base/floatfuncs.jl      434 fma_llvm
   424       424 @Base/int.jl              83 <
  1111      1111 @Base/int.jl             504 >>>
  1111         0 @Base/int.jl             512 >>>
  9047         0 @Base/loading.jl        1903 include_string(mapexpr::typeof(...
  8717         0 @Base/math.jl           1165 ^
   183       183 @Base/math.jl              ? pow_body(x::Float64, n::Int64)
    17        17 @Base/math.jl           1168 pow_body(x::Float64, n::Int64)
    23         0 @Base/math.jl           1171 pow_body(x::Float64, n::Int64)
   588       277 @Base/math.jl           1179 pow_body(x::Float64, n::Int64)
   168        55 @Base/math.jl           1180 pow_body(x::Float64, n::Int64)
  1205         0 @Base/math.jl           1181 pow_body(x::Float64, n::Int64)
  1482         0 @Base/math.jl           1182 pow_body(x::Float64, n::Int64)
  1367         0 @Base/math.jl           1183 pow_body(x::Float64, n::Int64)
  1054         0 @Base/math.jl           1185 pow_body(x::Float64, n::Int64)
   991         0 @Base/math.jl           1186 pow_body(x::Float64, n::Int64)
   104         0 @Base/math.jl           1187 pow_body(x::Float64, n::Int64)
  1111         0 @Base/math.jl           1188 pow_body(x::Float64, n::Int64)
    14         0 @Base/math.jl           1190 pow_body(x::Float64, n::Int64)
   428         0 @Base/math.jl           1191 pow_body(x::Float64, n::Int64)
  2147         0 @Base/math.jl             47 two_mul
   326         0 @Base/math.jl             48 two_mul
  1054         0 @Base/operators.jl       578 *
   424         0 @Base/operators.jl       369 >
   194         0 @Base/promotion.jl       411 *
   127       127 @Base/promotion.jl       499 ==
   104         0 @Base/range.jl           891 iterate
  9047         0 @Base/task.jl            514 (::IJulia.var"#15#18")()
  9047         0 ...ia/src/eventloop.jl     8 eventloop(socket::ZMQ.Socket)
  9047         0 .../execute_request.jl    67 execute_request(socket::ZMQ.Soc...
     1         1 ...ia/src/heartbeat.jl    22 heartbeat_thread(sock::Ptr{Noth...
  9047         0 .../SoftGlobalScope.jl    65 softscope_include_string(m::Mod...
Total snapshots: 9050. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task.

Na cvičení kontaktně podrobně rozebereme tento výpis. Zde textově by to bylo poměrně komplikované.

Úvodní část souvisí s IJulia komunikací s notebookem. Co jistě v obou výrazech vidíme, je vyčnívající řádek

 5237      5234 @Base/math.jl            920 ^

Tj. náš kód tráví hodně času umocňováním.

Vedle textového znázornění můžeme použít i grafické, například pomocí balíčku ProfileSVG.jl (tj. ] add ProfileSVG).

using ProfileSVG

[ Info: Precompiling ProfileSVG [132c30aa-f267-4189-9183-c8a63c7e05e6]

@profview poly_eval_v1(coeffs, -0.8945)

Opět vidíme, že umocňování je nejnáročnější.

3.2 Vylepšení I

Mohli bychom tedy nějak počet umocňování zmenšit? Zcela jistě! Můžeme se ho úplně zbavit!

function poly_eval_v2(coeffs::Vector{T}, x::T) where { T <: Number }
    value = zero(T)
    powx  = one(T)
    
    for k in 1:length(coeffs)
        value += coeffs[k] * powx
        powx  *= x
    end
    
    return value
end

poly_eval_v2 (generic function with 1 method)

poly_eval_v2([1., 2., 3.], 10.)

321.0

Profile.clear()
@profile poly_eval_v2(coeffs, -0.8945)

0.4620849149522243

Profile.print()

Overhead ╎ [+additional indent] Count File:Line; Function
=========================================================
    ╎12967 @Base/task.jl:514; (::IJulia.var"#15#18")()
    ╎ 12967 ...ia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket)
    ╎  12967 @Base/essentials.jl:816; invokelatest
    ╎   12967 @Base/essentials.jl:819; #invokelatest#2
    ╎    12967 ...execute_request.jl:67; execute_request(socket::ZMQ....
    ╎     12967 ...oftGlobalScope.jl:65; softscope_include_string(m::...
    ╎    ╎ 12967 @Base/loading.jl:1903; include_string(mapexpr::ty...
    ╎    ╎  12967 @Base/boot.jl:370; eval
    ╎    ╎   6245  In[11]:6; poly_eval_v2(coeffs::Vector...
  94╎    ╎    94    @Base/float.jl:410; *
6151╎    ╎    6151  @Base/float.jl:408; +
    ╎    ╎   602   In[11]:7; poly_eval_v2(coeffs::Vector...
 602╎    ╎    602   @Base/float.jl:410; *
6119╎    ╎   6120  In[11]:8; poly_eval_v2(coeffs::Vector...
Total snapshots: 12967. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task.

@profview poly_eval_v2(coeffs, -0.8945)

Teď vidíme, že ^ vůbec nevoláme a z algebraických operací používáme jen + a *.

3.3 Vylepšení II

Nemohli bychom snížit i počet algebraických operací nutných k vypočtení naší funkční hodnoty? Ano, mohli! Stačí k vyhodnocování použít tzv. Hornerovu metodu, což je pouze oficiální název pro využití distributivního zákona (jiné uzávorkování). Například $a x^3 + b x^2 + c x + d = \big((ax + b)x + c\big)x + d.$

"""
Vyhodnocování polynomu pomocí Hornerovy metody.
"""
function poly_eval_v3(coeffs::Vector{T}, x::T) where { T <: Number }
    value = coeffs[end]
    for k in (length(coeffs)-1):-1:1
        value = value * x + coeffs[k]
    end
    
    return value
end

poly_eval_v3

poly_eval_v3([1., 2., 3.], 10.)

321.0

poly_eval_v3([42.], -2.)

42.0

Profile.clear()
@profile poly_eval_v3(coeffs, -0.8945)

0.46208491495222453

Profile.print()

Overhead ╎ [+additional indent] Count File:Line; Function
=========================================================
   ╎267 @Base/task.jl:514; (::IJulia.var"#15#18")()
   ╎ 267 @IJulia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket)
   ╎  267 @Base/essentials.jl:816; invokelatest
   ╎   267 @Base/essentials.jl:819; #invokelatest#2
   ╎    267 .../execute_request.jl:67; execute_request(socket::ZMQ.So...
   ╎     267 .../SoftGlobalScope.jl:65; softscope_include_string(m::Mo...
   ╎    ╎ 267 @Base/loading.jl:1903; include_string(mapexpr::type...
   ╎    ╎  267 @Base/boot.jl:370; eval
   ╎    ╎   84  In[16]:7; poly_eval_v3(coeffs::Vector{...
 84╎    ╎    84  @Base/float.jl:408; +
   ╎    ╎   183 In[16]:8; poly_eval_v3(coeffs::Vector{...
   ╎    ╎    183 @Base/range.jl:891; iterate
182╎    ╎     183 @Base/promotion.jl:499; ==
Total snapshots: 267. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task.

@profview poly_eval_v3(coeffs, -0.8945)

Vidíme, že násobení i sčítání proporcionálně jistě zabírá méně času než v předchozím případě.

Pokud se podíváme na profil tak nás asi zaujmou poslední dva řádky:

   ╎    ╎    258 @Base/range.jl:674; iterate
257╎    ╎     258 @Base/promotion.jl:410; ==

Tyto má dle profilu na svědomí řádek 7 v buňce kde definujeme poly_eval_v3, tj.

value = value * x + coeffs[k]

Co se zde děje? Podle popisu půjde spíše o efekt související s iterátorem, resp. for cyklem. Nemohli bychom iterovat efektivněji?

"""
Vyhodnocování polynomu pomocí Hornerovy metody, alternativní iterování.
"""
function poly_eval_v4(reversed_coeffs::Vector{T}, x::T) where { T <: Number }
    value = reversed_coeffs[1]
    for k = 2:length(reversed_coeffs)
        value = value * x + reversed_coeffs[k]
    end
    
    return value
end

poly_eval_v4

poly_eval_v4([3., 2., 1.], 10.)

321.0

Profile.clear()
reversed_coeffs = reverse(coeffs)
@profile poly_eval_v4(reversed_coeffs, -0.8945)

0.46208491495222453

Profile.print()

Overhead ╎ [+additional indent] Count File:Line; Function
=========================================================
   ╎253 @Base/task.jl:514; (::IJulia.var"#15#18")()
   ╎ 253 @IJulia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket)
   ╎  253 @Base/essentials.jl:816; invokelatest
   ╎   253 @Base/essentials.jl:819; #invokelatest#2
   ╎    253 .../execute_request.jl:67; execute_request(socket::ZMQ.So...
   ╎     253 .../SoftGlobalScope.jl:65; softscope_include_string(m::Mo...
   ╎    ╎ 253 @Base/loading.jl:1903; include_string(mapexpr::type...
   ╎    ╎  253 @Base/boot.jl:370; eval
   ╎    ╎   104 In[22]:7; poly_eval_v4(reversed_coeffs...
103╎    ╎    104 @Base/float.jl:408; +
149╎    ╎   149 In[22]:8; poly_eval_v4(reversed_coeffs...
Total snapshots: 253. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task.

@profview poly_eval_v4(coeffs, -0.8945)

Zdá se, že jsme lehce ubrali. Zkusme jednotlivé implementace nyní porovnat.

using BenchmarkTools

coeffs = rand(10^6);

@benchmark poly_eval_v1($coeffs, -0.5)

BenchmarkTools.Trial: 33 samples with 1 evaluation.
 Range (min … max):  151.220 ms … 167.929 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     153.338 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   154.340 ms ±   3.044 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

     █▁▄▁█ ▁▁▁▁    ▁                                             
  ▆▁▆█████▁████▆▁▁▆█▁▆▁▁▁▁▁▆▁▁▁▁▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆ ▁
  151 ms           Histogram: frequency by time          168 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark poly_eval_v2($coeffs, -0.5)

BenchmarkTools.Trial: 3004 samples with 1 evaluation.
 Range (min … max):  1.404 ms …   3.080 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.573 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.650 ms ± 234.732 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ██▃▁                                                        
  ▄████▇██▇▆▅▆▆▅▅▅▄▄▄▄▄▃▃▃▄▃▃▂▃▂▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁ ▃
  1.4 ms          Histogram: frequency by time        2.43 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark poly_eval_v3($coeffs, -0.5)

BenchmarkTools.Trial: 1876 samples with 1 evaluation.
 Range (min … max):  2.409 ms …   3.993 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.573 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.651 ms ± 241.257 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▂▂▂▁▁ ▁                                                     
  ████████▇██▇▅▆▆▅▅▄▄▃▄▄▃▄▄▄▃▃▃▃▃▃▃▃▃▃▃▃▂▃▃▃▃▂▂▃▃▂▂▂▂▂▂▂▂▂▂▁▂ ▃
  2.41 ms         Histogram: frequency by time        3.46 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark poly_eval_v4($coeffs, -0.5)

BenchmarkTools.Trial: 1917 samples with 1 evaluation.
 Range (min … max):  2.385 ms …   3.572 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.529 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.594 ms ± 209.906 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▆█                                                           
  ██▆███▅▆▅▆▆▅▄▆▄▅▅▄▃▃▃▃▃▃▂▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▂▁▂▂▂▁▁▁▂▂▁▁▁▁▁▁ ▂
  2.39 ms         Histogram: frequency by time        3.25 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

"""
Vyhodnocování polynomu pomocí Hornerovy metody, alternativní iterování.
"""
function poly_eval_final(reversed_coeffs::Vector{T}, x::T) where { T <: Number }
    value = reversed_coeffs[1]
    for k = 2:length(reversed_coeffs)
        value = muladd(value, x, reversed_coeffs[k])
    end
    
    return value
end

poly_eval_final

@benchmark poly_eval_final($coeffs, -0.5)

BenchmarkTools.Trial: 3550 samples with 1 evaluation.
 Range (min … max):  1.226 ms …   2.440 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.336 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.396 ms ± 167.426 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▂█▅▁                                                        
  ▅████▆█▇▇▆▅▅▄▄▄▄▄▄▄▄▃▃▄▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▂▁▂▁▂▁▁▁▁▂▁▁▁▁ ▂
  1.23 ms         Histogram: frequency by time        1.94 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

"""
Vyhodnocování polynomu pomocí Hornerovy metody, alternativní iterování.
"""
function poly_eval_v6(reversed_coeffs::Vector{T}, x::T) where { T <: Number }
    value = reversed_coeffs[1]
    for k = 2:length(reversed_coeffs)
        value *= x
        value += reversed_coeffs[k]
    end
    
    return value
end

poly_eval_v6

@benchmark poly_eval_v6($coeffs, -0.5)

BenchmarkTools.Trial: 1849 samples with 1 evaluation.
 Range (min … max):  2.394 ms …   3.925 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.602 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.687 ms ± 271.034 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   █▁▂▃                                                        
  ▆█████▇▆█▇▇▅▆▅▅▅▅▄▅▅▃▄▃▃▃▃▃▃▃▃▂▂▂▂▂▂▃▃▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▂▁▁▂ ▃
  2.39 ms         Histogram: frequency by time        3.56 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

3.4 Poznámka

Nástrojů pro analýzu profilu je více. Zajímavou alternativou je PProf.jl založený na externím programu pprof od Google.

Řešení některých úkolů

"""
Vyhodnocování polynomu bez umocňování.
"""
function poly_eval_v2(coeffs::Vector{T}, x::T) where { T <: Number }
    value = zero(T)
    power = one(T)
    for k in 1:length(coeffs)
        value += coeffs[k] * power
        power *= x
    end
    
    return value
end

"""
Vyhodnocování polynomu pomocí Hornerovy metody.
"""
function poly_eval_v3(coeffs::Vector{T}, x::T) where { T <: Number }
    value = coeffs[end]
    for k in (length(coeffs)-1):-1:1
        value = value * x + coeffs[k]
    end
    
    return value
end

@everywhere function collatz(n)
    k = 0
    while n != 1
        if isodd(n)
            n = 3n+1
        else
            n = div(n, 2)
        end
        
        k += 1
    end
    
    return k
end

Reference

Announcing composable multi-threaded parallelism in Julia

Profilování se věnuje samostatná stránka Julia dokumentace Profiling.

O distribuovaných výpočtech se více dočtene na stránce Multi-processing and Distributed Computing.

10: Paralelní výpočty a profilování

1. Úvod

2. Paralelní výpočty

2.1: Více vláknové výpočty

Threads.@threads makro

Threads.@spawn makro, fetch a wait

Cvičení, resp. poznámka k 1. domácímu úkolu (z roku 2021)

Cvičení: Paralelizace výpočtu Julia fraktálu

2.2 Distribuované výpočty

Demonstrační ukázka: náhodná procházka

Cvičení: Collatz conjecture

3. Profilování

3.1 Ukázka

3.2 Vylepšení I

3.3 Vylepšení II

3.4 Poznámka

Řešení některých úkolů

Reference

`Threads.@threads` makro

`Threads.@spawn` makro, `fetch` a `wait`