- - Download this notebook and run it locally on your machine [recommended]. Click here. + Download this notebook and run it locally on your machine [recommended]. Click here.
-
You can also run this notebook in the cloud using Binder. Click here
@@ -25,7 +25,7 @@ EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/SCR
```
```@raw html
-
+
\n";
count = 1
)
- content = replace_colors(content)
- open( html_name, "w" ) do html_file
+ open( html_filepath, "w" ) do html_file
write( html_file, content )
end
return nothing
end
-# Replace colors to match Documenter.jl
-function replace_colors(content)
- #content = replace( content, "--jp-layout-color0: #111111;" => "--jp-layout-color0: #1f2424;")
- #content = replace(content, "--md-grey-900: #212121;" => "--md-grey-900: #282f2f;")
- return content
-end
# Loop over notebooks and generate html and markdown
-notebook_files = glob("*.ipynb", "docs/src/notebooks/")
+notebook_files = glob("*.ipynb", "notebooks/")
for filepath in notebook_files
convert_embedded_img_to_base64(filepath)
- create_md_nb_file(filepath)
filename_with_ext = splitpath(filepath)[end]
filename = splitext(filename_with_ext)[1]
- convert_notebook_to_html(filepath, output_name = filename)
- modify_notebook_html("docs/src/notebook-output/$(filename).html")
+ create_md_nb_file(filename)
+ convert_notebook_to_html("docs/src/notebooks/$filename_with_ext", output_name = filename)
+ modify_notebook_html("docs/src/notebook-html/$(filename).html")
end
makedocs(;
@@ -122,12 +114,12 @@ makedocs(;
canonical="https://fverdugo.github.io/XM_40017",
edit_link="main",),
pages=["Home" => "index.md","Getting started"=>"getting_started_with_julia.md", "Notebooks"=>[
- "Julia Basics" => "julia_basics.md",
- "Tasks and channels" => "julia_async.md",
- "Remote calls and remote channels" => "julia_distributed.md",
- "Matrix Multiplication"=>"matrix_matrix.md",
- "Jacobi" => "jacobi_method.md",
- "Solutions" => "sol_matrix_matrix.md"
+ "Julia Basics" => "notebooks/julia_basics.md",
+ "Tasks and channels" => "notebooks/julia_async.md",
+ "Remote calls and remote channels" => "notebooks/julia_distributed.md",
+ "Matrix Multiplication"=>"notebooks/matrix_matrix.md",
+ "Jacobi" => "notebooks/jacobi_method.md",
+ "Solutions" => "notebooks/sol_matrix_matrix.md"
]],
)
diff --git a/docs/src/getting_started_with_julia.md b/docs/src/getting_started_with_julia.md
index f0886bc..48d52aa 100644
--- a/docs/src/getting_started_with_julia.md
+++ b/docs/src/getting_started_with_julia.md
@@ -1,4 +1,3 @@
-
# Getting started
@@ -381,4 +380,3 @@ We have learned the basics of how to work with Julia. If you want to further dig
- [Package manager](https://pkgdocs.julialang.org/v1/getting-started/)
-
diff --git a/docs/src/jacobi_2D.md b/docs/src/jacobi_2D.md
deleted file mode 100644
index 209f638..0000000
--- a/docs/src/jacobi_2D.md
+++ /dev/null
@@ -1,30 +0,0 @@
-```@meta
-EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/jacobi_2D.ipynb"
-```
-
-```@raw html
---``` - -```@raw html - - -``` diff --git a/docs/src/julia_async.md b/docs/src/julia_async.md deleted file mode 100644 index a6754b8..0000000 --- a/docs/src/julia_async.md +++ /dev/null @@ -1,30 +0,0 @@ -```@meta -EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/julia_async.ipynb" -``` - -```@raw html -
Tip - ---``` - -```@raw html - - -``` diff --git a/docs/src/julia_basics.md b/docs/src/julia_basics.md deleted file mode 100644 index 7bc50b5..0000000 --- a/docs/src/julia_basics.md +++ /dev/null @@ -1,30 +0,0 @@ -```@meta -EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/julia_basics.ipynb" -``` - -```@raw html -Tip - ---``` - -```@raw html - - -``` diff --git a/docs/src/julia_distributed.md b/docs/src/julia_distributed.md deleted file mode 100644 index 5334078..0000000 --- a/docs/src/julia_distributed.md +++ /dev/null @@ -1,30 +0,0 @@ -```@meta -EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/julia_distributed.ipynb" -``` - -```@raw html -Tip - ---``` - -```@raw html - - -``` diff --git a/docs/src/julia_distributed_test.md b/docs/src/julia_distributed_test.md deleted file mode 100644 index 35655b0..0000000 --- a/docs/src/julia_distributed_test.md +++ /dev/null @@ -1,30 +0,0 @@ -```@meta -EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/julia_distributed_test.ipynb" -``` - -```@raw html -Tip - ---``` - -```@raw html - - -``` diff --git a/docs/src/julia_intro.md b/docs/src/julia_intro.md deleted file mode 100644 index 9042f66..0000000 --- a/docs/src/julia_intro.md +++ /dev/null @@ -1,30 +0,0 @@ -```@meta -EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/julia_intro.ipynb" -``` - -```@raw html -Tip - ---``` - -```@raw html - - -``` diff --git a/docs/src/julia_jacobi.md b/docs/src/julia_jacobi.md deleted file mode 100644 index 5dbc94d..0000000 --- a/docs/src/julia_jacobi.md +++ /dev/null @@ -1,30 +0,0 @@ -```@meta -EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/julia_jacobi.ipynb" -``` - -```@raw html -Tip - ---``` - -```@raw html - - -``` diff --git a/docs/src/julia_tutorial.md b/docs/src/julia_tutorial.md deleted file mode 100644 index 497778d..0000000 --- a/docs/src/julia_tutorial.md +++ /dev/null @@ -1,30 +0,0 @@ -```@meta -EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/julia_tutorial.ipynb" -``` - -```@raw html -Tip - ---``` - -```@raw html - - -``` diff --git a/docs/src/matrix_matrix.md b/docs/src/matrix_matrix.md deleted file mode 100644 index 6fc6208..0000000 --- a/docs/src/matrix_matrix.md +++ /dev/null @@ -1,30 +0,0 @@ -```@meta -EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/matrix_matrix.ipynb" -``` - -```@raw html -Tip - ---``` - -```@raw html - - -``` diff --git a/docs/src/notebook-hello.md b/docs/src/notebook-hello.md deleted file mode 100644 index cdcaea3..0000000 --- a/docs/src/notebook-hello.md +++ /dev/null @@ -1,30 +0,0 @@ -```@meta -EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/notebook-hello.ipynb" -``` - -```@raw html -Tip - ---``` - -```@raw html - - -``` diff --git a/docs/src/notebook-output/Jacobi_2D.html b/docs/src/notebook-output/Jacobi_2D.html deleted file mode 100644 index 6787d1f..0000000 --- a/docs/src/notebook-output/Jacobi_2D.html +++ /dev/null @@ -1,16628 +0,0 @@ - - - - - -Tip - -jacobi_2D - - - - - - - - - - - - - - - - - - - - - - - -----------Successive Over-relaxation (SOR)¶
In this section, we want to examine another parallel algorithm called Successive Over-relaxation (SOR).
-The SOR algorithm is an iterative method used to solve Laplace equations. The underlying data structure of the SOR algorithm is a two-dimensional grid, whose elements are updated iteratively through some weighted function that considers the old value as well as the values of neighbouring cells.
-
-
This algorithm is applied, for instance, in physics to simulate the climate or temperature of some object.
-- -
-- - - -Climate Simulation -Temperature of a metal plate -- - -
Sequential algorithm¶
The sequential SOR algorithm is as follows (where
-fis some update function andN,Mare the grid sizes):-grid = zeros(N, M) -grid_new = zeros(N, M) -for step in 1:NSTEPS - for i in 2:N #update grid - for j in 2:M - grid_new[i,j] = f(grid[i,j], grid[i-1,j], grid[i+1,j], grid[i,j-1], grid[i,j+1]) - end - end -end -
Diffusion equation on grid¶
Consider the diffusion of a chemical substance on a two-dimensional grid. The concentration of the chemical is given as $c(x,y)$, a function of the coordinates $x$ and $y$. We will consider a square grid with $0 \leq x,y \leq 1$ and the boundary conditions $c(x,y=1) = 1$ and $c(x,y=0) = 0$. That is, the concentration at the top of the grid is always 1 and the concentration at the very bottom is always 0. Furthermore, in the x-direction we will assume periodic boundary conditions, i.e. $c(x=0, y) = c(x=1,y)$.
-
-
We will take the initial condition $c(x,y) = 0$ for $0 \leq x \leq 1, 0 \leq y < 1$.
-The stable state of the diffusion, that is, when the concentraion does not change anymore, can be described by the Laplace equation
-$$ -\nabla^2 c = 0. -$$Numerically, we can approximate the solution with the Jacobi iteration
-$$ -\frac{1}{4}(c^k_{i,j+1} + c^k_{i,j-1} + c^k_{i+1,j} + c^k_{i-1,j}) = c^{k+1}_{i,j}. -$$The superscript $k$ denotes the $k$-th iteration. The algorithm stepwise updates the cells of the grid, until a steady state is reached. To determine the end of the algorithm, we use the stopping condition
-$$ -\max_{i,j} \lvert c^{k+1}_{i,j} - c^{k}_{i,j} \rvert < \epsilon. -$$That is, we stop when all changes to cell values are smaller than some small number, say $\epsilon = 10^{-5}$.
-Furthermore, for this set of initial and boundary conditions, there exists an analytical solution for the stable state, namely -$$ -c(x,y) = y. -$$ -That is, the concentration profile is the identity function of the y-coordinate.
- ----------Serial Program¶
---- -----In [1]:------ -# Update grid until threshold reached -function update_grid!(new_grid, grid, M, N) - @assert grid[1,:] == grid[M-1,:] - @assert grid[M,:] == grid[2,:] - @assert grid[:,1] == fill(0,M) == new_grid[:,1] - @assert grid[:,N] == fill(1,M) == new_grid[:,N] - # Remove ghost cells - g_left = grid[1:M-2, 2:N-1] - g_right = grid[3:M, 2:N-1] - g_up = grid[2:M-1, 3:N] - g_down = grid[2:M-1, 1:N-2] - # Jacobi iteration - new_grid[2:M-1,2:N-1] = 0.25 * (g_up + g_down + g_left + g_right) - # Update ghost cells - new_grid[1,:] = new_grid[M-1,:] - new_grid[M,:] = new_grid[2,:] - return new_grid -end - -function SOR!(new_grid, grid, ϵ, M, N) - grid_old = true - while true - if grid_old - update_grid!(new_grid, grid, M, N) - else - update_grid!(grid, new_grid, M, N) - end - # Flip boolean value - grid_old = !grid_old - diffs = abs.(new_grid[2:M-1, :] - grid[2:M-1, :]) - if maximum(diffs) < ϵ - break - end - end -end -
-- --- - --- -- -- -Out[1]:- - - - --- -SOR! (generic function with 1 method)
---- -----In [2]:------ -using Test, Plots -N = 100 -M = N + 2 -ϵ = 10^-5 - -# initialize grid -grid = zeros(M,N) -grid[:,N] .= 1 -new_grid = zeros(M,N) -new_grid[:,N] .= 1 - -update_grid!(new_grid, grid, M, N) - -# Test if first iteration successful -@test all(new_grid[:,N-1] .≈ fill(0.25, M)) -@test all(new_grid[:,N-2] .≈ fill(0,M)) -@test all(new_grid[:,N] == fill(1,M)) -@test all(new_grid[:,1] == fill(0,M)) -
-- --- - --- -- -- -Out[2]:- - - - --- -Test Passed---- -----In [3]:------ -N = 100 -M = N + 2 -ϵ = 10^-5 - -# initialize grid -grid = zeros(M,N) -grid[:,N] .= 1 -new_grid = zeros(M,N) -new_grid[:,N] .= 1 - -SOR!(new_grid, grid, ϵ, M, N) -plt = heatmap(transpose(new_grid[2:M-1,:])) -display(plt) -
-- --- - --- -- - - - - -- -- - - -- ----- -----In [4]:------ -# analytical solution -function analytical_solution(N) - # Returns the analytical solution as a square grid of size N - grid = zeros(N,N) - for i in 2:N - grid[:,i] .= (i-1)/(N-1) - end - return grid -end - -# Test if solution is identical with analytical solution -sol = analytical_solution(N) -@test maximum(abs.(sol - new_grid[2:M-1,:])) < 0.01 * N -
-- --- - --- -- -- -Out[4]:- - - - --- -Test Passed----------Parallel Program with MPI¶
---- -----In [5]:------ -#] add MPIClusterManagers DelimitedFiles ---- -----In [6]:------ -# to import MPIManager -using MPIClusterManagers - -# need to also import Distributed to use addprocs() -using Distributed - -# specify, number of mpi workers, launch cmd, etc. -manager=MPIWorkerManager(9) - -# start mpi workers and add them as julia workers too. -addprocs(manager) - -@mpi_do manager begin - -function calculate_partition(p, N, nrows, ncols) - # Calculates the row and column indices of this processor - # Get row and column number for processor p - if mod(p,ncols) == 0 - i = div(p,ncols) - else - i = floor(div(p,ncols)) + 1 - end - j = p - (i-1)*ncols - # Rows - if mod(N,nrows) == 0 - prows = div(N, nrows) - row_range =((i-1) * prows + 1) : (i*prows) - else - # nlower processors get the smaller partition - nlower = nrows - (mod(N,nrows)) - n_floor = floor(div(N,nrows)) - if i <= nlower - row_range =((i-1) * n_floor + 1) : (i*n_floor) - else - row_range = ((i-1) * n_floor + (i - nlower)) : (i*n_floor + (i-nlower)) - end - end - # Columns - if mod(N,ncols) == 0 - prows = div(N, ncols) - col_range =((j-1) * prows + 1) : (j*prows) - else - nlower = ncols - (mod(N,ncols)) - n_floor = floor(div(N,ncols)) - if j <= nlower - col_range =((j-1) * n_floor + 1) : (j*n_floor) - else - col_range = ((j-1) * n_floor + (j - nlower)) : (j*n_floor + (j-nlower)) - end - end - # Add 1 to each column index because of ghost cells - col_range = col_range .+ 1 - return row_range, col_range -end - - -function update_grid(grid) - # Returns the updated grid as M-2 x N-2 matrix where M and N are sizes of grid - M = size(grid, 1) - N = size(grid, 2) - # Remove ghost cells - g_left = grid[2:M - 1,1:N-2] - g_right = grid[2:M - 1, 3:N] - g_up = grid[1:M-2, 2:N-1] - g_down = grid[3:M, 2:N-1] - # Jacobi iteration - return 0.25 * (g_up + g_down + g_left + g_right) -end - - using MPI - comm=MPI.COMM_WORLD - id = MPI.Comm_rank(comm) + 1 - - - M = 50 - N = M + 2 - ϵ = 10^-5 # Stopping threshold - nrows = 3 # Number of grid rows - ncols = 3 # Number of grid columns - n_procs = nrows * ncols - @assert n_procs == MPI.Comm_size(comm) - max_diffs = ones(n_procs) # Differences between iterations - max_diff_buf = MPI.UBuffer(max_diffs,1) # Buffer to store maximum differences - - - # initialize grid - if id == 1 - grid_a = zeros(M,N) - grid_a[1,:] .= 1 - grid_b = zeros(M,N) - grid_b[1,:] .= 1 - else - grid_a = nothing - grid_b = nothing - end - - # Broadcast matrix to other processors - grid_a = MPI.bcast(grid_a, 0, comm) - grid_b = MPI.bcast(grid_b, 0, comm) - - # Determine if processor is in top or bottom row of grid - top_pos = id <= ncols - bottom_pos = id > ((nrows-1) * ncols) - local grid_a_old = false # Grid a is the source grid for the first update - - # Get local partition - ind_rows, ind_cols = calculate_partition(id, M, nrows, ncols) - println("Proc $(id) gets rows $(ind_rows) and columns $(ind_cols)") - - # Determine neighbors - n_left = id - 1 - n_right = id + 1 - n_down = id + ncols - n_up = id - ncols - if mod(id, ncols) == 1 - # Left neighbor is last in row - n_left = id + ncols - 1 - end - if mod(id, ncols) == 0 - # Right neighbor is first in row - n_right = id - ncols + 1 - end - - #println("Proc $(id) has neighbors left $(n_left) and right $(n_right) and up $(n_up) and down $(n_down)") - - local finished = false - - #Perform SOR - while !finished - # Flip old and new grid - grid_a_old = !grid_a_old - - # Determine which grid is updated - if grid_a_old - old_grid = grid_a - new_grid = grid_b - else - old_grid = grid_b - new_grid = grid_a - end - - # send left and right columns - left_ind = first(ind_cols) - right_ind = last(ind_cols) - left_col = old_grid[ind_rows, left_ind] - right_col = old_grid[ind_rows, right_ind] - slreq = MPI.Isend(left_col, comm; dest=n_left-1) - srreq = MPI.Isend(right_col, comm; dest=n_right-1) - - # Send bottom row if not bottom - bottom_ind = last(ind_rows) - if !bottom_pos - bottom_col = old_grid[bottom_ind, ind_cols] - sbreq = MPI.Isend(bottom_col, comm; dest=n_down -1) - end - - # Send top row if not at the top - top_ind = first(ind_rows) - if !top_pos - top_row = old_grid[top_ind, ind_cols] - streq = MPI.Isend(top_row, comm; dest = n_up -1) - end - - # Receive left and right column - left_buf = Array{Float64,1}(undef, length(ind_rows)) - right_buf = Array{Float64,1}(undef, length(ind_rows)) - rlreq = MPI.Irecv!(left_buf, comm; source=n_left -1) - rrreq = MPI.Irecv!(right_buf, comm; source=n_right -1) - - # Receive top row if not at the top - if !top_pos - top_buf = Array{Float64,1}(undef, length(ind_cols)) - rtreq = MPI.Irecv!(top_buf, comm; source=n_up -1) - end - - # Receive bottom row if not at the bottom - if !bottom_pos - bottom_buf = Array{Float64,1}(undef, length(ind_cols)) - rbreq = MPI.Irecv!(bottom_buf, comm; source=n_down -1) - end - - # Wait for results - statlr = MPI.Waitall([rlreq, rrreq], MPI.Status) - old_grid[ind_rows, left_ind - 1] = left_buf - old_grid[ind_rows, right_ind + 1] = right_buf - #println("Proc $(id) received left $(old_grid[ind_rows, left_ind - 1]) and right $(old_grid[ind_rows, right_ind + 1])") - - if !top_pos - statt = MPI.Wait(rtreq) - old_grid[top_ind - 1, ind_cols] = top_buf - #println("Proc $(id) received top $(old_grid[top_ind - 1, ind_cols])") - end - - if !bottom_pos - statb = MPI.Wait(rbreq) - old_grid[bottom_ind + 1, ind_cols] = bottom_buf - #println("Proc $(id) received bottom $(old_grid[bottom_ind + 1, ind_cols])") - end - - # Get local subgrid - if !top_pos & !bottom_pos - local_with_ghosts = old_grid[top_ind - 1 : bottom_ind + 1, left_ind - 1 : right_ind + 1] - lb_row = top_ind - ub_row = bottom_ind - elseif top_pos - local_with_ghosts = old_grid[top_ind : bottom_ind + 1, left_ind - 1 : right_ind + 1] - lb_row = top_ind + 1 - ub_row = bottom_ind - elseif bottom_pos - local_with_ghosts = old_grid[top_ind - 1 : bottom_ind, left_ind - 1 : right_ind + 1] - lb_row = top_ind - ub_row = bottom_ind - 1 - end - - # Perform one step of Jacobi iteration - new_grid[lb_row : ub_row, left_ind : right_ind] = update_grid(local_with_ghosts) - - # Calculate max difference - diffs = abs.(new_grid[lb_row : ub_row, left_ind : right_ind] - old_grid[lb_row : ub_row, left_ind : right_ind]) - maxdiff = maximum(diffs) - - # Gather maxdiffs in processor 1 - MPI.Gather!(maxdiff, max_diff_buf, comm; root=0) - - # First processor determines if threshold is exeeded - if id == 1 - if all(max_diffs .< ϵ) - finished = true - println("THRESHOLD SUBCEEDED - TERMINATE SOR") - end - end - - finished = MPI.bcast(finished, 0, comm) - - if finished - # Set ghost cells to zero again so MPI.Reduce gives correct output - new_grid[ind_rows, left_ind - 1] .= 0.0 - new_grid[ind_rows, right_ind + 1] .= 0.0 - if !bottom_pos - new_grid[bottom_ind + 1, ind_cols] .= 0.0 - end - if !top_pos - new_grid[top_ind - 1, ind_cols] .= 0.0 - end - end - end - - using DelimitedFiles - - # Reduce matrix & store result - if !grid_a_old - sor_result = grid_a - else - sor_result = grid_b - end - - MPI.Reduce!(sor_result, +, comm, root = 0) - sor_result[1,:] .= 1.0 - if id == 1 - writedlm("SOR_result.txt", sor_result) - end - - MPI.Finalize() -end -
-- --- - --- -- - - - -- ---From worker 5: Proc 4 gets rows 17:33 and columns 2:17 - From worker 8: Proc 7 gets rows 34:50 and columns 2:17 - From worker 3: Proc 2 gets rows 1:16 and columns 18:34 - From worker 6: Proc 5 gets rows 17:33 and columns 18:34 - From worker 9: Proc 8 gets rows 34:50 and columns 18:34 - From worker 4: Proc 3 gets rows 1:16 and columns 35:51 - From worker 10: Proc 9 gets rows 34:50 and columns 35:51 - From worker 7: Proc 6 gets rows 17:33 and columns 35:51 - From worker 2: Proc 1 gets rows 1:16 and columns 2:17 - From worker 2: THRESHOLD SUBCEEDED - TERMINATE SOR -
---- -----In [7]:------ -rmprocs(manager); -
--- -----In [8]:------ -using Plots, DelimitedFiles -final_grid = readdlm("SOR_result.txt") -M = size(final_grid, 1) -N = size(final_grid, 2) -plt = heatmap(final_grid[:, 2:N-1]) -display(plt) -
-- --- - --- -- - - - - -- -- - - -- ---- -----In [9]:------ -# Test if solution is identical with analytical solution -sol = analytical_solution(M) -# Bring solution into correct form -sol = reverse(transpose(sol), dims = 1) -@test maximum(abs.(sol - final_grid[:,2:N-1])) < 0.01 * M -
-- --- - --- -- -- -Out[9]:- - - - --- -Test Passed--- - - - - - - - - diff --git a/docs/src/notebook-output/julia_Jacobi.html b/docs/src/notebook-output/julia_Jacobi.html deleted file mode 100644 index 58d3ef6..0000000 --- a/docs/src/notebook-output/julia_Jacobi.html +++ /dev/null @@ -1,15936 +0,0 @@ - - - - - --- -----In [ ]:------ --
julia_jacobi - - - - - - - - - - - - - - - - - - - - - - - -- ------------Contents¶
-
-
- How to parallelize the Jacobi method in Julia -
- ------------Analytical solution¶
$u(x) = (x-L)/L$
- ------------- -
-----------A 3D Laplace equation in a complex domain¶
------------- -
-----------The jacobi method¶
------------ -
---- -----In [1]:------ -function jacobi(n,nsteps) - u = zeros(n+2) - u_new = zeros(n+2) - u[1] = -1 - u[end] = 1 - for istep in 1:nsteps - for i in 2:(n+1) - u_new[i] = 0.5*(u[i-1]+u[i+1]) - end - u[2:end-1] = u_new[2:end-1] - end - u -end -
-- --- - --- -- -- -Out[1]:- - - - --- -jacobi (generic function with 1 method)
----- -----In [2]:------ -jacobi(5,1000) -
-- --- - --- -- -- -Out[2]:- - - - --- -7-element Vector{Float64}: - -1.0 - -0.6666666666666666 - -0.3333333333333333 - 0.0 - 0.3333333333333333 - 0.6666666666666666 - 1.0-----------Parallelization strategy¶
-----------Data partition and dependencies¶
- -u_new[i] = 0.5*(u[i-1]+u[i+1])------------ -
-----------"Ghost" entries¶
------------- -
-----------Remote channels¶
------------- -
-----------Implementation¶
-----------Disclaimer¶
The most natural way of implementing a distributed Jacobi method is to use a library that provides distributed arrays:
- -We will use low-level Julia functions as an academic exercise.
- ----------Adding some procs¶
---- -----In [ ]:------ -using Distributed -
---- -----In [ ]:------ -addprocs(3) -
---------Helper functions¶
---- -----In [ ]:------ -@everywhere function initialize_vectors(n) - @assert mod(n,nworkers()) == 0 - local_n = div(n,nworkers()) - u = zeros(local_n+2) - u_new = zeros(local_n+2) - u[1] = -1.0 - u[end] = 1.0 - u,u_new -end -
---- -----In [ ]:------ -@everywhere function compute_error(u,n) - local_n = length(u)-2 - u_exact = copy(u) - p = myid() - 1 - for local_i in 0:(local_n+1) - global_i = local_i + (p-1)*local_n + 1 - Δ = 1/(n+1) - u_exact[local_i+1] = 2*(Δ*global_i-Δ)-1 - end - maximum(abs.(u[2:end-1].-u_exact[2:end-1])) -end -
---------Distributed Jacobi¶
----- -----In [ ]:------ -function jacobi_dist(n,nsteps) - nw = nworkers() - ftrs_prev_snd = Vector{Future}(undef,nw) - ftrs_next_snd = Vector{Future}(undef,nw) - fun = ()->Channel{Float64}(1) - for w in workers() - p = w-1 - ftrs_prev_snd[p] = @spawnat w RemoteChannel(fun,myid()) - ftrs_next_snd[p] = @spawnat w RemoteChannel(fun,myid()) - end - ftrs = Vector{Future}(undef,nw) - @sync for w in workers() - ftrs[w-1] = @spawnat w jacobi_worker( - n,nsteps,ftrs_prev_snd,ftrs_next_snd) - end - reduce(max,fetch.(ftrs)) -end -
---------Exercise: complete the cell below¶
----- -----In [ ]:------ -@everywhere function jacobi_worker( - n,nsteps,ftrs_prev_snd,ftrs_next_snd) - u, u_new = initialize_vectors(n) - w = myid(); p = w-1 - # TODO: Setup channels - #chn_prev_snd - #chn_prev_rcv - #chn_next_snd - #chn_next_rcv - for step in 1:nsteps - # TODO: Communicate - # Local computations - for i in 2:(length(u)-1) - u_new[i] = 0.5*(u[i-1]+u[i+1]) - end - u[2:end-1] = u_new[2:end-1] - end - @show u - err = compute_error(u,n) -end -
---------Test-driven development¶
----- -----In [ ]:------ -err = jacobi_dist(6,1000) -
----------Remember: remote channels¶
------------- -
-----------Possible enhancements¶
-----------Is it possible to overlap computation and communication?¶
------------- -
----------2nd exercise:¶
----- -----In [ ]:------ -@everywhere function jacobi_worker( - n,nsteps,ftrs_prev_snd,ftrs_next_snd) - u, u_new = initialize_vectors(n) - p = myid()-1 - chn_prev_snd = fetch(ftrs_prev_snd[p]) - chn_next_snd = fetch(ftrs_next_snd[p]) - if myid()!=2 - chn_prev_rcv = fetch(ftrs_next_snd[p-1]) - end - if myid() != nprocs() - chn_next_rcv = fetch(ftrs_prev_snd[p+1]) - end - for step in 1:nsteps - # TODO overlap communication and computation - if myid() != 2 - put!(chn_prev_snd,u[2]) - u[1] = take!(chn_prev_rcv) - end - if myid() != nprocs() - put!(chn_next_snd,u[end-1]) - u[end] = take!(chn_next_rcv) - end - for i in 2:(length(u)-1) - u_new[i] = 0.5*(u[i-1]+u[i+1]) - end - u[2:end-1] = u_new[2:end-1] - end - @show u - err = compute_error(u,n) -end -
----------Solution of the exercises¶
----------1st exercise¶
----- -----In [ ]:------ -@everywhere function jacobi_worker( - n,nsteps,ftrs_prev_snd,ftrs_next_snd) - u, u_new = initialize_vectors(n) - p = myid()-1 - chn_prev_snd = fetch(ftrs_prev_snd[p]) - chn_next_snd = fetch(ftrs_next_snd[p]) - if myid()!=2 - chn_prev_rcv = fetch(ftrs_next_snd[p-1]) - end - if myid() != nprocs() - chn_next_rcv = fetch(ftrs_prev_snd[p+1]) - end - for step in 1:nsteps - if myid() != 2 - put!(chn_prev_snd,u[2]) - u[1] = take!(chn_prev_rcv) - end - if myid() != nprocs() - put!(chn_next_snd,u[end-1]) - u[end] = take!(chn_next_rcv) - end - # Local computations - for i in 2:(length(u)-1) - u_new[i] = 0.5*(u[i-1]+u[i+1]) - end - u[2:end-1] = u_new[2:end-1] - end - @show u - err = compute_error(u,n) -end -
---------2nd exercise¶
--- - - - - - - - - diff --git a/docs/src/notebook-output/julia_Tutorial.html b/docs/src/notebook-output/julia_Tutorial.html deleted file mode 100644 index 0c2a483..0000000 --- a/docs/src/notebook-output/julia_Tutorial.html +++ /dev/null @@ -1,15547 +0,0 @@ - - - - - --- -----In [ ]:------ -@everywhere function jacobi_worker( - n,nsteps,ftrs_prev_snd,ftrs_next_snd) - u, u_new = initialize_vectors(n) - p = myid()-1 - chn_prev_snd = fetch(ftrs_prev_snd[p]) - chn_next_snd = fetch(ftrs_next_snd[p]) - if myid()!=2 - chn_prev_rcv = fetch(ftrs_next_snd[p-1]) - end - if myid() != nprocs() - chn_next_rcv = fetch(ftrs_prev_snd[p+1]) - end - for step in 1:nsteps - if myid() != 2 - put!(chn_prev_snd,u[2]) - t1 = @async u[1] = take!(chn_prev_rcv) - end - if myid() != nprocs() - put!(chn_next_snd,u[end-1]) - t2 = @async u[end] = take!(chn_next_rcv) - end - # Local computations (interior) - for i in 3:(length(u)-2) - u_new[i] = 0.5*(u[i-1]+u[i+1]) - end - # Wait for ghost values - if myid() != 2 - wait(t1) - end - if myid() != nprocs() - wait(t2) - end - # Update near boundary values - for i in [2,length(u)-1] - u_new[i] = 0.5*(u[i-1]+u[i+1]) - end - u[2:end-1] = u_new[2:end-1] - end - @show u - err = compute_error(u,n) -end -
julia_tutorial - - - - - - - - - - - - - - - - - - - - - - - -- ------------Contents¶
-
-
- Using the Julia REPL -
- Running serial and parallel code -
- Installing and managing packages -
----------Introduction¶
Julia has its own way of running code and using packages. Many educational sources about Julia assume that you have this basic knowledge, which can be confusing to new users. In this lesson, we will learn these basic skills so that you can start learning more on Julia.
- -----------Preliminaries¶
This is a tutorial. To follow it:
--
-
- Download and install Julia from https://julialang.org/ -
- Download and install VSCode and its Julia extension. How to: https://www.julia-vscode.org/docs/dev/gettingstarted/ -
- Continue reading this document -
----------The Julia REPL¶
There are several ways of opening Julia depending on your operating system and your IDE, but it is usually as simple as launching the Julia app. With VSCode, open a folder (File > Open Folder). Then, press
-Ctrl+Shift+Pto open the command bar, and executeJulia: Start REPL. If this does not work, make sure you have the Julia extension for VSCode installed. Independently of the method you use, opening Julia results in a window with some text ending with:
-julia> -You have just opened the Julia read-evaluate-print loop, or simply the Julia REPL. Congrats! You will spend most of time using the REPL, when working in Julia. The REPL is a console waiting for user input. Just as in other consoles, the string of text right before the input area (
-julia>in the case) is called the command prompt or simply the prompt.The usage is as follows:
--
-
- You write some input -
- press enter -
- you get the output -
For instance, try this
-- -julia> 1 + 1 -
----------Hello world!¶
A "Hello world" example looks like this in Julia
--julia> println("Hello, world!") -
Try to run it in the REPL.
- -----------Help mode¶
Curious about what function
-printlndoes? Enter into help mode to look into the documentation. This is done by typing a question mark (?) into the inut field:-julia> ? -
After typing
-?, the command prompt changes tohelp?>. It means we are in help mode. Now, we can type a function name to see its documentation.- -help?> println -
----------Package and shell modes¶
The REPL comes with two more modes, namely package and shell modes. To enter package mode type
--julia> ] -
Package mode is used to install and manage packages. We are going to discuss the package mode in greater detail later. To return back to normal mode press the backspace key several times.
-To enter shell mode type semicolon (
-;)-julia> ; -
The prompt should have changed to
-shell>indicating that we are in shell mode. Now you can type commands that you would normally do on your system command line. For instance,-shell> ls -
will display the contents of the current folder in Mac or Linux. Using shell mode in Windows is not straightforward, and thus not recommended for beginners.
- -----------Running more complex code¶
Real-world Julia programs are not typed in the REPL in practice. They are written in one or more files and included in the REPL. To try this, create a new file called
-hello.jl, write the code of the "Hello world" example above, and save it. If you are using VSCode, you can create the file using File > New File > Julia File. Once the file is saved with the namehello.jl, execute it as follows-julia> include("hello.jl") -
-Warning: Make sure that the file-"hello.jl"is located in the current working directory of your Julia session. You can query the current directory with functionpwd(). You can change to another directory with functioncd()if needed. Also, make sure that the file extension is.jl. -The recommended way of running Julia code is using the REPL as we did. But it is also possible to run code directly from the system command line. To this end, open a terminal and call Julia followed buy the path to the file containing the code you want to execute.
-
-$ julia hello.jl -Previous line assumes that you have Julia properly installed in the system and that is usable from the terminal. In UNIX systems (Linux and Mac), the Julia binary needs to be in one of the directories listed in the
-PATHenvironment variable. To check that Julia is properly installed, you can use
-$ julia --version -If this runs without error and you see a version number, you are good to go!
--Tip: In this tutorial, when a code snipped starts with-$, it should be run in the terminal. Otherwise, the code is to be run in the Julia REPL. --Tip: Avoid calling Julia code from the terminal, use the Julia REPL instead! Each time you call Julia from the terminal, you start a fresh Julia session and Julia will need to compile your code from scratch. This can be time consuming for large projects. In contrast, if you execute code in the REPL, Julia will compile code incrementally, which is much faster. Running code in a cluster (like in DAS-5 for the Julia assignment) is among the few situations you need to run Julia code from the terminal. -- -----------Running parallel code¶
Since we are in a parallel computing course, let's run a parallel "hello world" example in Julia. Open a Julia REPL and write
--julia> using Distributed -julia> @everywhere println("Hello, world! I am proc $(myid()) from $(nprocs())") -
Here, we are using the
-Distributedpackage, which is part of the Julia standard library that provides distributed memory parallel support. The code prints the process id and the number of processes in the current Julia session.You will provably only see output from 1 proces. We need to add more processes to run the example in parallel. This is done with the
-addprocsfunction.-julia> addprocs(3) -
We have added 3 new processes, plus the old one, we have 4 processes. Run the code again.
--julia> @everywhere println("Hello, world! I am proc $(myid()) from $(nprocs())") -
Now, you should see output from 4 processes.
-It is possible to specify the number of processes when starting Julia from the terminal with the
--pargument (useful, e.g., when running in a cluster). If you launch Julia from the terminal as
-$ julia -p 3 -and then run
--julia> @everywhere println("Hello, world! I am proc $(myid()) from $(nprocs())") -
You should get output from 4 processes as before.
- -----------Installing packages¶
One of the most useful features of Julia is its package manager. It allows one to install Julia packages in a straightforward and platform independent way. To illustrate this, let us consider the following parallel "Hello world" example.
-Copy the following block of code into a new file named
-"hello_mpi.jl"-# file hello_mpi.jl -using MPI -MPI.Init() -comm = MPI.COMM_WORLD -rank = MPI.Comm_rank(comm) -nranks = MPI.Comm_size(comm) -println("Hello world, I am rank $rank of $nranks") -
As you can see from this example, one can access MPI from Julia in a clean way, without type annotations and other complexities of C/C++ code.
-Now, run the file from the REPL
--julia> incude("hello_mpi.jl") -
It provably didn't work, right? Read the error message and note that the MPI package needs to be installed to run this code.
-To install a package, we need to enter package mode. Remember that we entered into help mode by typing
-?. Package mode is activated by typing]-julia> ] -
At this point, the promp should have changed to
-(@v1.8) pkg>indicating that we are in package mode. The text between parenthesis indicates which is the active project, i.e., where packages are going to be installed. In this case, we are working with the global project associated with our Julia installation (which is Julia 1.8 in this example, but it can be another version in your case).To install the MPI package, type
--(@v1.8) pkg> add MPI -
Congrats, you have installed MPI!
-- Tip: Many Julia package names end with-.jl. This is just a way of signaling that a package is written in Julia. When using such packages, the.jlneeds to be ommited. In this case, we have isntalled theMPI.jlpackage even though we have only typedMPIin the REPL. -- Tip: -The package you have installed it is the Julia interface to MPI, called-MPI.jl. Note that it is not a MPI library by itself. It is just a thin wrapper between MPI and Julia. To use this interface, you need an actual MPI library installed in your system such as OpenMPI or MPICH. Julia downloads and installs a MPI library for you, but it is also possible to use a MPI library already available in your system. This is useful, e.g., when running on HPC clusters. See the documentation ofMPI.jlfor further details. -To check that the package was installed properly, exit package mode by pressing the backspace key several times, and run it again
--julia> incude("hello_mpi.jl") -
Now, it should work, but you provably get output from a single MPI rank only.
- -----------Running MPI code¶
To run MPI applications in parallel, you need a launcher like
-mpiexec. MPI codes written in Julia are not an exception to this rule. From the system terminal, you can run
-$ mpiexec -np 4 julia hello_mpi.jl -But it will provably don't work since the version of
-mpiexecneeds to match with the MPI version we are using from Julia. You can find the path to thempiexecbinary you need to use with these commands-julia> using MPI -julia> MPI.mpiexec_path -
and then try again
-
-$ /path/to/my/mpiexec -np 4 julia hello_mpi.jl -with your particular path.
-However, this is not very convenient. Don't worry if you could not make it work! A more elegant way to run MPI code is from the Julia REPL directly, by using these commands:
--julia> using MPI -julia> mpiexec(cmd->run(`$cmd -np 4 julia hello_mpi.jl`)) -
Now, you should see output from 4 ranks.
- -----------Installing packages locally¶
We have installed the
-MPIpackage globally and it will be available in all Julia sessions. However, in some situations, we want to work with different versions of the same package or to install packages in an isolated way to avoid potential conflicts with other packages. This can be done by using local projects.A project is simply a folder in the hard disk. To use a particular folder as your project, you need to activate it. This is done by entering package mode and using the
-activatecommand followed by the path to the folder you want to activate.-(@v1.8) pkg> activate . -
Previous command will activate the current working directory. Note that the dot
-.is indeed the path to the current folder.The prompt has changed to
-(lessons) pkg>indicating that we are in the project within thelessonsfolder. The particular folder name can be different in your case.- Tip: You can activate a project directly when opening Julia from the terminal using the---projectflag. The command$ julia --project=.will open Julia and activate a project in the current directory. You can also achieve the same effect by setting the environment variableJULIA_PROJECTwith the path of the folder you want to activate. -- Tip: The active project folder and the current working directory are two independent concepts! For instance,-(@v1.8) pkg> activate folderBand thenjulia> cd("folderA"), will activate the project infolderBand change the current working directory tofolderA. -At this point all package-related operations will be local to the new project. For instance, install the
-DataFramespackage.-(lessons) pkg> add DataFrames -
Use the package to check that it is installed
--julia> using DataFrames -julia> DataFrame(a=[1,2],b=[3,4]) -
Now, we can return to the global project to check that
-DataFrameshas not been installed there. To return to the global environment, useactivatewithout a folder name.-(lessons) pkg> activate -
The prompt is again
-(@v1.8) pkg>Now, try to use
-DataFrames.-julia> using DataFrames -julia> DataFrame(a=[1,2],b=[3,4]) -
You should get an error or a warning unless you already had
- -DataFramesinstalled globally.----------Project and Manifest files¶
The information about a project is stored in two files
-Project.tomlandManifest.toml.-
-
-Project.tomlcontains the packages explicitly installed (the direct dependencies)
-
-Manifest.tomlcontains direct and indirect dependencies along with the concrete version of each package.
-
In other words,
-Project.tomlcontains the packages relevant for the user, whereasManifest.tomlis the detailed snapshot of all dependencies. TheManifest.tomlcan be used to reproduce the same envinonment in another machine.You can see the path to the current
-Project.tomlfile by using thestatusoperator (orstin its short form) while in package mode-(@v1.8) pkg> status -
The information about the
-Manifest.tomlcan be inspected by passing the-mflag.- -(@v1.8) pkg> status -m -
----------Installing packages from a project file¶
Project files can be used to install lists of packages defined by others. E.g., to install all the dependencies of a Julia application.
-Assume that a colleague has sent to you a
-Project.tomlfile with this content:
-[deps] -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" -Copy the contents of previous code block into a file called
-Project.tomland place it in an empty folder namednewproject. It is important that the file is namedProject.toml. You can create a new folder from the REPL with-julia> mkdir("newproject") -
To install all the packages registered in this file you need to activate the folder containing your
-Project.tomlfile-(@v1.8) pkg> activate newproject -
and then instantiating it
--(newproject) pkg> instantiate -
The instantiate command will download and install all listed packages and their dependencies in just one click.
- -----------Getting help in package mode¶
You can get help about a particular package operator by writing
-helpin front of it-(@v1.8) pkg> help activate -
You can get an overview of all package commands by typing
-helpalone- -(@v1.8) pkg> help -
----------Package operations in Julia code¶
In some situations it is required to use package commands in Julia code, e.g., to automatize installation and deployment of Julia applications. This can be done using the
-Pkgpackage. For instance-julia> using Pkg -julia> Pkg.status() -
is equivalent to call
-statusin package mode.- -(@v1.8) pkg> status -
-- - - - - - - - - diff --git a/docs/src/notebook-output/julia_async.html b/docs/src/notebook-output/julia_async.html deleted file mode 100644 index 746f8a2..0000000 --- a/docs/src/notebook-output/julia_async.html +++ /dev/null @@ -1,15669 +0,0 @@ - - - - - ---------Conclusion¶
We have learned the basics of how to work with Julia. Now, you should be ready to start learning more on the language. If you want to further dig into the topics we have covered here, you can take a look and the following links
--
-
- Julia Manual https://docs.julialang.org/en/v1/manual/getting-started/ -
- Package manager https://pkgdocs.julialang.org/v1/getting-started/ -
If you want to interact with the Julia community on discourse, sign in at https://discourse.julialang.org/
- -julia_async - - - - - - - - - - - - - - - - - - - - - - - - -----------Contents¶
In this notebook, we will learn the basics of asynchronous programming in Julia. In particular, we will learn about:
--
-
- Tasks -
- Channels -
Understanding these concepts is important to learn later distributed computing.
- ----------Tasks¶
Creating a task¶
A task is a piece of computation work that can be run asynchronously (i.e., that can be run in the background). To create a task, we first need to create a function that represents the work to be done in the task. In next cell, we generate a task that generates and sums two matrices.
- ---- -----In [ ]:------ -function work() - println("Starting work") - sleep(7) - a = rand(3,3) - b = rand(3,3) - r = a + b - println("Finishing work") - r -end -
---- -----In [ ]:------ -t = Task(work) -
---------Scheduling a task¶
The task has been created, but the corresponding work has not started. Note that we do not see any output from function
- -workyet. To run the task we need to schedule it.---- -----In [ ]:------ -schedule(t) -
---------Fetching the task result¶
The task has been executed, but we do not see the result. To get the result we need to fetch it.
- ----- -----In [ ]:------ -fetch(t) -
---------Tasks run asynchronously¶
It is important to note that tasks run asynchronously. To illustrate this let's create and schedule a new task.
- ---- -----In [ ]:------ -t = Task(work) -
---- -----In [ ]:------ -schedule(t) -
---------Note that while the task is running we can execute Julia code. To check this, execute the next two cells while the task is running.
- ---- -----In [ ]:------ -sin(4π)*exp(-0.1) -
---- -----In [ ]:------ -1 + 1 -
----------How is this possible? Tasks run in the brackground and this particular task is sleeping for most of the time. Thus, it is possible to use the current Julia process for other operations while the tasks is sleeping.
- ----------Tasks do not run in parallel¶
It is also important to note that tasks do not run in parallel. We were able to run code while previous tasks was running because the task was idling most of the time. If the task does actual work, the current process will be busy running this task and it is likely that we cannot run other code at the same time. Let's illustrate this with an example. The following code computes an approximation of $\pi$ using Leibniz formula. The quality of the approximation increases with the value of
- -n.---- -----In [ ]:------ -function compute_π(n) - s = 1.0 - for i in 1:n - s += (isodd(i) ? -1 : 1) / (i*2+1) - end - 4*s -end -
---------Call this function with a large number. Note that it will take some time.
- ----- -----In [ ]:------ -compute_π(4_000_000_000) -
---------Create a task that performs this computation.
- ----- -----In [ ]:------ -fun = () -> compute_π(4_000_000_000) -t = Task(fun) -
---------Schedule the tasks and then try to execute the 2nd cell bellow. Note that the current process will be busy running the task.
- ---- -----In [ ]:------ -schedule(t) -
---- -----In [ ]:------ -1+1 -
---------Tasks take a function with no arguments¶
This function needs to have zero arguments, but it can capture variables if needed. If we try to create a task with a function that has arguments, it will result in an error when we schedule it.
- ---- -----In [ ]:------ -add(a,b) = a + b -
--- -----In [ ]:------ -t = Task(add) -
---- -----In [ ]:------ -schedule(t) -
---------If we need, we can capture variables in the function to be run by the task as shown in the next cells.
- ---- -----In [ ]:------ -a = rand(3,3) -b = rand(3,3); -
---- -----In [ ]:------ -fun = () -> a + b -t = Task(fun) -schedule(t) -
---------Useful macro:
@async¶So far, we have created tasks using low-level functions, but there are more convenient ways of creating and scheduling tasks. For instance using the
- -@asyncmacro. This macro is used to run a piece of code asynchronously. Under the hood it puts the code in an anonymous function, creates a task, and schedules it. For instance, the next cell is equivalent to previous one.---- -----In [ ]:------ -@async a + b -
---------Another useful macro:
@sync¶This macro is used to wait for all the tasks created with
- -@asyncin a given block of code.---- -----In [ ]:------ -@sync begin - @async sleep(3) - @async sleep(4) -end -
---------Channels¶
Sending data between tasks¶
Julia provides channels as a way to send data between tasks. A channel is like a FIFO queue in which tasks can put and take values from. In next example, we create a channel and a task that puts five values into the channel. Finally, the task closes the channel.
- ---- -----In [ ]:------ -chnl = Channel{Int}() -
---- -----In [ ]:------ -@async begin - for i in 1:5 - put!(chnl,i) - end - close(chnl) -end -
---------By executing next cell several times, we will get the values from the channel. We are indeed communicating values from two different tasks. If we execute the cell more than 5 times, it will raise an error since the channel is closed.
- ----- -----In [ ]:------ -take!(chnl) -
---------Channels are iterable¶
Instead of taking values from a channel until an error occurs, we can also iterate over the channel in a for loop until the channel is closed.
- ---- -----In [ ]:------ -chnl = Channel{Int}() -
--- -----In [ ]:------ -@async begin - for i in 1:5 - put!(chnl,i) - end - close(chnl) -end -
---- -----In [ ]:------ -for i in chnl - @show i -end -
---------Calls to
put!andtake!are blocking¶Note that
- -put!andtake!are blocking operations. Callingput!blocks the tasks until another task callstake!and viceversa. Thus, we need at least 2 tasks for this to work. If we callput!andtake!from the same task, it will result in a dead lock. We have added a print statement to previous example. Run it again and note howput!blocks until we calltake!.--- -----In [ ]:------ -chnl = Channel{Int}() -
--- -----In [ ]:------ -@async begin - for i in 1:5 - put!(chnl,i) - println("I have put $i") - end - close(chnl) -end -
---- -----In [ ]:------ -take!(chnl) -
---------Buffered channels¶
We can be a bit more flexible and use a buffered channel. In this case,
- -put!will block only if the channel is full andtake!will block if the channel is empty. We repeat previous example, but with a buffered channel of size 2. Note that we can callput!until the channel is full. At this point, we need to wait to until we calltake!which removes an item from the channel, making room for a new item.--- -----In [ ]:------ -buffer_size = 2 -chnl = Channel{Int}(buffer_size) -
--- -----In [ ]:------ -@async begin - for i in 1:5 - put!(chnl,i) - println("I have put $i") - end - close(chnl) -end -
---- -----In [ ]:------ -take!(chnl) -
---------Questions¶
----- -----In [ ]:------ -t = @elapsed compute_π(100_000_000) -
----------Question (Q1): How long will the compute time of next cell be? --
- -a) 10*t -b) t -c) 0.1*t -d) near 0*t ----- -----In [ ]:------ -@time for i in 1:10 - compute_π(100_000_000) -end -
----------Question (Q2): How long will the compute time of next cell be? --
- -a) 10*t -b) t -c) 0.1*t -d) near 0*t - ----- -----In [ ]:------ -@time for i in 1:10 - @async compute_π(100_000_000) -end -
----------Question (Q3): How long will the compute time of next cell be? --
- -a) 10*t -b) t -c) 0.1*t -d) near 0*t ----- -----In [ ]:------ -@time @sync for i in 1:10 - @async compute_π(100_000_000) -end -
----------Question (Q4): How long will the compute time of the 2nd cell be? --
- -a) infinity -b) 1 second -c) near 0 seconds -d) 3 seconds ---- -----In [ ]:------ -buffer_size = 4 -chnl = Channel{Int}(buffer_size) -
---- -----In [ ]:------ -@time begin - put!(chnl,3) - i = take!(chnl) - sleep(i) -end -
----------Question (Q5): How long will the compute time of the 2nd cell be? --
- -a) infinity -b) 1 second -c) near 0 seconds -d) 3 seconds ---- -----In [ ]:------ -chnl = Channel{Int}() -
-- - - - - - - - - diff --git a/docs/src/notebook-output/julia_basics.html b/docs/src/notebook-output/julia_basics.html deleted file mode 100644 index 2145075..0000000 --- a/docs/src/notebook-output/julia_basics.html +++ /dev/null @@ -1,16860 +0,0 @@ - - - - - --- -----In [ ]:------ -@time begin - put!(chnl,3) - i = take!(chnl) - sleep(i) -end -
julia_basics - - - - - - - - - - - - - - - - - - - - - - - - -----------Contents¶
In this notebook, we will cover the basic parts of Julia needed later to learn parallel computing. In particular, we will learn about:
--
-
- Variables -
- Functions -
- Arrays -
For a more general introduction to Julia see the nice tutorials made available by JuliaAcademy here or the official Julia educational resources here.
- -----------Using Jupyter notebooks in Julia¶
We are going to use Jupyter notebooks in this and other lectures. You provably have worked with notebooks (in Python). If not, here are the basic concepts you need to know to follow the lessons.
--Tip: Did you know that Jupyter stands for Julia, Python and R? --How to start a Jupyter nootebook in Julia¶
To run a Julia Jupyther notebook, open a Julia REPL and type
--julia> ] -pkg> add IJulia -julia> using IJulia -julia> notebook() -
A new browser window will open. Navigate to the corresponding notebook and open it.
--Warning: Make sure that the notebook is using the same Julia version as the one you used to launch `IJulia`. If it is not the same, go to Kernel > Change Kernel and choose the right version. -- ----------Running a cell¶
To run a cell, click on a cell and press
- -Shift+Enter. You can also use the "Run" button in the toolbar above.---- -----In [ ]:------ -1+3 -4*5 -
---------As you can see from the output of previous cell, the value of the last line is displayed. We can suppress the output with a semicolon. Try it. Execute next cell.
- ----- -----In [ ]:------ -1+3 -4*5; -
---------Cell order is important¶
Running the two cells below in reverse order won't work (try it).
- ---- -----In [ ]:------ -foo() = "Well done!" -
---- -----In [ ]:------ -foo() -
---------REPL modes¶
This is particular to Julia notebooks. You can use package, help, and shell mode just like in the Julia REPL.
- ---- -----In [ ]:------ -] add MPI -
--- -----In [ ]:------ -? print -
---- -----In [ ]:------ -; ls -
------------ -----In [ ]:------ -a = 1 -
----------When assigning a variable, the value on the right hand side is not copied into the variable. It is just an association of a name with a value (much like in Python).
- ----------Re-assign a variable¶
We can re-assign a variable, even with a value of another type. However, avoid changing the variable type for performance reasons.
- ---- -----In [ ]:------ -a = 2 -
--- -----In [ ]:------ -a = 1.0 -
---- -----In [ ]:------ -a = "Hi!" -
---------Unreachable objects¶
When an object is not associated with a variable any more, it cannot be reached by the user. This can happen, e.g., when we re-assign a variable. Another case is when local variables, e.g in a function, go out of scope. The following line allocates a large array and assigns it to variable a
- ----- -----In [ ]:------ -a = zeros(300000000); -
---------If we re-assign the variable to another value, the large array will be inaccessible.
- ----- -----In [ ]:------ -a = nothing -
---------Garbage collector¶
Luckily, Julia has a garbage collector that deallocates unreachable objects. You don't need to bother about manual deallocation! Julia is not constantly looking for unreachable objects. Thus, garbage collection does not happen instantaneously, but it will happen at some point. You can also explicitly call the garbage collector, but it is almost never done in practice.
- ----- -----In [ ]:------ -GC.gc() -
---------Type declarations are optional¶
Julia knows the type of the object associated with a variable.
- ----- -----In [ ]:------ -a = 1 -typeof(a) -
---------We can annotate types if we want, but this will not improve performance (except in very special situations). Thus, annotating types is not done in practice.
- ----- -----In [ ]:------ -c::Int = 1 -typeof(c) -
---------If you annotate a variable with a type, then it cannot refer to objects of other types.
- ----- -----In [ ]:------ -c = "I am a string" -
- -----------We can use Unicode (UTF-8 encoding) characters in variables and function names.
- ----- -----In [ ]:------ -🐱 = "I am a cat" -🐶 = "I am a dog" -🐱 == 🐶 -
---------We can also use Greek letters and other mathematical symbols. Just write the corresponding LaTeX command and press
- -Tab. For example:ωis written\omega+Tab---- -----In [ ]:------ -ω = 1.234 -sin(ω) -
---------In fact, some useful mathematical constants are predefined in Julia with math Greek letters.
- ----- -----In [ ]:------ -sin(π/2) -
----------Question: Which will be the value of `x` in the last line ? (Think your answer before executing next cell to find out the result) -- ----- -----In [ ]:------ -x = 1 -y = x -y = 2 -x -
---------Functions¶
Julia is very much a functional programming language. In consequence, Julia is more centered on functions than on types. This is in contrast to object-oriented languages, which are more centered on types (classes). For instance, you don't need to know the details of the Julia type system to learn parallel programming in Julia, but you need to have a quite advanced knowledge of how Julia functions work.
-Defining functions¶
Functions are defined as shown in next cell. The closing
- -endis necessary. Do to forget it! However, thereturnis optional. The value of last line is returned by default. Indentation is recommended, but it is also optional. That's why the closingendis needed.---- -----In [ ]:------ -function add(a,b) - return a + b -end -
---------Once defined, a function can be called using bracket notation as you would expect.
- ----- -----In [ ]:------ -add(1,3) -
---------Broadcast syntax¶
We can apply functions to arrays element by element using broadcast (dot) syntax.
- ----- -----In [ ]:------ -a = [1,2,3] -b = [4,5,6] -add.(a,b) -
---------Mathematical operators can also be broadcasted (like in Matlab). The following cell won't work. If we want to multiply element by element, we can use the broadcasted version below.
- ----- -----In [ ]:------ -a .* b -
----------Question: Which will be the value of `x` in the last line ? -- ---- -----In [ ]:------ -function q(x) - x = 2 - x -end -
---- -----In [ ]:------ -x = 1 -y = q(x) -x -
---------Defining functions (shorter way)¶
For short functions, we can skip the
- -functionandendkeywords as follows.--- -----In [ ]:------ -add_short(a,b) = a+b -
---- -----In [ ]:------ -add_short(1,3) -
---------Anonymous (lambda) functions¶
Since we can assign function to variables, it is not needed for a function to have a function name in many cases. We can simply create an anonymous function (i.e., a function without name) and assign it to a variable.
- ----- -----In [ ]:------ -add_anonymous = (a,b) -> a+b -
---------We can call the function by using the variable name.
- ----- -----In [ ]:------ -add_anonymous(2.0,3.5) -
---------Note that
- -add_anonymousis not a function name. It is just a variable associated with a function with no function name (well, it has a name technically, but with an arbitrary value).--- -----In [ ]:------ -nameof(add_anonymous) -
---- -----In [ ]:------ -nameof(add) -
---------Functions are first-class objects¶
We can work with Julia functions like with any other type of object. For instance, we can assign functions to variables.
- ----- -----In [ ]:------ -a = add -
---------Now, we can call the function using the variable name.
- ----- -----In [ ]:------ -a(4,5) -
---------We can also create an array of functions (this will not work in Python).
- ----- -----In [ ]:------ -funs = [+,-,*] -
---------To call a specific function in the array, we index the array and then call the returned function
- ----- -----In [ ]:------ -funs[2](2,3) -
---------Higher-order functions¶
Higher order functions are functions that take and/or return other functions. And example is the
-countfunction in Julia.For instance, we can pass a user-defined function to count the number of even elements in an array.
- ----- -----In [ ]:------ -func = i->i%2==0 -a = [1,2,3,5,32,2,4] -count(func,a) -
---------Do-blocks¶
There is yet another way to define anonymous functions. If a function takes a function in its first argument (like
- -count) we can skip the first argument, when calling the function, and define the function we want to pass in a do-block. This is useful, e.g., if we want to define a multi-line anonymous function. The two next cells are equivalent.--- -----In [ ]:------ -function f(i) - m = i%2 - m != 0 -end -count(f,[1,2,3,5,32,2,4]) -
---- -----In [ ]:------ -count([1,2,3,5,32,2,4]) do i - m = i%2 - m != 0 -end -
---------Returning multiple values¶
Julia functions always return a single variable. To return multiple values, we can wrap them in a tuple.
- ----- -----In [ ]:------ -function divrem(a,b) - α = div(a,b) - β = rem(a,b) - (α,β) -end -
---------The output is a tuple as expected, but we can recover the individual values by unpacking the tuple.
- ---- -----In [ ]:------ -d,r = divrem(10,3) -
--- -----In [ ]:------ -d ----- -----In [ ]:------ -r ----------Variable number of input arguments¶
Functions with multiple arguments are also supported. The following example iterates over the given arguments and prints them.
- -argsis just a tuple with all arguments.--- -----In [ ]:------ -function showargs(args...) - for (i,arg) in enumerate(args) - println("args[$i] = $arg") - end -end -
--- -----In [ ]:------ -showargs(1,"Hi!",π) -
---- -----In [ ]:------ -showargs(6) -
---------Positional and keyword arguments¶
Functions can combine positional and keyword arguments much like in Python, but keyword arguments start with semicolon
- -;in Julia.--- -----In [ ]:------ -function foo(a,b;c,d) - println("Positional: a=$a, b=$b. Keyword: c=$c, d=$d") -end -
---- -----In [ ]:------ -foo(3,4,d=2,c=1) -
---------Optional arguments¶
We can provide default values to arguments to make them optional.
- ---- -----In [ ]:------ -function bar(a,b=0;c,d=1) - println("Positional: a=$a, b=$b. Keyword: c=$c, d=$d") -end -
---- -----In [ ]:------ -bar(1,c=2) -
----------Question: Which will be the value of `x` below? -- ---- -----In [ ]:------ -function hofun(x) - y -> x*y -end -
--- -----In [ ]:------ -f2 = hofun(2) -
---- -----In [ ]:------ -x = f2(3) -x -
---------Arrays¶
Julia supports multi-dimensional arrays. They are very similar to Numpy arrays in Python. Let's learn the basics of Julia arrays.
-Array literals¶
We can create (small) arrays from the given values using array literals.
-Next cell creates a vector with 3 integers. Note for Python users: there is no difference between vectors and lists in Julia.
- ----- -----In [ ]:------ -vec = [1,2,3] -
---------We can create a matrix as follows.
- ----- -----In [ ]:------ -mat = [1 2 3 4 - 5 6 7 8 - 9 10 11 12] -
---------Array initialization¶
We can create arrays with all the entries equal to zero, to one, or to a specific given value. The value can be any Julia object, even a function!
- ---- -----In [ ]:------ -zeros(4) -
--- -----In [ ]:------ -zeros(Int,4) -
--- -----In [ ]:------ -ones(2,3) -
--- -----In [ ]:------ -fill(5.0,3,4) -
---- -----In [ ]:------ -fill(add,3,4) -
---------Array comprehensions¶
We can also create the items in the array using a loop within an array comprehension.
- ----- -----In [ ]:------ -squares = [ i^2 for i in 1:8 ] -
---------Indexing¶
We can get and set the items of an array by indexing the array.
- ---- -----In [ ]:------ -squares[3] -
--- -----In [ ]:------ -squares[end] -
--- -----In [ ]:------ -squares[2:4] -
--- -----In [ ]:------ -squares[4] = 16 -
---- -----In [ ]:------ -squares[2:3] = [4,9] -
----------Immutable element type¶
----------Note that once set, the type of the elements in the array cannot be changed. If we try to set an item with an object of a different type, Julia will try to do a conversion, which can fail depending on the passed value.
- ---- -----In [ ]:------ -a = [10,11,12,13] -
---- -----In [ ]:------ -a[2] = "Hi!" -
---------Arrays of any element type¶
Arrays of fixed element type seem to be very rigid, right? Python list have not this limitation. However, we can use arrays of
- -Anytype, which are as flexible as Python lists, or even more since they can also contain functions.--- -----In [ ]:------ -a = Any[10,11,12,13] -
--- -----In [ ]:------ -a[3] = "HI!" -
---- -----In [ ]:------ -a ----------Loops¶
The loop in next cell visits the elements in
- -aone after the other.---- -----In [ ]:------ -a = [10,20,30,40] -for ai in a - @show ai -end -
---------This loop visits the integers from 1 to the length of the array and indexes the array at each of these integers.
- ----- -----In [ ]:------ -for i in 1:length(a) - ai = a[i] - @show (i,ai) -end -
---------This loop "enumerates" the items in the array.
- ----- -----In [ ]:------ -for (i,ai) in enumerate(a) - @show (i,ai) -end -
---------Arrays indices are 1-based by default¶
Be aware of this if you are a C or Python user.
- ----- -----In [ ]:------ -a = [10,20,30,40] -a[0] -
---------Slicing allocates a new array¶
This is also different from Numpy in Python.
- ---- -----In [ ]:------ -a = [1 2 3 - 4 5 6] -
--- -----In [ ]:------ -s = a[:,2] -
--- -----In [ ]:------ -s[2] = 0 -
---- -----In [ ]:------ -a ----------Array views¶
If you want to modify the original array, use
- -viewinstead.--- -----In [ ]:------ -v = view(a,:,2) -
--- -----In [ ]:------ -v[1] = 0 -
---- -----In [ ]:------ -a -----------Exercises¶
----------Exercise 1¶
Implement a function
- -ex1(a)that finds the largest item in the arraya. It should return the largest item and its corresponding position in the array. If there are multiple maximal elements, then the first one will be returned. Assume that the array is not empty. Implement the function in the next cell. Test your implementation with the other one.--- -----In [ ]:------ -# Implement here ----- -----In [ ]:------ -using Test -arr = [3,4,7,3,1,7,2] -@test ex1(arr) == (7,3) -
---------Exercise 2¶
Implement a function
- -ex2(f,g)that takes two functionsf(x)andg(x)and returns a new functionh(x)representing the sum offandg, i.e.,h(x)=f(x)+g(x).--- -----In [ ]:------ -# Implement here ----- -----In [ ]:------ -h = ex2(sin,cos) -xs = LinRange(0,2π,100) -@test all(x-> h(x) == sin(x)+cos(x), xs) -
---------Exercise 3 (hard)¶
Function
- -mandelestimates if a given point(x,y)in the complex plane belongs to the Mandelbrot set.---- -----In [ ]:------ -function mandel(x,y,max_iters) - z = Complex(x,y) - c = z - threshold=2 - for n in 1:max_iters - if abs(z)>threshold - return n-1 - end - z = z^2 +c - end - max_iters -end -
---------If the value of
-mandelis less thanmax_iters, the point is provably outside the Mandelbrot set. Ifmandelis equal tomax_iters, then the point is provably inside the set. The largermax_iters, the better the quality of the estimate (the nicer will be your plot).Plot the value of function
-$$(-1.7,0.7)\times(-1.2,1.2).$$mandelfor each pixel in a 2D grid of the box.Use a grid resolution of at least 1000 points in each direction and
- -max_itersat least 10. You can increase these values to get nicer plots. To plot the values use functionheatmapfrom the Julia packageGLMakie. UseLinRangeto divide the horizontal and vertical axes into pixels. See the documentation of these functions for help.GLMakieis a GPU-accelerated plotting back-end for Julia. It is a large package and it can take some time to install and to generate the first plot. Be patient.-- - - - - - - - - diff --git a/docs/src/notebook-output/julia_distributed.html b/docs/src/notebook-output/julia_distributed.html deleted file mode 100644 index 9509928..0000000 --- a/docs/src/notebook-output/julia_distributed.html +++ /dev/null @@ -1,16411 +0,0 @@ - - - - - --- -----In [ ]:------ -# Implement here -julia_distributed - - - - - - - - - - - - - - - - - - - - - - - -----------
-Programming large-scale parallel systems¶
-----------Distributed computing in Julia¶
-----------Contents¶
In this notebook, we will learn the basics of distributed computing in Julia. In particular, we will focus on the Distributed module available in the Julia standard library. The main topics we are going to cover are:
--
-
- How to create Julia processes -
- How to execute code remotely -
- How to send and receive data -
With this knowledge you will be able to implement simple and complex parallel algorithms in Julia.
- -----------Note: Do not forget to execute the next cell before starting this notebook! -- ----- -----In [ ]:------ -using Printf -function answer_checker(answer,solution) - if answer == solution - "🥳 Well done! " - else - "It's not correct. Keep trying! 💪" - end |> println -end -q_1_check(answer) = answer_checker(answer,"a") -q_2_check(answer) = answer_checker(answer,"b") -
----------How to create Julia processes¶
First of all, we need several processes in order to run parallel algorithms in parallel. In this section, we discuss different ways to create new processes in Julia.
- ----------Adding processes locally¶
The simplest way of creating processes for parallel computing is to add them locally in the current Julia session. This is done by using the following commands.
- ---- -----In [ ]:------ -using Distributed -
---- -----In [ ]:------ -addprocs(3) -
----------Last cell created 3 new Julia processes. By default, they run locally in the same computer as the current Julia session, using multiple cores if possible. However, it is also possible to start the new processes in other machines as long as they are interconnected (more details on this later).
--Tip: We can also start new processes when launching Julia from the command line by suing the `-p` command-line argument. E.g., `$ julia -p 3 ` would launch Julia with 3 extra processes. -- -----------Each process runs a separated Julia instance¶
When adding the new processes, you can imagine that 3 new Julia REPLs have started under the hood (see figure below). The main point of the Distributed module is to provide a way of coordinating all these Julia processes to run code in parallel. It is important to note that each process runs in a separated Julia instance. This means that each process has its own memory space and therefore they do not share memory. This results in distributed-memory parallelism, and allows one to run processes in different machines.
--- -
----------Basic info about processes¶
The following functions provide basic information about the underlying processes. If more than one process is available, the first process is called the main or master and the other the workers. If only a single process is available, it is the master and the first worker simultaneously.
- ---- -----In [ ]:------ -procs() -
--- -----In [ ]:------ -workers() -
--- -----In [ ]:------ -nprocs() -
--- -----In [ ]:------ -nworkers() -
--- -----In [ ]:------ -myid() -
---- -----In [ ]:------ -@everywhere println(myid()) -
----------In previous cell, we have used the macro
- -@everywherethat evaluates the given code on all processes. As a result, each process will print its own process id.----------Creating workers in other machines¶
For large parallel computations, one typically needs to use different computers in parallel. Function
-addprocsalso provides a low-level method to start workers in other machines. Next code example would create 3 workers inserver1and 4 new workers in serverserver2(see figure below). Under the hood, Julia connects via ssh to the other machines and starts the new processes there. In order this to work, the local computer and the remote servers need to be properly configured (see the Julia manual for details).- -using Distributed -machines = [("user@server1",3),("user@server2",4)] -addprocs(machines) -
------------ -
-----------Adding workers with ClusterManagers.jl¶
Previous way of starting workers in other machines is very low-level. Happily, there is a Julia package called ClusterManagers.jl that helps to create workers remotely in number of usual scenarios. For instance, when running the following code from the login node in a computer cluster, it will submit a job to the cluster queue allocating 128 threads. A worker will be generated for each one of these threads. If the compute node have 64 cores, 2 compute nodes will be used to create to contain the 128 workers (see below).
-- -using Distributed -using ClusterManagers -addprocs(SlurmManager(128), partition="debug", t="00:5:00") -
------------ -
-----------Executing code remotely¶
We have added new processes to our Julia session. Let's start using them!
- ----------Function
remotecall¶The most basic thing we can do with a remote processor is to execute a given function on it. This is done by using function
- -remotecall. To make clear how local and remote executions compare, let's call a function locally and then remotely. Next cell uses functiononesto create a matrix locally.---- -----In [ ]:------ -a = ones(2,3) -
---------The next cell does the same operation, but remotely on process 2. Note that
- -remotecalltakes the function we want to execute remotely, the process id where we want to execute it and, finally, the function arguments.---- -----In [ ]:------ -proc = 2 -ftr = remotecall(ones,proc,2,3) -
---------Note that
- -remotecalldoes not return the result of the underlying function, but aFuture. This object represents a reference to a task running on the remote process. To move a copy of the result to the current process we can usefetch.---- -----In [ ]:------ -fetch(ftr) -
---------remotecallis asynchronous¶It is important to note that
- -remotecalldoes not wait for the remote process to finish. It turns immediately. This can be checked be calling remotely the following function that sleeps for 10 secods and then generates a matrix.---- -----In [ ]:------ -fun = (m,n) -> (sleep(10); ones(m,n)) -
---------When running next cell, it will return immediately, event though the remote process will sleep for 10 seconds. We can even run code in parallel. To try this execute the second next cell while the remote call is running in the worker.
- ---- -----In [ ]:------ -proc = 2 -ftr = remotecall(fun,proc,2,3) -
---- -----In [ ]:------ -1+1 -
---------However, when fetching the result, the current process blocks waiting until the result is available in the remote process and arrives to its destination.
- ----- -----In [ ]:------ -fetch(ftr) -
---------Useful macro:
@spawnat¶You have provably realized that in order to use
- -remotecallwe have written auxiliary anonymous functions. They are needed to wrap the code we want to execute remotely. Writing these functions can be tedious. Happily, the macro@spawnatgenerates an auxiliary function from the given block of code and callsremotecallfor us. For instance, the two following cells are equivalent.--- -----In [ ]:------ -@spawnat proc ones(2,3) -
---- -----In [ ]:------ -fun = () -> ones(2,3) -remotecall(fun,proc) -
---------@asyncvs@spawnat¶The relation between
- -@asyncand@pawnatis obvious. From the user perspective they work almost in the same way. However,@asyncgenerates a task that runs asynchronously in the current process, whereas@spawnatexecutes a task in a remote process in parallel. In both cases, the result is obtained usingfetch.--- -----In [ ]:------ -tsk = @async begin - sleep(3) - zeros(2) -end -fetch(tsk) -
---- -----In [ ]:------ -ftr = @spawnat :any begin - sleep(3) - zeros(2) -end -fetch(ftr) -
---------Another usefull macro:
@fetchfrom¶Macro
- -@fetchfromis the blocking version of@spawnat. It blocks and returns the corresponding result instead of aFutureobject.---- -----In [ ]:------ -a = @fetchfrom proc begin - sleep(3) - zeros(2) -end -
----------Data movement¶
Data movement is a crucial part in distributed-memory computations and it is usually one of its main computational -bottlenecks. Being aware of the data we are moving when using functions such as
- -remotecallis important to write efficient distributed algorithms in Julia. Julia also provides a special type of channel, called remote channel, to send and receive data between processes.---------Data movement in
remotecall/fetch¶When usig
- -remotecallwe send to the remote process a function and its arguments. In this example, we send function name+and matricesaandbto proc 4. When fetching the result we receive a copy of the matrix from proc 4.---- -----In [ ]:------ -proc = 4 -a = rand(10,10) -b = rand(10,10) -ftr = remotecall(+,proc,a,b) -fetch(ftr); -
---------Implicit data movement¶
Be aware that data movements can be implicit. This usually happens when we execute remotely functions that capture variables. In the following example, we are also sending matrices
- -aandbto proc 4, even though they do not appear as arguments in the remote call. These variables are captured by the anonymous function and will be sent to proc 4.---- -----In [ ]:------ -proc = 4 -a = rand(10,10) -b = rand(10,10) -fun = () -> a+b -ftr = remotecall(fun,proc) -fetch(ftr); -
---------Data movement with remote channels¶
Another way of moving data between processes is to use remote channels. Their usage is very similar to conventional channels for moving data between tasks, but there are some important differences. In the next cell, we create a remote channel. Process 4 puts several values and closes the channel. Like for conventional channels, calls to
- -put!are blocking, but next cell is not blocking the master process since the call toput!runs asynchronously on process 4.--- -----In [ ]:------ -fun = ()->Channel{Int}() -chnl = RemoteChannel(fun) -
---- -----In [ ]:------ -@spawnat 4 begin - for i in 1:5 - put!(chnl,i) - end - close(chnl) -end; -
---------We can take values from the remote channel form any process using
- -take!. Run next cell several times. The sixth time it should raise and error since the channel was closed.---- -----In [ ]:------ -take!(chnl) -
---------Remote channels can be buffered¶
Just like conventional channels, remote channels can be buffered. The buffer is stored in the process that owns the remote channel. By default this corresponds to process that creates the remote channel, but it can be a different one. For instance, process 3 will be the owner in the following example.
- ---- -----In [ ]:------ -buffer_size = 2 -owner = 3 -fun = ()->Channel{Int}(buffer_size) -chnl = RemoteChannel(fun,owner) -
---- -----In [ ]:------ -@spawnat 4 begin - println("start") - for i in 1:5 - put!(chnl,i) - println("I have put $i") - end - close(chnl) - println("stop") -end; -
---------Note that since the channel is buffered, worker 4 can start putting values into it before any call to
- -take!. Run next cell several times until the channel is closed.---- -----In [ ]:------ -take!(chnl) -
---------Remote channels are not iterable¶
One main difference with respect to conventional channels is that remote channels cannot be iterated. Let's repeat the example above.
- ---- -----In [ ]:------ -fun = ()->Channel{Int}() -chnl = RemoteChannel(fun) -
---- -----In [ ]:------ -@spawnat 4 begin - for i in 1:5 - put!(chnl,i) - end - close(chnl) -end; -
---------Now, try to iterate over the channel in a for loop. It will result in an error since channels are not iterable.
- ----- -----In [ ]:------ -for j in chnl - @show j -end -
---------If we want to take values form a remote channel and stop automatically when the channel is closed, we can combine a while loop and a try-catch statement. This works since
- -take!raises an error if the channel is closed, which will execute thecatchblock and breaks the loop.---- -----In [ ]:------ -while true - try - j = take!(chnl) - @show j - catch - break - end -end -
----------Questions¶
-----------Question (Q1): How many integers are transferred between master and worker? Including both directions. ---a = rand(Int,4,4) -proc = 4 -@fetchfrom proc sum(a^2) -
- -a) 17 -b) 32 -c) 16^2 -d) 65 ----- -----In [ ]:------ -answer = "x" #Replace x with a, b, c, or d -q_1_check(answer) -
----------Question (Q2): How many integers are transferred between master and worker? Including both directions. ---a = rand(Int,4,4) -proc = 4 -@fetchfrom proc sum(a[2,2]^2) -
- -a) 2 -b) 17 -c) 5 -d) 32 ----- -----In [ ]:------ -answer = "x" #Replace x with a, b, c, or d -q_2_check(answer) -
----------Question (Q3): Which value will be the value of `x` ? -- ----- -----In [ ]:------ -a = zeros(Int,3) -proc = 3 -@sync @spawnat proc a[2] = 2 -x = a[2] -x -
----------Question (Q4): Which value will be the value of `x` ? --Which value will be the value of
- -x?---- -----In [ ]:------ -a = zeros(Int,3) -proc = myid() -@sync @spawnat proc a[2] = 2 -x = a[2] -x -
----------Remember: each process runs in a separated Julia instance¶
In particular, this means that each process can load different functions or packages. In consequence, it is important to make sure that the code we run is defined in the corresponding process.
- ----------Functions are defined in a single process¶
This is a very common pitfall when running parallel code. If we define a function in a process, it is not automatically available in the other processes. This is illustrated in the next example. The remote call in the last line in next cell will fail since the function
- -sleep_onesis only being defined in the local process.--- -----In [ ]:------ -function sleep_ones(m,n) - sleep(4) - ones(m,n) -end -
---- -----In [ ]:------ -proc = 3 -remotecall_fetch(sleep_ones,proc,3,4) -
---------To fix this, we can define the function on all processes with the
- -@everywheremacro.--- -----In [ ]:------ -@everywhere function sleep_ones(m,n) - sleep(4) - ones(m,n) -end -
---- -----In [ ]:------ -proc = 3 -remotecall_fetch(sleep_ones,proc,3,4) -
---------Anonymous functions are available everywhere¶
If a function has a name, Julia only sends the function name to the corresponding process. Then, Julia looks for the corresponding function code in the remote process and executes it. This is why the function needs to be defined also in the remote process. However, if a function is anonymous, Julia needs to send the complete function definition to the remote process. This is why anonymous functions do not need to be defined with the macro
- -@everywhereto work in a remote call.--- -----In [ ]:------ -fun = (m,n) -> (sleep(4);ones(m,n)) -
---- -----In [ ]:------ -proc = 3 -remotecall_fetch(fun,proc,3,4) -
---------Each proc uses packages independently¶
When using a package in a process, it is not available in the other ones. For instance, if we load the
- -LinearAlgebrapackage in the current process and use one of its exported functions in another process, we will get an error.--- -----In [ ]:------ -using LinearAlgebra -
---- -----In [ ]:------ -@fetchfrom 3 norm([1,2,3]) -
---------To fix this, we can load the package on all processors with the
- -@everywheremacro.--- -----In [ ]:------ -@everywhere using LinearAlgebra -
---- -----In [ ]:------ -@fetchfrom 3 norm([1,2,3]) -
---------Each process has its own active package environment¶
This is another very common source of errors. You can check that if you activate the current directory, this will have no effect in the other processes.
- ----- -----In [ ]:------ -] activate . -
---------We have activated the current folder. Now let's see which is the active project in another process, say process 2. You will see that process 2 is provably still using the global package environment.
- ---- -----In [ ]:------ -@everywhere using Pkg -
---- -----In [ ]:------ -@spawnat 2 Pkg.status(); -
---------To fix this, you need to activate the current directory on all processes.
- ---- -----In [ ]:------ -@everywhere Pkg.activate(".") -
---- -----In [ ]:------ -@spawnat 2 Pkg.status(); -
----------Easy ways of parallelizing code¶
A part from the low-level parallel routines we have seen so-far, Julia also provides much more simple ways to parallelizing loops and maps.
- ----------Useful macro: @distributed¶
This macro is used when we want to perform a very large for loops made of independent small iterations. To illustrate this, let's consider again the function that computes $\pi$ with Leibniz formula.
- ----- -----In [ ]:------ -function compute_π(n) - s = 1.0 - for i in 1:n - s += (isodd(i) ? -1 : 1) / (i*2+1) - end - 4*s -end -
---------Paralelizing this function might require some work with low-level functions like
- -remotecall, but it is trivial using the macro@distributed. This macro runs the for loop using the available processes and optionally reduces the result using a given reduction function (+in this case).---- -----In [ ]:------ -function compute_π_dist(n) - s = 1.0 - r = @distributed (+) for i in 1:n - (isodd(i) ? -1 : 1) / (i*2+1) - end - 4*(s+r) -end -
---------Run next cell to measure the performance of the serial function for a large value of
- -n. Run it at least 2 times to get rid of compilation times.---- -----In [ ]:------ -@time compute_π(4_000_000_000) -
---------Run next cell to measure the performance of the parallel function.
- ----- -----In [ ]:------ -@time compute_π_dist(4_000_000_000) -
---------Useful function:
pmap¶This function is used when we want to call a very expensive function a small number of evaluations and we want to distribute these evaluations over the available processes. To illustrate the usage of
- -pmapconsider the following example. Next cell generates sixty 30x30 matrices. The goal is to compute the singular value decomposition of all of them. This operation is known to be expensive for large matrices. Thus, this is a perfect scenario forpmap.---- -----In [ ]:------ -a = [ rand(300,300) for i in 1:60]; -
---------First, lets measure the serial performance
- ---- -----In [ ]:------ -using LinearAlgebra -
---- -----In [ ]:------ -@time svd.(a); -
---------If we use
- -pmapinstead of broadcast, the different calls tosvdwill be distributed over the available processes.---- -----In [ ]:------ -@time pmap(svd,a); -
---------Summary¶
We have seen the basics of distributed computing in Julia. The programming model is essentially an extension of tasks and channels to parallel computations on multiple machines. The low-level functions are
- -remotecallandRemoteChannel, but there are other functions and macros likepmapand@distributedthat simplify the implementation of parallel algorithms.-- - - - - - - - - diff --git a/docs/src/notebook-output/julia_distributed_test.html b/docs/src/notebook-output/julia_distributed_test.html deleted file mode 100644 index 1533c1c..0000000 --- a/docs/src/notebook-output/julia_distributed_test.html +++ /dev/null @@ -1,16868 +0,0 @@ - - - - - --- -----In [ ]:------ --
julia_distributed_test - - - - - - - - - - - - - - - - - - - - - - - -- ------------Contents¶
In this notebook, we will learn the basics of distributed computing in Julia. In particular, we will focus on the Distributed module available in the Julia standard library. The main topics we are going to cover are:
--
-
- How to create Julia processes -
- How to execute code remotely -
- How to send and receive data -
With this knowledge you will be able to implement simple and complex parallel algorithms in Julia.
- -----------How to create Julia processes¶
First of all, we need several processes in order to run parallel algorithms in parallel. In this section, we discuss different ways to create new processes in Julia.
- ----------Adding processes locally¶
The simplest way of creating processes for parallel computing is to add them locally in the current Julia session. This is done by using the following commands.
- ---- -----In [ ]:------ -using Distributed -
---- -----In [ ]:------ -addprocs(3) -
----------Last cell created 3 new Julia processes. By default, they run locally in the same computer as the current Julia session, using multiple cores if possible. However, it is also possible to start the new processes in other machines as long as they are interconnected (more details on this later).
--Tip: We can also start new processes when launching Julia from the command line by suing the `-p` command-line argument. E.g., `$ julia -p 3 ` would launch Julia with 3 extra processes. -- -----------Each process runs a separated Julia instance¶
When adding the new processes, you can imagine that 3 new Julia REPLs have started under the hood (see figure below). The main point of the Distributed module is to provide a way of coordinating all these Julia processes to run code in parallel. It is important to note that each process runs in a separated Julia instance. This means that each process has its own memory space and therefore they do not share memory. This results in distributed-memory parallelism, and allows one to run processes in different machines.
--- -
----------Basic info about processes¶
The following functions provide basic information about the underlying processes. If more than one process is available, the first process is called the main or master and the other the workers. If only a single process is available, it is the master and the first worker simultaneously.
- ---- -----In [ ]:------ -procs() -
--- -----In [ ]:------ -workers() -
--- -----In [ ]:------ -nprocs() -
--- -----In [ ]:------ -nworkers() -
--- -----In [ ]:------ -myid() -
---- -----In [ ]:------ -@everywhere println(myid()) -
----------In previous cell, we have used the macro
- -@everywherethat evaluates the given code on all processes. As a result, each process will print its own process id.----------Creating workers in other machines¶
For large parallel computations, one typically needs to use different computers in parallel. Function
-addprocsalso provides a low-level method to start workers in other machines. Next code example would create 3 workers inserver1and 4 new workers in serverserver2(see figure below). Under the hood, Julia connects via ssh to the other machines and starts the new processes there. In order this to work, the local computer and the remote servers need to be properly configured (see the Julia manual for details).- -using Distributed -machines = [("user@server1",3),("user@server2",4)] -addprocs(machines) -
------------ -
-----------Adding workers with ClusterManagers.jl¶
Previous way of starting workers in other machines is very low-level. Happily, there is a Julia package called ClusterManagers.jl that helps to create workers remotely in number of usual scenarios. For instance, when running the following code from the login node in a computer cluster, it will submit a job to the cluster queue allocating 128 threads. A worker will be generated for each one of these threads. If the compute node have 64 cores, 2 compute nodes will be used to create to contain the 128 workers (see below).
-- -using Distributed -using ClusterManagers -addprocs(SlurmManager(128), partition="debug", t="00:5:00") -
------------ -
-----------Executing code remotely¶
We have added new processes to our Julia session. Let's start using them!
- ----------Function
remotecall¶The most basic thing we can do with a remote processor is to execute a given function on it. This is done by using function
- -remotecall. To make clear how local and remote executions compare, let's call a function locally and then remotely. Next cell uses functiononesto create a matrix locally.---- -----In [ ]:------ -a = ones(2,3) -
---------The next cell does the same operation, but remotely on process 2. Note that
- -remotecalltakes the function we want to execute remotely, the process id where we want to execute it and, finally, the function arguments.---- -----In [ ]:------ -proc = 2 -ftr = remotecall(ones,proc,2,3) -
---------Note that
- -remotecalldoes not return the result of the underlying function, but aFuture. This object represents a reference to a task running on the remote process. To move a copy of the result to the current process we can usefetch.---- -----In [ ]:------ -fetch(ftr) -
---------remotecallis asynchronous¶It is important to note that
- -remotecalldoes not wait for the remote process to finish. It turns immediately. This can be checked be calling remotely the following function that sleeps for 10 secods and then generates a matrix.---- -----In [ ]:------ -fun = (m,n) -> (sleep(10); ones(m,n)) -
---------When running next cell, it will return immediately, event though the remote process will sleep for 10 seconds. We can even run code in parallel. To try this execute the second next cell while the remote call is running in the worker.
- ---- -----In [ ]:------ -proc = 2 -ftr = remotecall(fun,proc,2,3) -
---- -----In [ ]:------ -1+1 -
---------However, when fetching the result, the current process blocks waiting until the result is available in the remote process and arrives to its destination.
- ----- -----In [ ]:------ -fetch(ftr) -
---------Useful macro:
@spawnat¶You have provably realized that in order to use
- -remotecallwe have written auxiliary anonymous functions. They are needed to wrap the code we want to execute remotely. Writing these functions can be tedious. Happily, the macro@spawnatgenerates an auxiliary function from the given block of code and callsremotecallfor us. For instance, the two following cells are equivalent.--- -----In [ ]:------ -@spawnat proc ones(2,3) -
---- -----In [ ]:------ -fun = () -> ones(2,3) -remotecall(fun,proc) -
---------@asyncvs@spawnat¶The relation between
- -@asyncand@pawnatis obvious. From the user perspective they work almost in the same way. However,@asyncgenerates a task that runs asynchronously in the current process, whereas@spawnatexecutes a task in a remote process in parallel. In both cases, the result is obtained usingfetch.--- -----In [ ]:------ -tsk = @async begin - sleep(3) - zeros(2) -end -fetch(tsk) -
---- -----In [ ]:------ -ftr = @spawnat :any begin - sleep(3) - zeros(2) -end -fetch(ftr) -
---------Another usefull macro:
@fetchfrom¶Macro
- -@fetchfromis the blocking version of@spawnat. It blocks and returns the corresponding result instead of aFutureobject.---- -----In [ ]:------ -a = @fetchfrom proc begin - sleep(3) - zeros(2) -end -
----------Data movement¶
Data movement is a crucial part in distributed-memory computations and it is usually one of its main computational -bottlenecks. Being aware of the data we are moving when using functions such as
- -remotecallis important to write efficient distributed algorithms in Julia. Julia also provides a special type of channel, called remote channel, to send and receive data between processes.---------Data movement in
remotecall/fetch¶When usig
- -remotecallwe send to the remote process a function and its arguments. In this example, we send function name+and matricesaandbto proc 4. When fetching the result we receive a copy of the matrix from proc 4.---- -----In [ ]:------ -proc = 4 -a = rand(10,10) -b = rand(10,10) -ftr = remotecall(+,proc,a,b) -fetch(ftr); -
---------Implicit data movement¶
Be aware that data movements can be implicit. This usually happens when we execute remotely functions that capture variables. In the following example, we are also sending matrices
- -aandbto proc 4, even though they do not appear as arguments in the remote call. These variables are captured by the anonymous function and will be sent to proc 4.---- -----In [ ]:------ -proc = 4 -a = rand(10,10) -b = rand(10,10) -fun = () -> a+b -ftr = remotecall(fun,proc) -fetch(ftr); -
---------Data movement with remote channels¶
Another way of moving data between processes is to use remote channels. Their usage is very similar to conventional channels for moving data between tasks, but there are some important differences. In the next cell, we create a remote channel. Process 4 puts several values and closes the channel. Like for conventional channels, calls to
- -put!are blocking, but next cell is not blocking the master process since the call toput!runs asynchronously on process 4.--- -----In [ ]:------ -fun = ()->Channel{Int}() -chnl = RemoteChannel(fun) -
---- -----In [ ]:------ -@spawnat 4 begin - for i in 1:5 - put!(chnl,i) - end - close(chnl) -end; -
---------We can take values from the remote channel form any process using
- -take!. Run next cell several times. The sixth time it should raise and error since the channel was closed.---- -----In [ ]:------ -take!(chnl) -
---------Remote channels can be buffered¶
Just like conventional channels, remote channels can be buffered. The buffer is stored in the process that owns the remote channel. By default this corresponds to process that creates the remote channel, but it can be a different one. For instance, process 3 will be the owner in the following example.
- ---- -----In [ ]:------ -buffer_size = 2 -owner = 3 -fun = ()->Channel{Int}(buffer_size) -chnl = RemoteChannel(fun,owner) -
---- -----In [ ]:------ -@spawnat 4 begin - println("start") - for i in 1:5 - put!(chnl,i) - println("I have put $i") - end - close(chnl) - println("stop") -end; -
---------Note that since the channel is buffered, worker 4 can start putting values into it before any call to
- -take!. Run next cell several times until the channel is closed.---- -----In [ ]:------ -take!(chnl) -
---------Remote channels are not iterable¶
One main difference with respect to conventional channels is that remote channels cannot be iterated. Let's repeat the example above.
- ---- -----In [ ]:------ -fun = ()->Channel{Int}() -chnl = RemoteChannel(fun) -
---- -----In [ ]:------ -@spawnat 4 begin - for i in 1:5 - put!(chnl,i) - end - close(chnl) -end; -
---------Now, try to iterate over the channel in a for loop. It will result in an error since channels are not iterable.
- ----- -----In [ ]:------ -for j in chnl - @show j -end -
---------If we want to take values form a remote channel and stop automatically when the channel is closed, we can combine a while loop and a try-catch statement. This works since
- -take!raises an error if the channel is closed, which will execute thecatchblock and breaks the loop.---- -----In [ ]:------ -while true - try - j = take!(chnl) - @show j - catch - break - end -end -
----------Questions¶
----------Q1 (Distributed.ipynb)¶
How many integers are transferred between master and worker? Including both directions.
- ----- -----In [ ]:------ -a = rand(Int,4,4) -proc = 4 -@fetchfrom proc sum(a^2) -
---------Q2 (Distributed.ipynb)¶
How many integers are transferred between master and worker? Including both directions.
- ----- -----In [ ]:------ -a = rand(Int,4,4) -proc = 4 -@fetchfrom proc sum(a[2,2]^2) -
---------Q3 (Distributed.ipynb)¶
Which value will be the value of
- -x?---- -----In [ ]:------ -a = zeros(Int,3) -proc = 3 -@sync @spawnat proc a[2] = 2 -x = a[2] -x -
---------Q4 (Distributed.ipynb)¶
Which value will be the value of
- -x?---- -----In [ ]:------ -a = zeros(Int,3) -proc = myid() -@sync @spawnat proc a[2] = 2 -x = a[2] -x -
----------Remember: each process runs in a separated Julia instance¶
In particular, this means that each process can load different functions or packages. In consequence, it is important to make sure that the code we run is defined in the corresponding process.
- ----------Functions are defined in a single process¶
This is a very common pitfall when running parallel code. If we define a function in a process, it is not automatically available in the other processes. This is illustrated in the next example. The remote call in the last line in next cell will fail since the function
- -sleep_onesis only being defined in the local process.--- -----In [ ]:------ -function sleep_ones(m,n) - sleep(4) - ones(m,n) -end -
---- -----In [ ]:------ -proc = 3 -remotecall_fetch(sleep_ones,proc,3,4) -
---------To fix this, we can define the function on all processes with the
- -@everywheremacro.--- -----In [ ]:------ -@everywhere function sleep_ones(m,n) - sleep(4) - ones(m,n) -end -
---- -----In [ ]:------ -proc = 3 -remotecall_fetch(sleep_ones,proc,3,4) -
---------Anonymous functions are available everywhere¶
If a function has a name, Julia only sends the function name to the corresponding process. Then, Julia looks for the corresponding function code in the remote process and executes it. This is why the function needs to be defined also in the remote process. However, if a function is anonymous, Julia needs to send the complete function definition to the remote process. This is why anonymous functions do not need to be defined with the macro
- -@everywhereto work in a remote call.--- -----In [ ]:------ -fun = (m,n) -> (sleep(4);ones(m,n)) -
---- -----In [ ]:------ -proc = 3 -remotecall_fetch(fun,proc,3,4) -
---------Each proc uses packages independently¶
When using a package in a process, it is not available in the other ones. For instance, if we load the
- -LinearAlgebrapackage in the current process and use one of its exported functions in another process, we will get an error.--- -----In [ ]:------ -using LinearAlgebra -
---- -----In [ ]:------ -@fetchfrom 3 norm([1,2,3]) -
---------To fix this, we can load the package on all processors with the
- -@everywheremacro.--- -----In [ ]:------ -@everywhere using LinearAlgebra -
---- -----In [ ]:------ -@fetchfrom 3 norm([1,2,3]) -
---------Each process has its own active package environment¶
This is another very common source of errors. You can check that if you activate the current directory, this will have no effect in the other processes.
- ----- -----In [ ]:------ -] activate . -
---------We have activated the current folder. Now let's see which is the active project in another process, say process 2. You will see that process 2 is provably still using the global package environment.
- ---- -----In [ ]:------ -@everywhere using Pkg -
---- -----In [ ]:------ -@spawnat 2 Pkg.status(); -
---------To fix this, you need to activate the current directory on all processes.
- ---- -----In [ ]:------ -@everywhere Pkg.activate(".") -
---- -----In [ ]:------ -@spawnat 2 Pkg.status(); -
----------Easy ways of parallelizing code¶
A part from the low-level parallel routines we have seen so-far, Julia also provides much more simple ways to parallelizing loops and maps.
- ----------Useful macro: @distributed¶
This macro is used when we want to perform a very large for loops made of independent small iterations. To illustrate this, let's consider again the function that computes $\pi$ with Leibniz formula.
- ----- -----In [ ]:------ -function compute_π(n) - s = 1.0 - for i in 1:n - s += (isodd(i) ? -1 : 1) / (i*2+1) - end - 4*s -end -
---------Paralelizing this function might require some work with low-level functions like
- -remotecall, but it is trivial using the macro@distributed. This macro runs the for loop using the available processes and optionally reduces the result using a given reduction function (+in this case).---- -----In [ ]:------ -function compute_π_dist(n) - s = 1.0 - r = @distributed (+) for i in 1:n - (isodd(i) ? -1 : 1) / (i*2+1) - end - 4*(s+r) -end -
---------Run next cell to measure the performance of the serial function for a large value of
- -n. Run it at least 2 times to get rid of compilation times.---- -----In [ ]:------ -@time compute_π(4_000_000_000) -
---------Run next cell to measure the performance of the parallel function.
- ----- -----In [ ]:------ -@time compute_π_dist(4_000_000_000) -
---------Useful function:
pmap¶This function is used when we want to call a very expensive function a small number of evaluations and we want to distribute these evaluations over the available processes. To illustrate the usage of
- -pmapconsider the following example. Next cell generates sixty 30x30 matrices. The goal is to compute the singular value decomposition of all of them. This operation is known to be expensive for large matrices. Thus, this is a perfect scenario forpmap.---- -----In [ ]:------ -a = [ rand(300,300) for i in 1:60]; -
---------First, lets measure the serial performance
- ---- -----In [ ]:------ -using LinearAlgebra -
---- -----In [ ]:------ -@time svd.(a); -
---------If we use
- -pmapinstead of broadcast, the different calls tosvdwill be distributed over the available processes.---- -----In [ ]:------ -@time pmap(svd,a); -
---------Summary¶
We have seen the basics of distributed computing in Julia. The programming model is essentially an extension of tasks and channels to parallel computations on multiple machines. The low-level functions are
- -remotecallandRemoteChannel, but there are other functions and macros likepmapand@distributedthat simplify the implementation of parallel algorithms.-- - - - - - - - - diff --git a/docs/src/notebook-output/julia_intro.html b/docs/src/notebook-output/julia_intro.html deleted file mode 100644 index 600f84f..0000000 --- a/docs/src/notebook-output/julia_intro.html +++ /dev/null @@ -1,15629 +0,0 @@ - - - - - --- -----In [ ]:------ --
julia_intro - - - - - - - - - - - - - - - - - - - - - - - -----------Is Julia fast?¶
-----------Contents¶
With this notebook, you will learn
--
-
- Two basic julia concepts related performance:
-
-
- type-inference -
- JIT compilation -
- - Some Julia syntax -
- Some useful Julia packages -
----------Using Jupyter notebooks in Julia¶
We are going to use Jupyter notebooks in this and other lectures. You provably have worked with notebooks (in Python). If not, here are the basic concepts you need to know to follow the lessons.
--Tip: Did you know that Jupyter stands for Julia, Python and R? --How to start a Jupyter nootebook in Julia¶
To run a Julia Jupyther notebook, open a Julia REPL and type
--julia> ] -pkg> add IJulia -julia> using IJulia -julia> notebook() -
A new browser window will open. Navigate to the corresponding notebook and open it.
--Warning: Make sure that the notebook is using the same Julia version as the one you used to launch `IJulia`. If it is not the same, go to Kernel > Change Kernel and choose the right version. -- ----------Running a cell¶
To run a cell, click on a cell and press
- -Shift+Enter. You can also use the "Run" button in the toolbar above.---- -----In [ ]:------ -1+3 -4*5 -
---------As you can see from the output of previous cell, the value of the last line is displayed. We can suppress the output with a semicolon. Try it. Execute next cell.
- ----- -----In [ ]:------ -1+3 -4*5; -
---------Cell order is important¶
Running the two cells below in reverse order won't work (try it).
- ---- -----In [ ]:------ -foo() = "Well done!" -
---- -----In [ ]:------ -foo() -
---------REPL modes¶
This is particular to Julia notebooks. You can use package, help, and shell mode just like in the Julia REPL.
- ---- -----In [ ]:------ -] add MPI -
--- -----In [ ]:------ -? print -
---- -----In [ ]:------ -; ls -
----------How fast is Julia code?¶
-----------NB. Most of the examples below are taken from the lecture by S.G. Johnson at MIT. See here: -https://github.com/mitmath/18S096/blob/master/lectures/lecture1/Boxes-and-registers.ipynb
- -- -----------Hand-written sum function¶
----- -----In [ ]:------ -function sum_hand(a) - s = zero(eltype(a)) - for ai in a - s += ai - end - s -end -
---------Test it¶
The Julia macro
- -@testwhich is provided in theTestpackage is useful to write (unit) tests in Julia.--- -----In [ ]:------ -using Test -
--- -----In [ ]:------ -a = rand(5) -
---- -----In [ ]:------ -@test sum_hand(a) ≈ sum(a) -
---------Benchmarking¶
In Julia, the most straight-forward way of measuring the computation time of a piece of code is with the macro
- -@time.--- -----In [ ]:------ -a = rand(10^7); -
---- -----In [ ]:------ -@time sum_hand(a) -
----------Note that
- -@timealso measures the compile time of a function if it's the first call to that function. So make sure to run@timetwice on a freshly compiled function in order to get a more meaningful result.---------A part of getting rid of compilation time, one typically wants to measure the runtime several times and compute sole. To do this we can call our code in a for-loop and gather the runtimes using the Julia macro
- -@elapsed. This measures the runtime of an expression in seconds, just as the@timemacro, only@elapseddiscards the result of the computation and returns the elapsed time instead.---- -----In [ ]:------ -@elapsed sum_hand(a) -
---------BenchmarkTools¶
The
- -BenchmarkToolsextension package provides useful macros for sampling runtimes automatically.---- -----In [ ]:------ -using BenchmarkTools -
---------First of all, the
- -@benchmarkmacro runs the code multiple times and gives out a lot of details: the minimum and maximum time, mean time, median time, number of samples taken, memory allocations, etc.---- -----In [ ]:------ -bch_sum_hand = @benchmark sum_hand($a) -
---------For quick sanity checks, one can use the
- -@btimemacro, which is a convenience wrapper around@benchmark. It returns only the minimum execution time and memory allocations.---- -----In [ ]:------ -@btime sum_hand($a) -
---------Similar to the
- -@elapsedmacro,BenchmarkTool's@belapseddiscards the return value of the function and instead returns the minimum runtime in seconds.---- -----In [ ]:------ -@belapsed sum_hand($a) -
----------As opposed to
- -@timeand@elapsed,@btimeand@belapsedrun the code several times and return the minimum runtime, thus eliminating possible compilation times from the measurement.---------Built-in sum function¶
----- -----In [ ]:------ -bch_sum = @benchmark sum($a) -
---------Hand-written sum in Python¶
---- -----In [ ]:------ -using PyCall -
--- -----In [ ]:------ -py""" -def sum_py_hand(A): - s = 0.0 - for a in A: - s += a - return s -""" -sum_py_hand = py"sum_py_hand" -
--- -----In [ ]:------ -@test sum(a) ≈ sum_py_hand(a) -
---- -----In [ ]:------ -bch_sum_py_hand = @benchmark sum_py_hand($a) -
---------Numpy sum¶
---- -----In [ ]:------ -using Conda -
--- -----In [ ]:------ -numpy = pyimport("numpy") -sum_numpy = numpy["sum"] -
--- -----In [ ]:------ -@test sum_numpy(a) ≈ sum(a) -
---- -----In [ ]:------ -bch_sum_numpy = @benchmark sum_numpy($a) -
---------Sumary of the results¶
---- -----In [ ]:------ -timings = [bch_sum_hand,bch_sum,bch_sum_py_hand,bch_sum_numpy] -
--- -----In [ ]:------ -methods = ["sum_hand","sum","sum_py_hand","sum_numpy"] -
--- -----In [ ]:------ -using DataFrames -
---- -----In [ ]:------ -df = DataFrame(method=methods,time=timings) -
---------Improving the hand-written sum in Julia¶
---- -----In [ ]:------ -# ✍️ Exercise 3 -function sum_hand_fast(a) - s = 0.0 - @simd for ai in a - s += ai - end - s -end -
--- -----In [ ]:------ -@test sum_hand_fast(a) ≈ sum(a) -
---- -----In [ ]:------ -@benchmark sum_hand_fast($a) -
----------Conlcusions so far¶
-
-
- Julia code (for loops) are much faster than in Python -
- Julia code can be as fast as optimized C code -
----------Why Julia is fast?¶
-
-
- Julia is a compiled language (like C, C++, Fortran) -
- Julia is JIT compiled (C, C++, Fortran are AOT compiled) -
- Type declarations are optional in Julia -
----------Conclusion: Why we use Julia in this course¶
-
-
- Julia code is fast (it can be as fast as C) -
- Julia is a high-level language with simpler syntax than C -
- Julia supports different parallel programming models -
We will look into the third point in a later section of this course.
- ----- -----In [ ]:------ -function sum_hand(a) - s = 0.0 - for ai in a - s += ai - end - s -end -
---------Solution to Exercise 2¶
---- -----In [ ]:------ -using Statistics - -a = rand(10^7) -num_it = 15 -runtimes = zeros(num_it) -for i in 1:num_it - runtimes[i] = @elapsed sum_hand(a) -end -@show mean(runtimes) -@show std(runtimes) -@show minimum(runtimes) -@show maximum(runtimes); -
-- - - - - - - - - diff --git a/docs/src/notebook-output/matrix_matrix.html b/docs/src/notebook-output/matrix_matrix.html deleted file mode 100644 index 20af38f..0000000 --- a/docs/src/notebook-output/matrix_matrix.html +++ /dev/null @@ -1,16058 +0,0 @@ - - - - - --- -----In [ ]:------ -# ✍️ Exercise 3 -function sum_hand_fast(a) - s = 0.0 - @simd for ai in a - s += ai - end - s -end -
matrix_matrix - - - - - - - - - - - - - - - - - - - - - - - -----------
-Programming large-scale parallel systems¶
-----------Matrix-matrix multiplication¶
-----------Contents¶
In this notebook, we will:
--
-
- Parallelize a simple algorithm -
- Study the performance of different parallelization strategies -
- Implement them using Julia -
----------Note: Do not forget to execute the cells below before starting this notebook! -- ---- -----In [ ]:------ -] add BenchmarkTools -
---- -----In [ ]:------ -using Distributed -using BenchmarkTools -using Printf -if procs() == workers() - addprocs(4) -end -function answer_checker(answer,solution) - if answer == solution - "🥳 Well done! " - else - "It's not correct. Keep trying! 💪" - end |> println -end -alg_1_deps_check(answer) = answer_checker(answer,"b") -alg_1_comm_overhead_check(answer) = answer_checker(answer, "c") -alg_1_comp_check(answer) = answer_checker(answer, "a") -alg_2_complex_check(answer) = answer_checker(answer, "b") -alg_2_deps_check(answer) = answer_checker(answer,"d") -alg_3_deps_check(answer) = answer_checker(answer, "c") -alg_3_complex_check(answer) = answer_checker(answer, "d") -
----------Problem Statement¶
Let us consider the (dense) matrix-matrix product
- -C=A*B.------------ -
-----------Goals¶
We want to
--
-
- compute the product in parallel using more than one process (distributed implementation) -
- study the performance of different parallelization alternatives -
- implement the algorithms using Julia -
----------Assumptions¶
-
-
- All matrices
A,B, andCare initially stored in the master process
- - The result will be overwritten in
C
-
------------ -
-----------Steps¶
To develop and study the parallel implementation, we will follow these steps:
--
-
- Identify the parts of the sequential algorithm that can be parallelized -
- Consider different parallelization strategies -
- Discuss the (theoretical) performance of these implementations -
---------Serial implementation¶
We start by considering the (naive) sequential algorithm:
- ----- -----In [ ]:------ -@everywhere function matmul_seq!(C,A,B) - m = size(C,1) - n = size(C,2) - l = size(A,2) - @assert size(A,1) == m - @assert size(B,2) == n - @assert size(B,1) == l - z = zero(eltype(C)) - for j in 1:n - for i in 1:m - Cij = z - for k in 1:l - @inbounds Cij += A[i,k]*B[k,j] - end - C[i,j] = Cij - end - end - C -end -
----------Note: The matrix-matrix multiplication naively implemented with 3 nested loops as above is known to be very inefficient (memory bound). Libraries such as BLAS provide much more efficient implementations, which are the ones used in practice (e.g., by the `*` operator in Julia). We consider, our hand-written implementation as a simple way of expressing the algorithm we are interested in. --Run the following cell to compare the performance of our hand-written function with respect to the built in function
- -mul!---- -----In [ ]:------ -using LinearAlgebra -N = 1000 -A = rand(N,N) -B = rand(N,N) -C = rand(N,N) -@btime matmul_seq!(C,A,B) -@btime mul!(C,A,B); -
----------Where do we can exploit parallelism?¶
Look at the three nested loops in the sequential implementation:
--for j in 1:n - for i in 1:m - Cij = z - for k in 1:l - @inbounds Cij += A[i,k]*B[k,j] - end - C[i,j] = Cij - end -end -
-
-
- Loops over
iandjare trivially parallelizable.
- - The loop over
kcan be parallelized but it requires a reduction.
-
----------Parallel algorithms¶
All the entries of matrix C can be potentially computed in parallel, but is it the most efficient solution to solve all these entries in parallel in a distributed system? To find this we will consider different parallelization strategies:
--
-
- Algorithm 1: each worker computes a single entry of C -
- Algorithm 2: each worker computes a single row of C -
- Algorithm 3: each worker computes a block rows of C -
------------ -
-----------Parallel algorithm 1¶
-----------Data dependencies¶
Moving data through the network is expensive and reducing data movement is one of the key points in a distributed algorithm. To this end, we determine which is the minimum data needed by a worker to perform its computations.
-In algorithm 1, each worker computes only an entry of the result matrix C.
- ------------- -
-----------Question: Which are the data dependencies of the computations done by the worker in charge of computing entry C[i,j] ? --
- -a) column A[:,i] and row B[j,:] -b) row A[i,:] and column B[:,j] -c) the whole matrices A and B -d) row A[i,:] and the whole matrix B ----- -----In [ ]:------ -answer = "x" # replace x with a, b, c, or d -alg_1_deps_check(answer) -
----------Implementation¶
Taking into account the data dependencies, the parallel algorithm 1 can be efficiently implemented following these steps from the worker perspective:
--
-
- The worker receives the corresponding row A[i,:] and column B[:,j] from the master process -
- The worker computes the dot product of A[i,:] and B[:,j] -
- The worker sends back the result of C[i,j] to the master process -
------------ -
----------A possible implementation of this algorithm in Julia is as follows:
- ----- -----In [ ]:------ -function matmul_dist_1!(C, A, B) - m = size(C,1) - n = size(C,2) - l = size(A,2) - @assert size(A,1) == m - @assert size(B,2) == n - @assert size(B,1) == l - z = zero(eltype(C)) - @assert nworkers() == m*n - iw = 0 - @sync for j in 1:n - for i in 1:m - Ai = A[i,:] - Bj = B[:,j] - iw += 1 - w = workers()[iw] - ftr = @spawnat w begin - Cij = z - for k in 1:l - @inbounds Cij += Ai[k]*Bj[k] - end - Cij - end - @async C[i,j] = fetch(ftr) - end - end - C -end -
---------You can execute the following cells to test this implementation.
- ----- -----In [ ]:------ -using Test -N = 2 -A = rand(N,N) -B = rand(N,N) -C = similar(A) -@test matmul_dist_1!(C,A,B) ≈ A*B -
----------Performance¶
Let us study the performance of this algorithm. To this end, we will analyze if algorithm 1 is able to achieve the optimal parallel speedup. The parallel speedup on $P$ processes is defined as
-$$ -S_P = \frac{T_1}{T_P}, -$$where $T_1$ denotes the runtime of the sequential algorithm on one node and $T_P$ denotes the runtime of the parallel algorithm on $P$ processes. If we run an optimal parallel algorithm with $P$ processes we expect it to run $p$ times faster than the sequential implementation. I.e., the optimal speedup of a parallel algorithm on $p$ processes is equal to $P$:
-$$ -S^{*}_p = P. -$$The ratio of the actual speedup over the optimal one is called the parallel efficiency
-$$ -E_p = \frac{S_p}{S^{*}_p} = \frac{T_1/T_P}{P}. -$$ ----------Experimental speedup¶
The following cell measures the speedup of parallel algorithm 1. Do we achieve the optimal speedup?
- ----- -----In [ ]:------ -N = 2 -A = rand(N,N) -B = rand(N,N) -C = similar(A) -T1 = @belapsed matmul_seq!(C,A,B) -C = similar(A) -TP = @belapsed matmul_dist_1!(C,A,B) -P = nworkers() -println("Speedup = ", T1/TP) -println("Optimal speedup = ", P) -println("Efficiency = ", 100*(T1/TP)/P, "%") -
---------Communication overhead¶
Since communication is usually the main bottleneck in a distributed algorithm, we want to reduce the amount of communication per unit of computation in a worker. Let us compute the (theoretical) communication overhead for algorithm 1. This will help us understand why the speedup of this algorithm was so bad.
-Remember, algorithm 1 consisted of these main steps:
--
-
- The worker receives the corresponding row A[i,:] and column B[:,j] from the master process -
- The worker computes the dot product of A[i,:] and B[:,j] -
- The worker sends back the result of C[i,j] to the master process -
-Question: How many scalars are communicated from and to a worker? Assume that matrices A, B, and C are N by N matrices. --
- -a) 3N -b) 2N + 2 -c) 2N + 1 -d) N² + 1 ----- -----In [ ]:------ -answer = "x" # replace x with a, b, c, or d -alg_1_comm_overhead_check(answer) -
----------Question: How many operations are done in a worker? --
- -a) O(N) -b) O(N²) -c) O(N³) ----- -----In [ ]:------ -answer = "x" # replace x with a, b, or c -alg_1_comp_check(answer) -
----------From these results we can conclude:
--
-
- The communication complexity is O(N) -
- The computation complexity is O(N) -
- The ratio communication over computation (the communication overhead) is O(1) -
In other words, the communication cost is of the same order of magnitude as the computation cost. Since, communication is orders of magnitude slower in real systems, the runtime in the worker will be dominated by communication. This explains why we obtained such a bad speedup.
- -----------Parallel algorithm 2¶
Let us study the next algorithm to see if we can improve the efficiency by augmenting the granularity (i.e. the amount of work) in each parallel task. In parallel algorithm 2, each worker computes an entire row of C.
- ------------- -
-----------Data dependencies¶
-----------Question: Which are the data dependencies of the computations done by the worker in charge of computing row C[i,:] ? --
- -a) column A[:,i] and row B[j,:] -b) row A[i,:] and column B[:,j] -c) the whole matrices A and B -d) row A[i,:] and the whole matrix B ----- -----In [ ]:------ -answer = "x" # replace x with a, b, c, or d -alg_2_deps_check(answer) -
----------Implementation¶
These are the main steps of the implementation of algorithm 2:
--
-
- The worker receives the corresponding row A[i,:] and matrix B from the master process -
- The worker computes the product of row A[i,:] times B -
- The worker sends back the result of row C[i,:] to the master process -
------------ -
----------A possible implementation of this algorithm in Julia is as follows:
- ----- -----In [ ]:------ -function matmul_dist_2!(C, A, B) - m = size(C,1) - n = size(C,2) - l = size(A,2) - @assert size(A,1) == m - @assert size(B,2) == n - @assert size(B,1) == l - z = zero(eltype(C)) - @assert nworkers() == m - iw = 0 - @sync for i in 1:m - Ai = A[i,:] - iw += 1 - w = workers()[iw] - ftr = @spawnat w begin - Ci = fill(z,l) - for j in 1:n - for k in 1:l - @inbounds Ci[j] += Ai[k]*B[k,j] - end - end - Ci - end - @async C[i,:] = fetch(ftr) - end - C - end -
---------Test it using next cell
- ----- -----In [ ]:------ -using Test -N = 4 -A = rand(N,N) -B = rand(N,N) -C = similar(A) -@test matmul_dist_2!(C,A,B) ≈ A*B -
---------Experimental speedup¶
----- -----In [ ]:------ -N = 4 -A = rand(N,N) -B = rand(N,N) -C = similar(A) -T1 = @belapsed matmul_seq!(C,A,B) -C = similar(A) -TP = @belapsed matmul_dist_2!(C,A,B) -P = nworkers() -println("Speedup = ", T1/TP) -println("Optimal speedup = ", P) -println("Efficiency = ", 100*(T1/TP)/P, "%") -
---------Complexity¶
The speedup is still far from the optimal one. Let us study the communication overhead for this algorithm. Remember, algorithm 2 consists in these main steps:
--
-
- The worker receives the corresponding row A[i,:] and matrix B from the master process -
- The worker computes the product of row A[i,:] times B -
- The worker sends back the result of row C[i,:] to the master process -
-Question: Which is the complexity of the communication and computations done by a worker in algorithm 2? --
- -a) O(N) communication and O(N^2) computation -b) O(N^2) communication and O(N^2) computation -c) O(N^3) communication and O(N^3) computation -d) O(N) communication and O(N) computation ----- -----In [ ]:------ -answer = "x" # replace x with a, b, c, or d -alg_2_complex_check(answer) -
----------The communication and computation cost are still of the same order of magnitude even though we have increased the grain size.
- -----------Parallel algorithm 3¶
Let us increase even more the granularity of the parallel tasks by computing several rows of C in a worker. Each worker computes N/P consecutive rows of C, where P is the number of workers.
- ------------- -
-----------Data dependencies¶
-----------Question: Which are the data dependencies of the computations done by the worker in charge of computing the range of rows C[rows,:] ? --
- -a) A[rows,:] and B[:,rows] -b) the whole matrix A and B[:,rows] -c) A[rows,:] and the whole matrix B -d) the whole matrices A and B ----- -----In [ ]:------ -answer = "x" # replace x with a, b, c, or d -alg_3_deps_check(answer) -
----------Implementation¶
These are the main steps of the implementation of algorithm 3:
--
-
- The worker receives the corresponding rows A[rows,:] and matrix B from the master process -
- The worker computes the product of A[rows,:] times B -
- The worker sends back the result of C[rows,:] to the master process -
------------ -
-----------The implementation of this variant is let as an exercise (see below).
- ----------Communication overhead¶
Let us analyze the (theoretical) communication overhead for algorithm 3.
--Question: Which is the complexity of the communication and computations done by a worker in algorithm 3? --
- -a) O(N²) communication and O(N³) computation -b) O(N²) communication and O(N³/P) computation -c) O(N²) communication and O(N²) computation -d) O(N²/P) communication and O(N³/P) computation ----- -----In [ ]:------ -answer = "x" # replace x with a, b, c, or d -alg_3_complex_check(answer) -
----------In this case, the ratio between communication and computation is O(P/N). If the matrix size N is much larger than the number of workers P, then the communication overhead O(P/N) would be negligible. This opens the door to an scalable implementation.
- -----------Summary¶
The table below compares the three parallel algorithms.
- -- -
-- - - -Algorithm -Parallelism -
(#workers)Communication -
per workerComputation -
per workerRatio communication/ -
computation- -1 -N² -2N + 1 -N -O(1) -- -2 -N -2N + N² -N² -O(1) -- - -3 -P -N² + 2N²/P -N³/P -O(P/N) --
-
- Matrix-matrix multiplication is trivially parallelizable (all entries in the result matrix can be computed in parallel, at least in theory) -
- However, we cannot exploit all the potential parallelism in a distributed system due to communication overhead -
- We need a sufficiently large grain size to obtain a near optimal speedup -
----------Exercises¶
----------Exercise 1¶
Implement algorithm 3 in the function below. For simplicity, assume that the number of rows of C is a multiple of the number of workers.
- ----- -----In [ ]:------ -function matmul_dist_3!(C,A,B) - m = size(C,1) - n = size(C,2) - l = size(A,2) - @assert size(A,1) == m - @assert size(B,2) == n - @assert size(B,1) == l - @assert mod(m,nworkers()) == 0 - # Implement here - - C -end -
---------Use test-driven development to implement the algorithm. Use this test:
- ----- -----In [ ]:------ -using Test -P = nworkers() -load = 100 -N = load*P -A = rand(N,N) -B = rand(N,N) -C = similar(A) -@test matmul_dist_3!(C,A,B) ≈ A*B -
---------Measure the performance of your implementation by running next cell. Do you get close to the optimal speedup?
- ----- -----In [ ]:------ -P = nworkers() -load = 100 -N = load*P -A = rand(N,N) -B = rand(N,N) -C = similar(A) -T1 = @belapsed matmul_seq!(C,A,B) -C = similar(A) -TP = @belapsed matmul_dist_3!(C,A,B) -println("Speedup = ", T1/TP) -println("Optimal speedup = ", P) -println("Efficiency = ", 100*(T1/TP)/P, "%") -
----------Exercise 2¶
----------The implementation of algorithm 1 is very impractical. One needs as many processors as entries in the result matrix C. For 1000 times 1000 matrix one would need a supercomputer with one million processes! We can easily fix this problem by using less processors and spawning the computation of an entry in any of the available processes. -See the following code:
- ----- -----In [ ]:------ -function matmul_dist_1_v2!(C, A, B) - m = size(C,1) - n = size(C,2) - l = size(A,2) - @assert size(A,1) == m - @assert size(B,2) == n - @assert size(B,1) == l - z = zero(eltype(C)) - @sync for j in 1:n - for i in 1:m - Ai = A[i,:] - Bj = B[:,j] - ftr = @spawnat :any begin - Cij = z - for k in 1:l - @inbounds Cij += Ai[k]*Bj[k] - end - Cij - end - @async C[i,j] = fetch(ftr) - end - end - C -end -
---------With this new implementation, we can multiply matrices of arbitrary size with a fixed number of workers. Test it:
- ----- -----In [ ]:------ -using Test -N = 50 -A = rand(N,N) -B = rand(N,N) -C = similar(A) -@test matmul_dist_1_v2!(C,A,B) ≈ A*B -
---------Run the next cell to check the performance of this implementation. Note that we are far away from the optimal speed up. Why? To answer this question compute the theoretical communication over computation ratio for this implementation and reason about the obtained result. Hint: the number of times a worker is spawned in this implementation is N^3/P on average.
- ---- -----In [ ]:------ -N = 100 -A = rand(N,N) -B = rand(N,N) -C = similar(A) -P = nworkers() -T1 = @belapsed matmul_seq!(C,A,B) -C = similar(A) -TP = @belapsed matmul_dist_1_v2!(C,A,B) -println("Speedup = ", T1/TP) -println("Optimal speedup = ", P) -println("Efficiency = ", 100*(T1/TP)/P, "%") -
-- - - - - - - - - diff --git a/docs/src/notebook-output/notebook-hello.html b/docs/src/notebook-output/notebook-hello.html deleted file mode 100644 index 5a41d45..0000000 --- a/docs/src/notebook-output/notebook-hello.html +++ /dev/null @@ -1,14832 +0,0 @@ - - - - - --- -----In [ ]:------ --
notebook-hello - - - - - - - - - - - - - - - - - - - - - - - ----------Hello, world!¶
---- -----In [ ]:------ -] add Distributed -
--- -----In [ ]:------ -] add https://github.com/fverdugo/W2bYx4.git -
--- -----In [ ]:------ -using W2bYx4 -
---- -----In [ ]:------ -W2bYx4.greet() -
----------This is an attached picture:
- ------------- -
-----------Another attached picture with markdown syntax
-![]():
- -----------A picture embedded with Base64:
--- -
-----------This is a linked picture:
- ------------- -
----------- -julia> 1 + 1 -
---- -----In [ ]:------ -1+1 -
-----------Warning: Make sure that the file- -"hello.jl"is located in the current working directory of your Julia session. -----------Tip: Avoid calling Julia code from the terminal -- --- - - - - - - - - diff --git a/docs/src/notebook-output/sol_matrix_matrix.html b/docs/src/notebook-output/sol_matrix_matrix.html deleted file mode 100644 index 4d8dc4a..0000000 --- a/docs/src/notebook-output/sol_matrix_matrix.html +++ /dev/null @@ -1,14652 +0,0 @@ - - - - - --- -----In [ ]:------ --
sol_matrix_matrix - - - - - - - - - - - - - - - - - - - - - - - --- - - - - - - - - diff --git a/docs/src/sol_matrix_matrix.md b/docs/src/sol_matrix_matrix.md deleted file mode 100644 index 770dc55..0000000 --- a/docs/src/sol_matrix_matrix.md +++ /dev/null @@ -1,30 +0,0 @@ -```@meta -EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/sol_matrix_matrix.ipynb" -``` - -```@raw html --- -----In [ ]:------ -function matmul_dist_3!(C,A,B) - m = size(C,1) - n = size(C,2) - l = size(A,2) - @assert size(A,1) == m - @assert size(B,2) == n - @assert size(B,1) == l - @assert mod(m,nworkers()) == 0 - # Implement here - nrows_w = div(m,nworkers()) - @sync for (i,w) in enumerate(workers()) - rows_w = (1:nrows_w) .+ (i-1)*nrows_w - Aw = A[rows_w,:] - ftr = @spawnat w begin - Cw = similar(Aw,nrows_w,n) - matmul_seq!(Cw,Aw,B) - Cw - end - @async C[rows_w,:] = fetch(ftr) - end - C -end -
--``` - -```@raw html - - -``` diff --git a/notebooks/jacobi_method.ipynb b/notebooks/jacobi_method.ipynb index 7e95ef0..019359c 100644 --- a/notebooks/jacobi_method.ipynb +++ b/notebooks/jacobi_method.ipynb @@ -33,6 +33,38 @@ "\n" ] }, + { + "cell_type": "code", + "execution_count": 5, + "id": "1dc78750", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "jacobi_2_check (generic function with 1 method)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using Printf\n", + "\n", + "function answer_checker(answer,solution)\n", + " if answer == solution\n", + " \"🥳 Well done! \"\n", + " else\n", + " \"It's not correct. Keep trying! 💪\"\n", + " end |> println\n", + "end\n", + "gauss_seidel_1_check(answer) = answer_checker(answer,\"c\")\n", + "jacobi_1_check(answer) = answer_checker(answer, \"d\")\n", + "jacobi_2_check(answer) = answer_checker(answer, \"b\")" + ] + }, { "cell_type": "markdown", "id": "d4cb59d5", @@ -244,12 +276,23 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "4edad93f", - "metadata": {}, - "outputs": [], + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "It's not correct. Keep trying! 💪\n" + ] + } + ], "source": [ - "#TODO answer (c)" + "answer = \"x\" # replace x with a, b, c or d\n", + "gauss_seidel_1_check(answer)" ] }, { @@ -304,7 +347,7 @@ "source": [ "### Ghost (aka halo) cells\n", "\n", - "A usual way of handling this type of data dependencies is using so-called ghost cells. Ghost cells represent the missing data dependencies in the data owned by each process. After importing the appropriate values from the neighbor processes one can perform the usual sequential jacoby update locally in the processes." + "A usual way of handling this type of data dependencies is using so-called ghost cells. Ghost cells represent the missing data dependencies in the data owned by each process. After importing the appropriate values from the neighbor processes one can perform the usual sequential Jacobi update locally in the processes." ] }, { @@ -324,17 +367,31 @@ "source": [ "Tip - -\n", "Question: Which is the communication and computation complexity in each process? N is the length of the vector and P the number of processes. \n", - "" + "