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.
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
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
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.
# analytical solution
-functionanalytical_solution(N)
-# Returns the analytical solution as a square grid of size N
-grid=zeros(N,N)
-foriin2:N
-grid[:,i].=(i-1)/(N-1)
-end
-returngrid
-end
-
-# Test if solution is identical with analytical solution
-sol=analytical_solution(N)
-@testmaximum(abs.(sol-new_grid[2:M-1,:]))<0.01*N
-
# to import MPIManager
-usingMPIClusterManagers
-
-# need to also import Distributed to use addprocs()
-usingDistributed
-
-# specify, number of mpi workers, launch cmd, etc.
-manager=MPIWorkerManager(9)
-
-# start mpi workers and add them as julia workers too.
-addprocs(manager)
-
-@mpi_domanagerbegin
-
-functioncalculate_partition(p,N,nrows,ncols)
-# Calculates the row and column indices of this processor
-# Get row and column number for processor p
-ifmod(p,ncols)==0
-i=div(p,ncols)
-else
-i=floor(div(p,ncols))+1
-end
-j=p-(i-1)*ncols
-# Rows
-ifmod(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))
-ifi<=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
-ifmod(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))
-ifj<=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
-returnrow_range,col_range
-end
-
-
-functionupdate_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
-return0.25*(g_up+g_down+g_left+g_right)
-end
-
-usingMPI
-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
-@assertn_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
-ifid==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)
-localgrid_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
-ifmod(id,ncols)==1
-# Left neighbor is last in row
-n_left=id+ncols-1
-end
-ifmod(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)")
-
-localfinished=false
-
-#Perform SOR
-while!finished
-# Flip old and new grid
-grid_a_old=!grid_a_old
-
-# Determine which grid is updated
-ifgrid_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
-elseiftop_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
-elseifbottom_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
-ifid==1
-ifall(max_diffs.<ϵ)
-finished=true
-println("THRESHOLD SUBCEEDED - TERMINATE SOR")
-end
-end
-
-finished=MPI.bcast(finished,0,comm)
-
-iffinished
-# 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
-
-usingDelimitedFiles
-
-# 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
-ifid==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
-
# Test if solution is identical with analytical solution
-sol=analytical_solution(M)
-# Bring solution into correct form
-sol=reverse(transpose(sol),dims=1)
-@testmaximum(abs.(sol-final_grid[:,2:N-1]))<0.01*M
-
The Jacobi method is a numerical tool to solve systems of linear algebraic equations. One of the main applications of the method is to solve boundary value problems (BVPs). I.e., given the values at the boundary (of a grid), the Jacobi method will find the interior values that fulfill a certain equation.
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
When solving a Laplace equation in 1D, the Jacobi method leads to this iterative scheme. The entry $i$ of vector $u$ at iteration $t+1$ is computed as:
-Note: In our version of the jacobi method, we return after a given number of iterations. Other stopping criteria are possible. For instance, iterate until the difference between u and u_new is below a tolerance.
-
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.
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-Question: Which is the communication and computation complexity in each process? N is the length of the vector and P the number of processes.
-
We consider the implementation using MPI. The programming model of MPI is generally better suited for data-parallel algorithms like this one than the task-based model provided by Distributed.jl. In any case, one can also implement it using Distributed, but it requires some extra effort to setup remote channels right for the communication between neighbor processes.
-
Take a look at the implementation below and try to understand it. Note that we have used MPIClustermanagers and Distributed just to run the MPI code on the notebook. When running it on a cluster, MPIClustermanagers and Distributed are not needed.
-Question: How many messages per iteration are sent from a process away from the boundary?
-
-
a) 1
-b) 2
-c) 3
-d) 4
-
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
answer="x"# replace x with a, b, c or d
-jacobi_2_check(answer)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-Question: After the end of the for-loop (line 43), ...
-
-
a) each worker holds the complete solution.
-b) the root process holds the solution.
-c) the ghost cells contain redundant values.
-d) all ghost cells contain the initial values -1 and 1.
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
answer="x"# replace x with a, b, c or d
-jacobi_3_check(answer)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-Question: In line 35 of the code, we wait for all receive and send requests. Is it possible to instead wait for just the receive requests?
-
-
a) No, because the send buffer might be overwritten if we don't wait for send requests.
-b) No, because MPI does not allow an asynchronous send without a Wait().
-c) Yes, because each send has a matching receive, so all requests are done when the receive requests return.
-d) Yes, because there are no writes to the send buffer in this iteration.
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
answer="x"# replace x with a, b, c or d.
-jacobi_4_check(answer)
-
Note that we only need communications to update the values at the boundary of the portion owned by each process. The other values (the one in green in the figure below) can be updated without communications. This provides the opportunity of overlapping the computation of the interior values (green cells in the figure) with the communication of the ghost values. This technique is called latency hiding, since we are hiding communication latency by overlapping it with communications that we need to do anyway.
-
The modification of the implementation above to include latency hiding is leaved as an exercise (see below).
For the Laplace equation in 2D, the interior values in the computational grid (represented by a matrix $u$) are computed with this iterative scheme. The entry $(i,j)$ of matrix $u$ at iteration $t+1$ is computed as:
Transform the following parallel implementation of the 1d Jacobi method (it is copied from above) to use latency hiding (overlap between computation of interior values and communication)
Compute the complexity of the communication and computation of the three data partition strategies (1d block partition, 2d block partition, and 2d cyclic partition) when computing a single iteration of the Jacobi method in 2D. Assume that the grid is of size $N \times N$ and the number of processes $P$ is a perfect square number, i.e. $\sqrt{P}$ is an integer. Hint: For the complexity analysis, you can ignore the effect of the boundary conditions.
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.
The task has been created, but the corresponding work has not started. Note that we do not see any output from function work yet. To run the task we need to schedule it.
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.
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.
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.
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 @async macro. 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.
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.
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.
Note that put! and take! are blocking operations. Calling put! blocks the tasks until another task calls take! and viceversa. Thus, we need at least 2 tasks for this to work. If we call put! and take! 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 how put! blocks until we call take!.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
chnl=Channel{Int}()
-
-
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
@asyncbegin
-foriin1:5
-put!(chnl,i)
-println("I have put $i")
-end
-close(chnl)
-end
-
We can be a bit more flexible and use a buffered channel. In this case, put! will block only if the channel is full and take! will block if the channel is empty. We repeat previous example, but with a buffered channel of size 2. Note that we can call put! until the channel is full. At this point, we need to wait to until we call take! which removes an item from the channel, making room for a new item.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
buffer_size=2
-chnl=Channel{Int}(buffer_size)
-
-
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
@asyncbegin
-foriin1:5
-put!(chnl,i)
-println("I have put $i")
-end
-close(chnl)
-end
-
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?
-
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.
-
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.
The usage of variables in Julia is pretty similar to Python and quite different from C/C++ and Fortran. However, there are also some differences with Python.
A variable is a name associated (bound) to a value. We associate variables with values with = as usual.
-
-
-
-
-
-
-
-
-
-
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).
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.
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.
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.
There is a great flexibility to choose variable names in Julia. See all the rules in the manual if you are interested.
-
-
-
-
-
-
-
-
-
-
-
-
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)
-
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.
Functions are defined as shown in next cell. The closing end is necessary. Do to forget it! However, the return is optional. The value of last line is returned by default. Indentation is recommended, but it is also optional. That's why the closing end is needed.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
functionadd(a,b)
-returna+b
-end
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Once defined, a function can be called using bracket notation as you would expect.
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 ?
-
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_anonymous is 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).
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.
Functions with multiple arguments are also supported. The following example iterates over the given arguments and prints them. args is just a tuple with all arguments.
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.
Arrays of fixed element type seem to be very rigid, right? Python list have not this limitation. However, we can use arrays of Any type, which are as flexible as Python lists, or even more since they can also contain functions.
Implement a function ex1(a) that finds the largest item in the array a. 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.
Implement a function ex2(f,g) that takes two functions f(x) and g(x) and returns a new function h(x) representing the sum of f and g, i.e., h(x)=f(x)+g(x).
If the value of mandel is less than max_iters, the point is provably outside the Mandelbrot set. If mandel is equal to max_iters, then the point is provably inside the set. The larger max_iters, the better the quality of the estimate (the nicer will be your plot).
-
Plot the value of function mandel for each pixel in a 2D grid of the box.
-$$(-1.7,0.7)\times(-1.2,1.2).$$
Use a grid resolution of at least 1000 points in each direction and max_iters at least 10. You can increase these values to get nicer plots. To plot the values use function heatmap from the Julia package GLMakie. Use LinRange to divide the horizontal and vertical axes into pixels. See the documentation of these functions for help. GLMakie is 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.
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 [ ]:
-
-
-
usingPrintf
-functionanswer_checker(answer,solution)
-ifanswer==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")
-
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.
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 [ ]:
-
-
-
usingDistributed
-
-
-
-
-
-
-
-
-
-
-
-
-
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.
-
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.
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 [ ]:
-
-
-
@everywhereprintln(myid())
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
In previous cell, we have used the macro @everywhere that evaluates the given code on all processes. As a result, each process will print its own process id.
For large parallel computations, one typically needs to use different computers in parallel. Function addprocs also provides a low-level method to start workers in other machines. Next code example would create 3 workers in server1 and 4 new workers in server server2 (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).
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).
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 function ones to create a matrix locally.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
a=ones(2,3)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
The next cell does the same operation, but remotely on process 2. Note that remotecall takes 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 remotecall does not return the result of the underlying function, but a Future. 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 use fetch.
It is important to note that remotecall does 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.
You have provably realized that in order to use remotecall we 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 @spawnat generates an auxiliary function from the given block of code and calls remotecall for us. For instance, the two following cells are equivalent.
The relation between @async and @pawnat is obvious. From the user perspective they work almost in the same way. However, @async generates a task that runs asynchronously in the current process, whereas @spawnat executes a task in a remote process in parallel. In both cases, the result is obtained using fetch.
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 remotecall is 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.
When usig remotecall we send to the remote process a function and its arguments. In this example, we send function name + and matrices a and b to proc 4. When fetching the result we receive a copy of the matrix from proc 4.
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 a and b to 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.
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 to put! runs asynchronously on process 4.
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.
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.
@spawnat4begin
-println("start")
-foriin1: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.
Now, try to iterate over the channel in a for loop. It will result in an error since channels are not iterable.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
forjinchnl
-@showj
-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 the catch block and breaks the loop.
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.
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_ones is only being defined in the local process.
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 @everywhere to work in a remote call.
When using a package in a process, it is not available in the other ones. For instance, if we load the LinearAlgebra package in the current process and use one of its exported functions in another process, we will get an error.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
usingLinearAlgebra
-
-
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
@fetchfrom3norm([1,2,3])
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
To fix this, we can load the package on all processors with the @everywhere macro.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
@everywhereusingLinearAlgebra
-
-
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
@fetchfrom3norm([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 [ ]:
-
-
-
@everywhereusingPkg
-
-
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
@spawnat2Pkg.status();
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
To fix this, you need to activate the current directory on all processes.
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.
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).
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 pmap consider 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 for pmap.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
a=[rand(300,300)foriin1:60];
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
First, lets measure the serial performance
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
usingLinearAlgebra
-
-
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
@timesvd.(a);
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
If we use pmap instead of broadcast, the different calls to svd will be distributed over the available processes.
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 remotecall and RemoteChannel, but there are other functions and macros like pmap and @distributed that simplify the implementation of parallel algorithms.
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?
-
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.
-
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 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 [ ]:
-
-
-
@timesum_hand(a)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Note that @time also measures the compile time of a function if it's the first call to that function. So make sure to run @time twice 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 @time macro, only @elapsed discards the result of the computation and returns the elapsed time instead.
The BenchmarkTools extension package provides useful macros for sampling runtimes automatically.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
usingBenchmarkTools
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
First of all, the @benchmark macro 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=@benchmarksum_hand($a)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
For quick sanity checks, one can use the @btime macro, which is a convenience wrapper around @benchmark. It returns only the minimum execution time and memory allocations.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
@btimesum_hand($a)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Similar to the @elapsed macro, BenchmarkTool's @belapsed discards the return value of the function and instead returns the minimum runtime in seconds.
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
@belapsedsum_hand($a)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
As opposed to @time and @elapsed, @btime and @belapsed run the code several times and return the minimum runtime, thus eliminating possible compilation times from the measurement.
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.
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+P to open the command bar, and execute Julia: 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.
Curious about what function println does? 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 to help?>. It means we are in help mode. Now, we can type a function name to see its documentation.
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.
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 name hello.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 function pwd(). You can change to another directory with function cd() 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 PATH environment 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.
-
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>usingDistributed
-julia>@everywhereprintln("Hello, world! I am proc $(myid()) from $(nprocs())")
-
-
Here, we are using the Distributed package, 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 addprocs function.
-
julia>addprocs(3)
-
-
We have added 3 new processes, plus the old one, we have 4 processes. Run the code again.
-
julia>@everywhereprintln("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 -p argument (useful, e.g., when running in a cluster). If you launch Julia from the terminal as
-
$ julia -p 3
-
-
and then run
-
julia>@everywhereprintln("Hello, world! I am proc $(myid()) from $(nprocs())")
-
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
-usingMPI
-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>addMPI
-
-
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 .jl needs to be ommited. In this case, we have isntalled the MPI.jl package even though we have only typed MPI in 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 of MPI.jl for 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.
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 mpiexec needs to match with the MPI version we are using from Julia. You can find the path to the mpiexec binary you need to use with these commands
-
julia>usingMPI
-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>usingMPI
-julia>mpiexec(cmd->run(`$cmd -np 4 julia hello_mpi.jl`))
-
We have installed the MPI package 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 activate command 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 the lessons folder. 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 --project flag. 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 variable JULIA_PROJECT with 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 folderB and then julia> cd("folderA"), will activate the project in folderB and change the current working directory to folderA.
-
-
At this point all package-related operations will be local to the new project. For instance, install the DataFrames package.
Now, we can return to the global project to check that DataFrames has not been installed there. To return to the global environment, use activate without a folder name.
The information about a project is stored in two files Project.toml and Manifest.toml.
-
-
Project.toml contains the packages explicitly installed (the direct dependencies)
-
-
Manifest.toml contains direct and indirect dependencies along with the concrete version of each package.
-
-
-
In other words, Project.toml contains the packages relevant for the user, whereas Manifest.toml is the detailed snapshot of all dependencies. The Manifest.toml can be used to reproduce the same envinonment in another machine.
-
You can see the path to the current Project.toml file by using the status operator (or st in its short form) while in package mode
-
(@v1.8)pkg>status
-
-
The information about the Manifest.toml can be inspected by passing the -m flag.
Copy the contents of previous code block into a file called Project.toml and place it in an empty folder named newproject. It is important that the file is named Project.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.toml file
-
(@v1.8)pkg>activatenewproject
-
-
and then instantiating it
-
(newproject)pkg>instantiate
-
-
The instantiate command will download and install all listed packages and their dependencies in just one click.
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 Pkg package. For instance
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
-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!
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
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)
-
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
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.
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.
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.
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.
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.
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:
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.
Message Passing Interface (MPI) is a standardized and portable library specification for communication between parallel processes in distributed memory systems. Julia offers a convenient way to work with MPI for creating efficient parallel and distributed applications. In this tutorial, you will learn how to use MPI from Julia to perform parallel computing tasks.
When you run an MPI-enabled Julia script, MPI takes care of spawning multiple instances of the Julia executable, each acting as a separate process. These workers can communicate with each other using MPI communication functions. This enables parallel processing and distributed computation. Here's a summary of how it works:
-
-- TODO: insert picture here --
-
-
MPI Spawns Processes: The mpiexec command launches multiple instances of the Julia executable, creating separate worker processes. In this example, 4 Julia workers are spawned.
-
-
Worker Communication: These workers can communicate with each other using MPI communication functions, allowing them to exchange data and coordinate actions.
-
-
Parallel Tasks: The workers execute parallel tasks simultaneously, working on different parts of the computation to potentially speed up the process.
-
-
-
Installing MPI.jl and MPIClusterManagers Packages¶
To use MPI in Julia, you'll need the MPI.jl package, and if you intend to run MPI programs in a Jupyter Notebook, you'll also need the MPIClusterManagers package. These packages provide the necessary bindings to the MPI library and cluster management capabilities. To install the packages, open a terminal and run the following commands:
- Tip:
-The package MPI.jl is the Julia interface to MPI. 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 of MPI.jl for further details.
-
Let's start by creating a simple MPI program that prints a message along with the rank of each worker.
-
Create a new Julia script, for example, mpi_hello_world.jl:
-
usingMPI
-
-# Initialize MPI
-MPI.Init()
-
-# Get the default communicator (MPI_COMM_WORLD) for all processes
-comm=MPI.COMM_WORLD
-
-# Get the number of processes in this communicator
-nranks=MPI.Comm_size(comm)
-
-# Get the rank of the current process within the communicator
-rank=MPI.Comm_rank(comm)
-
-# Print a message with the rank of the current process
-println("Hello, I am process $rank of $nranks processes!")
-
-# Finalize MPI
-MPI.Finalize()
-
In MPI, a communicator is a context in which a group of processes can communicate with each other. MPI_COMM_WORLD is one of the MPI standard communicators, it represents all processes in the MPI program. Custom communicators can also be created to group processes based on specific requirements or logical divisions.
-
The rank of a processor is a unique identifier assigned to each process within a communicator. It allows processes to distinguish and address each other in communication operations.
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 mpi_hello_world.jl
-
-
In this command, -np 4 specifies the desired number of processes.
-But it will probably not work since the version of mpiexec needs to match with the MPI version we are using from Julia. You can find the path to the mpiexec binary you need to use with these commands
-
julia>usingMPI
-julia>MPI.mpiexec_path
-
-
and then try again
-
$ /path/to/my/mpiexec -np 4 julia mpi_hello_world.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>usingMPI
-julia>mpiexec(cmd->run(`$cmd -np 4 julia mpi_hello_world.jl`))
-
-
Now, you should see output from 4 ranks.
-
-
-
-
-
-
-
-
-
-
-
-
Running MPI Programs in Jupyter Notebook with MPIClusterManagers¶
If you want to run your MPI code from a Jupyter Notebook, you can do so using the MPIClusterManagers package.
-
-
Load the packages and start an MPI cluster with the desired number of workers:
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
usingMPIClusterManagers
-# Distributed package is needed for addprocs()
-usingDistributed
-
-manager=MPIWorkerManager(4)
-addprocs(manager)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Run your MPI code inside a @mpi_do block to execute it on the cluster workers:
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
@mpi_domanagerbegin
-usingMPI
-comm=MPI.COMM_WORLD
-rank=MPI.Comm_rank(comm)
-println("Hello from process $rank")
-end
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
MPI is automatically initialized and finalized within the @mpi_do block.
MPI provides point-to-point communication using blocking send and receiving functions MPI.send, MPI.recv; or their non-blocking versions MPI.Isend, and MPI.Irecv!. These functions allow individual processes to send and receive data between each other.
To demonstrate asynchronous communication, let's modify the example using MPI.Isend and MPI.Irecv!:
-
usingMPI
-
-MPI.Init()
-
-comm=MPI.COMM_WORLD
-rank=MPI.Comm_rank(comm)
-
-# Asynchronous communication using MPI.Isend and MPI.Irecv!
-ifrank==0
-data="Hello from process $rank !"
-request=MPI.Isend(data,comm,dest=1)
-# Other computation can happen here
-MPI.Wait(request)
-elseifrank==1
-received_data=Array{UInt8}(undef,50)# Preallocate buffer
-request=MPI.Irecv!(received_data,comm,source=0)
-# Other computation can happen here
-MPI.Wait(request)
-println("Process $rank received: $(String(received_data))")
-end
-
-MPI.Finalize()
-
-
-
-
-
-
-
-
-
-
-
-
-
In this example, process 0 uses MPI.Isend to send the message asynchronously. This function returns immediately, allowing the sender process to continue its execution. However, the actual sending of data is done asynchronously in the background. Similar to MPI.Isend, MPI.Irecv! returns immediately, allowing the receiver process to continue executing.
-
-Important: In asynchronous communication, always use MPI.Wait() to ensure the communication is finished before accessing the send or receive buffer.
-
MPI provides collective communication functions for communication involving multiple processes. Let's explore some of these functions:
-
-
MPI.Gather: Gathers data from all processes to a single process.
-
MPI.Scatter: Distributes data from one process to all processes.
-
MPI.Bcast: Broadcasts data from one process to all processes.
-
MPI.Barrier: Synchronizes all processes.
-
-
Let's illustrate the use of MPI.Gather and MPI.Scatter with an example:
-
# TODO: check if this runs correctly
-usingMPI
-usingRandom
-
-MPI.Init()
-
-comm=MPI.COMM_WORLD
-rank=MPI.Comm_rank(comm)
-size=MPI.Comm_size(comm)
-
-# Root processor generates random data
-data=rand(rank==0?size*2:0)
-
-# Scatter data to all processes
-local_data=Vector{Float64}(undef,2)
-MPI.Scatter!(data,local_data,comm,root=0)
-
-# Compute local average
-local_average=sum(local_data)/length(local_data)
-
-# Gather local averages at the root processor
-gathered_averages=Vector{Float64}(undef,size)
-MPI.Gather!(local_average,gathered_averages,comm,root=0)
-
-ifrank==0
-# Compute global average of sub-averages
-global_average=sum(gathered_averages)/size
-println("Global average: $global_average")
-end
-
-MPI.Finalize()
-
-
-
-
-
-
-
-
-
-
-
In [ ]:
-
-
-
usingMPI
-usingRandom
-
-# TODO: check if this runs correctly
-
-MPI.Init()
-
-comm=MPI.COMM_WORLD
-rank=MPI.Comm_rank(comm)
-size=MPI.Comm_size(comm)
-
-# Root processor generates random data
-data=rand(rank==0?size*2:0)
-
-# Scatter data to all processes
-local_data=Vector{Float64}(undef,2)
-MPI.Scatter!(data,local_data,comm,root=0)
-
-# Compute local average
-local_average=sum(local_data)/length(local_data)
-
-# Gather local averages at the root processor
-gathered_averages=Vector{Float64}(undef,size)
-MPI.Gather!(local_average,gathered_averages,comm,root=0)
-
-ifrank==0
-# Compute global average of sub-averages
-global_average=sum(gathered_averages)/size
-println("Global average: $global_average")
-end
-
-MPI.Finalize()
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
In this example, the root processor generates random data and then scatters it to all processes using MPI.Scatter. Each process calculates the average of its local data, and then the local averages are gathered using MPI.Gather. The root processor computes the global average of all sub-averages and prints it.
\n",
- "Note: In our version of the jacobi method, we return after a given number of iterations. Other stopping criteria are possible. For instance, iterate until the difference between u and u_new is below a tolerance. \n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d6918c31",
- "metadata": {},
- "source": [
- "\n",
- "### Where can we exploit parallelism?\n",
- "\n",
- "Look at the two nested loops in the sequential implementation:\n",
- "\n",
- "```julia\n",
- "for t in 1:nsteps\n",
- " for i in 2:(n+1)\n",
- " u_new[i] = 0.5*(u[i-1]+u[i+1])\n",
- " end\n",
- " u, u_new = u_new, u\n",
- "end\n",
- "```\n",
- "\n",
- "- The outer loop cannot be parallelized. The value of `u` at step `t+1` depends on the value at the previous step `t`.\n",
- "- The inner loop can be parallelized.\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "798968b1",
- "metadata": {},
- "source": [
- "### The Gauss-Seidel method\n",
- "\n",
- "The usage of `u_new` seems a bit unnecessary at first sight, right?. If we remove it, we get another method called Gauss-Seidel.\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "0f77b547",
- "metadata": {},
- "outputs": [],
- "source": [
- "function gauss_seidel(n,niters)\n",
- " u = zeros(n+2)\n",
- " u[1] = -1\n",
- " u[end] = 1\n",
- " for t in 1:niters\n",
- " for i in 2:(n+1)\n",
- " u[i] = 0.5*(u[i-1]+u[i+1])\n",
- " end\n",
- " end\n",
- " u\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0dbc5358",
- "metadata": {},
- "source": [
- "Note that the final solution is the same (up to machine precision)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ca31518b",
- "metadata": {},
- "outputs": [],
- "source": [
- "gauss_seidel(5,1000)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c92e9c73",
- "metadata": {},
- "source": [
- "
\n",
- "Question: Which of the two loops in the Gauss-Seidel method are trivially parallelizable?\n",
- "
\n",
- "\n",
- "```julia\n",
- "for t in 1:niters\n",
- " for i in 2:(n+1)\n",
- " u[i] = 0.5*(u[i-1]+u[i+1])\n",
- " end\n",
- "end\n",
- "```\n",
- "\n",
- " a) Both of them\n",
- " b) The outer, but not the inner\n",
- " c) None of them\n",
- " d) The inner, but not the outer\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4edad93f",
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, c or d\n",
- "gauss_seidel_1_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9df06442",
- "metadata": {},
- "source": [
- "## Parallelization of the Jacobi method\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "97a2d7d5",
- "metadata": {},
- "source": [
- "### Parallelization strategy\n",
- "\n",
- "- Each worker updates a consecutive section of the array `u_new` \n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c4e2d000",
- "metadata": {},
- "source": [
- "### Data dependencies\n",
- "\n",
- "Recall:\n",
- "\n",
- "`u_new[i] = 0.5*(u[i-1]+u[i+1])`\n",
- "\n",
- "\n",
- "Thus, each process will need values from the neighboring processes to perform the update of its boundary values."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "428bce86",
- "metadata": {},
- "source": [
- "
\n",
- "\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1b3c8c05",
- "metadata": {},
- "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 Jacobi update locally in the processes."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "cd1e3760",
- "metadata": {},
- "source": [
- "
\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",
- "
\n",
- "\n",
- " a) Communication: O(P), computation: O(N/P)\n",
- " b) Communication: O(1), computation: O(N)\n",
- " c) Communication: O(P), computation: O(N)\n",
- " d) Communication: O(1), computation: O(N/P)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3a03fc4c",
- "metadata": {},
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, c or d\n",
- "jacobi_1_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8ed4129c",
- "metadata": {},
- "source": [
- "## Implementation\n",
- "\n",
- "We consider the implementation using MPI. The programming model of MPI is generally better suited for data-parallel algorithms like this one than the task-based model provided by Distributed.jl. In any case, one can also implement it using Distributed, but it requires some extra effort to setup remote channels right for the communication between neighbor processes.\n",
- "\n",
- "Take a look at the implementation below and try to understand it. Note that we have used MPIClustermanagers and Distributed just to run the MPI code on the notebook. When running it on a cluster, MPIClustermanagers and Distributed are not needed.\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e15082fb",
- "metadata": {},
- "outputs": [],
- "source": [
- "] add MPI MPIClusterManagers"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a66cbf9a",
- "metadata": {},
- "outputs": [],
- "source": [
- "using MPIClusterManagers \n",
- "using Distributed"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e0d63c6b",
- "metadata": {},
- "outputs": [],
- "source": [
- "if procs() == workers()\n",
- " nw = 3\n",
- " manager = MPIWorkerManager(nw)\n",
- " addprocs(manager)\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a0923606",
- "metadata": {},
- "outputs": [],
- "source": [
- "# Test cell, remove me\n",
- "u = [-1, 0, 0, 0, 0, 1]\n",
- "view(u, 6:6)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "68851107",
- "metadata": {
- "code_folding": []
- },
- "outputs": [],
- "source": [
- "@mpi_do manager begin\n",
- " using MPI\n",
- " comm = MPI.Comm_dup(MPI.COMM_WORLD)\n",
- " nw = MPI.Comm_size(comm)\n",
- " iw = MPI.Comm_rank(comm)+1\n",
- " function jacobi_mpi(n,niters)\n",
- " if mod(n,nw) != 0\n",
- " println(\"n must be a multiple of nw\")\n",
- " MPI.Abort(comm,1)\n",
- " end\n",
- " n_own = div(n,nw)\n",
- " u = zeros(n_own+2)\n",
- " u[1] = -1\n",
- " u[end] = 1\n",
- " u_new = copy(u)\n",
- " for t in 1:niters\n",
- " reqs = MPI.Request[]\n",
- " # Exchange cell values with neighbors\n",
- " if iw != 1\n",
- " neig_rank = (iw-1)-1\n",
- " req = MPI.Isend(view(u,2:2),comm,dest=neig_rank,tag=0)\n",
- " push!(reqs,req)\n",
- " req = MPI.Irecv!(view(u,1:1),comm,source=neig_rank,tag=0)\n",
- " push!(reqs,req)\n",
- " end\n",
- " if iw != nw\n",
- " neig_rank = (iw+1)-1\n",
- " s = n_own+1\n",
- " r = n_own+2\n",
- " req = MPI.Isend(view(u,s:s),comm,dest=neig_rank,tag=0)\n",
- " push!(reqs,req)\n",
- " req = MPI.Irecv!(view(u,r:r),comm,source=neig_rank,tag=0)\n",
- " push!(reqs,req)\n",
- " end\n",
- " MPI.Waitall(reqs)\n",
- " for i in 2:(n_own+1)\n",
- " u_new[i] = 0.5*(u[i-1]+u[i+1])\n",
- " end\n",
- " u, u_new = u_new, u\n",
- " end\n",
- " u\n",
- " @show u\n",
- " # Gather results in root process\n",
- " results = zeros(n+2)\n",
- " results[1] = -1\n",
- " results[n+2] = 1\n",
- " MPI.Gather!(view(u,2:n_own+1), view(results, 2:n+1), root=0, comm)\n",
- " if iw == 1\n",
- " @show results\n",
- " end \n",
- " end\n",
- " niters = 100\n",
- " load = 4\n",
- " n = load*nw\n",
- " jacobi_mpi(n,niters)\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "eff25246",
- "metadata": {},
- "source": [
- "
\n",
- "Question: How many messages per iteration are sent from a process away from the boundary?\n",
- "
\n",
- "\n",
- " a) 1\n",
- " b) 2\n",
- " c) 3\n",
- " d) 4\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "98bd9b5e",
- "metadata": {},
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, c or d\n",
- "jacobi_2_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "075dd6d8",
- "metadata": {},
- "source": [
- "
\n",
- "Question: After the end of the for-loop (line 43), ...\n",
- "
\n",
- "\n",
- " a) each worker holds the complete solution.\n",
- " b) the root process holds the solution. \n",
- " c) the ghost cells contain redundant values. \n",
- " d) all ghost cells contain the initial values -1 and 1. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c3b58002",
- "metadata": {},
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, c or d\n",
- "jacobi_3_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "4537661d",
- "metadata": {},
- "source": [
- "
\n",
- "Question: In line 35 of the code, we wait for all receive and send requests. Is it possible to instead wait for just the receive requests?\n",
- "
\n",
- "\n",
- " \n",
- " a) No, because the send buffer might be overwritten if we don't wait for send requests.\n",
- " b) No, because MPI does not allow an asynchronous send without a Wait().\n",
- " c) Yes, because each send has a matching receive, so all requests are done when the receive requests return. \n",
- " d) Yes, because there are no writes to the send buffer in this iteration."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e16ea5eb",
- "metadata": {},
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, c or d.\n",
- "jacobi_4_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c9aa2901",
- "metadata": {},
- "source": [
- "### Latency hiding\n",
- "\n",
- "Note that we only need communications to update the values at the boundary of the portion owned by each process. The other values (the one in green in the figure below) can be updated without communications. This provides the opportunity of overlapping the computation of the interior values (green cells in the figure) with the communication of the ghost values. This technique is called latency hiding, since we are hiding communication latency by overlapping it with communications that we need to do anyway.\n",
- "\n",
- "The modification of the implementation above to include latency hiding is leaved as an exercise (see below).\n"
- ]
- },
- {
- "attachments": {
- "fig.png": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAADyCAYAAABAvOgkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAewgAAHsIBbtB1PgAAABl0RVh0U29mdHdhcmUAd3d3Lmlua3NjYXBlLm9yZ5vuPBoAACAASURBVHic7d17fFxlgf/xz0nTO71ZKKWFUrByKbctlB+2ogJVi0q5rCCLrhVd7bqC6y6yK+yq4OruqoDu/lxXWryh+/uhFERIqSDIRS2lVKwCcoe2SIG29ELb9Jpk9o9n0pxMZ5I0JPM8ST7v12tePU/mzMy3STrNN+ec58nmzp1bQJIkSZKkyGpiB5AkSZIkCSyokiRJkqRE1JaMZ2ZZtiJGEEmSJElS31IoFCYCdzaPWxXULMtWzJkz5+lqh5IkSZIk9T3z5s2jUGiZFslTfCVJkiRJSbCgSpIkSZKSYEGVJEmSJCXBgipJkiRJSoIFVZIkSZKUBAuqJEmSJCkJFlRJkiRJUhIsqJIkSZKkJFhQJUmSJElJsKBKkiRJkpJgQZUkSZIkJcGCKkmSJElKggVVkiRJkpQEC6okSZIkKQkWVEmSJElSEiyokiRJkqQkWFAlSZIkSUmwoEqSJEmSkmBBlSRJkiQlwYIqSZIkSUqCBVWSJEmSlAQLqiRJkiQpCRZUSZIkSVISLKiSJEmSpCRYUCVJkiRJSbCgSpIkSZKSYEGVJEmSJCXBgipJkiRJSkJt7ACSJCk5Q4D3FrcXAvURs1TboOJtF33r7y1JSbCgSpKkUqOBG4vbbwSe7+Tz1ABvAqYCU4BhxY9fQrrl75+AzwP3ADMiZ5GkPseCKkmSusN/Ax8Ehpe573LSLagxvBc4E3gGuDpyFkmKyoIqSZJKbQQuK26v7+RzTKalnK4qPudRrzNXNfwJeBD4YxVf83hgDvArqldQDwJOyN0OAzLgJuCzVcogSXuwoEqSpFKbga++zue4DrgKeBh4Bfgw8IPX+ZzVcF3x1ptNBZZWuG+/agaRpFIWVEmS1B3+X+wAPch3CJNRbdnLxx0GHEe4RvjhTrzuBuC3xcd+CBjfieeQpC7lMjOSJKnUQUCheDs0cpZquxC4F7imiq/5MqEkPrWXj5tFmMxqzl4+7mnC5FdvAN5FuCZ43V4+hyR1C4+gSpIktTgUOAVoquJrTiZcn7sWuK8Kr7epeJOk5HgEVZIkKa73EY6EfjF2EEmKzYIqSZIkSUqCBVWSJEmSlASvQZUkSR01AphU4b5n8LpGSdLrZEGVJEkd9Tbgtgr3nQ7cWcUsfcWhVJ6l98TinycBX6mwzy3Akq4OJUndxYIqSZI6aiOV19v06Gn3mAB8tp19jiveylmOBVVSD2JBlSRJHfVrYGrsEH3MC8BXK9x3InAa8Afgjgr7/L47QklSd7GgSpIkpet54LIK932GUFCXtLGPJPUozuIrSZIkSUqCR1AlSVJ3OI4wcVKz43Pbnwa25sbfArZUI5QkKW0WVEmS1B3amln2CyXjH2JBlSRhQZUkSXvaAswrbnd2dt4ncs/Rnq3t71I1W4C1wIbYQbrZ94EDc+NDi3/OBO7Kffxe4N+qFUqSLKiSJKnUBuCvX+dz/Lp462m+Vrz1dtOAw8t8fFzx1mxtdeJIUmBBlSRJiutm4En2vgzWEZaheb4Tr/l3wLAO7PdCJ55bkjrNgipJkhTX48Xb3nq6eOuMSuumSlJULjMjSZIkSUqCBVWSJEmSlAQLqiRJkiQpCRZUSZIkSVISLKiSJEmSpCRYUCVJkiRJSXCZGSWvUChkdXV1D8XO0ZZrrrnm4OXLlw8B2Lx586sbNmxIbWHzfYH9itvbgBXxopQ1GJiYGz8FNMWJUt6ECRMOp/hLvQsvvPDFqVOnbo4cqaIsy26cNWvWVdV4rbq6uh8UCoWjqvFanTF//vwx999//2iAHTt21K9evTq1NR2HAhNy4ydiBWnDkbntF4D6WEEqmED4PAKsA9ZEzLKHkSNHjh4+fPgYgPHjx2+//PLLl8fO1JaBAweeMnPmzG7/Gt9+++1va2xsvKa7X6eztm/fXnPppZce3jxevXr1ih07dmyLmamMiYT/PyGsoftqvChl7Uf4+QNgK7AyYpY9DBo0aMiYMWMObh5/85vffCLLspiR2vOpM88888HYIarBgqrkffGLX8xOOOGEqbFztGX9+vW88MLun3uHAge3sXtsQ2n5DyNVx8cOUCr39aWxsfHwxP8TW1KtFyoUCkdlWZbsv89t27aV/tscEzFORyT7uSw6sv1doiot/NFt3LiRjRs3AjBgwIChWZaNjhypTTU1NVX52bCxsXFUyu8dWZa1et8Hkv1FXNFQWv+iNzVDaflFeRK2b99e+jWemvL/7VmWjYidoVo8xVeSJEmSlASPoKonui7Lsudih8jbtWvXFRRPsxk/fvzKVatWfTtypFJ/BbypuL0T+ELELOVMBc7Njb8MbImUpaza2tqvNDQ0ALB+/fp7syy7M3KkUrMLhcLkmAEKhcIdNTU198XMUGrz5s2fpHhEbfTo0dvWrVv3xciRSr0NeE9ufFmsIG34Ssvm+Adg/2XxopTz+49B08CwPeIVeOPNcfO0NmZM/aw1a56aANDYWGh4/vmRc2NnyhsyZOeYsWO3nhc5xq4syz4fOUMrO3bs2Af4XPN47NixN73yyiu/jRipnH8BBhS3nwW+EzFLOZ+g5ajuduDKaEnKGDdu3MkvvfTSGc3jhoaGfx44cGBjzEylCoXCvwPpHtbtJhZU9TiFQuHHs2bNuid2jrxJkyZdRrGgTpgwYdWqVau+GjlSqXfQUlB3Aanl+xitC+p/AasjZSmrX79+uwvq8uXLl5xxxhlJfQ4XLFgwDYhaULMs+1Vqn5fp06efTbGgjho1avu6deuSyke41jpfUL8GFCJlqSRXUCf/Dn7xo3hRyun/oVxBfRUe/mHcPK3tv//7TmgpqI0Nn/70yUnl++hHH598zjnPRy+oqb13XHDBBfuTK6gnnnjinXV1dakVwM/RUlBXkt7/7bNoXVCTynfCCSfszBfUH//4x1fPnz9/Z8xMperq6v6NPlhQPcVXkiRJkpQEC6okSZIkKQkWVEmSJElSEiyokiRJkqQkWFAlSZIkSUmwoEqSJEmSkmBBlSRJkiQlwYIqSZIkSUqCBVWSJEmSlAQLqiRJkiQpCRZUSZIkSVISLKiSJEmSpCRYUCVJkiRJSbCgSpIkSZKSYEGVJEmSJCXBgipJkiRJSoIFVZIkSZKUhNrYAdRtBgOTc+NlQFM7j6kBpuTGjwPbujiXJEmSJJVlQe29DgV+mxsPof2yOaDkMccAj3VxLkmSJEkqy1N8JUmSJElJsKBKkiRJkpJgQZUkSZIkJcGCKkmSJElKggVVkiRJkpQEC6okSZIkKQkWVEmSJElSEiyokiRJkqQkWFAlSZIkSUmwoEqSJEmSkmBB7b0KJeOsA4/ZpzuCSJIkSVJHWFB7ry0l446Uz4O7I4gkSZIkdYQFtfdaT+ujqJM68Jh3dFMWSZIkSWqXBbX32gIsz41Pb2f/YcBF3RdHkiRJktpmQe3d7s5tXwyMr7DfQOB64KBuTyRJkiRJFVhQe7d5QFNxexTwa+BcWq5H3Rc4H3gIOAf4fbUDSpIkSVKz2tgB1K0eBr4JfLo4PgSYX9xuovUvKF4ilNdnq5au8y6qq6s7K3aIvH/8x38c3Lz97LPPHgL8Z8Q45RyW2x5AevmOKhl/GdgaI0glDQ0Nu7cnTZr0jrq6uiER4+yhUCgcEzsDcEZdXd3Y2CHyvv3tb09s3l67du1g0vveP75k/B9RUnTY798Kk/aLnaK1pkEt2+sPgEmXxMuypxdfbJjYvN2vX03/6667N6l8/fs3vCF2BmBAXV1dUv82N23aNOSGG27YPV60aNFfACm8z+YNzG0fRnrvb4fktpN7/128ePGU/PgDH/jANbNnz26qtH8kHVmFo9exoPZ+nyn++SlaF9L89m+AvwDWVSvU65Fl2Z/HzlCqtrbln9LatWsPAP42Xpp29SftfAAfix2gVGNj4+7tESNGTAWmxkuTrOnFWzKGDGn5PcJrr702iPS/9xPPt/a4cEvVltGw5fzYKfI2bGjZ7tevpt/YsfVJ5UtELYl97w8YMKDVeP369TOAGXHSdMhBJPY5LDGQxPK9+uqrrcY1NTUXR4qiEhbU3q8R+Dvgu4QSejywH2GW3+cIR1TvIxxR7Qe8P/fYP1UzqCRJkqS+zYLadzxavLWlkZZTgJNxxRVXFBYsWHBb7BxtGTdu3BGNjY3DATZs2PDSyy+//GLsTCXGAQcWt7cAj0fMUs4+wOTc+Le0XD+dhCOOOOKEmpqafgD9+/d/BtjQzkOiybLskSq+1v2ESwSSNHz48AMnT548DmDbtm2bli9f/mTsTCWGA0fkxg/FCtKG/5PbfhLYFCtIBYcDI4rbL5PYL1f333//A0aPHn1Qcbse+GPkSG3avn37rmq8TpZlLwMp/99eM3ny5N1nyqxcufLx+vr60jXmY5tMy7wiL5Lee/GBhJ8/ILxvJPX+u88++wybMGHCkc3jQqGwlNZLNCalsbFxTewM1ZLNnTt39xciy7LD58yZ83TMQJIkSZKkvmHevHmHFQqFp5rHzuIrSZIkSUqCBVWSJCl9A9rfRZJ6PguqJElS+j5PmHV/39hBJKk7WVAlSZLS9xvgBOAR4LTIWSSp21hQJUmS0rcYeA04ALidsHxc/6iJJKkbWFAlSZLStwloKG4PAj4ALAMOjZZInTUodgApZRZUSZKknuGe3PYgwjqYi4G/jBNHnTCIcLr2XJz4SirLgipJktQz3AbU58YZMAb4r+J9Q2KE0l65lnAt8dnAfpGzSEmyoEqSJPUM9wNbynx8BDATeAI4tqqJtDc+CXwYaCQc9V4VN46UJguqJElSz7AW2FnhvgHABMJpwJdVLZE6ahrwjeL2Z4G7ImaRkmZBlSRJ6jkebuf+NwBfAm6pQhZ1zP7AjYRfItwCfD1uHClttbEDSJIkqcNuBd4NDCz5+FbgWeAm4HfA+irnUnm1wE+AA4GngAuBQsxAUuosqJIkST3HfcBGwlE5CKf8DiDMDnspnjqamq8Bbwc2A+cQlguS1AZP8ZUkSeo5VhAm2YEwo+/Pge8Sfqb7EXBAnFgq43zg7wlHTD9CmMRKUjssqJIkST3L44TS8wTwPuBTwKOEo6rX4893KZhK+MUBwL8DN0fMIvUovoFJkiT1LLcRjqKeW/xzG+FoXT3wTuAL8aIJmAjUAUOBhfj1kPaKBVWSJKlnWQj8X2Bl7mNPAJ8obn8BOLvaoQTAKMLXZyxhsqrzaTklW1IHWFAlSZJ6lueAz5T5+P8Q1trMCNejHl3NUGIAMB84ElgFnAVsiZpI6oEsqJIkSb3HPwB3APsAPwVGxo3TZ2TAd4AZhJl63wO8GDWR1ENZUCVJknqPRuAvgeXAm4D/jz/vVcOXgQ8Rlv05B3gkbhyp5/INS5IkqXdZR7gGtR54N/CtuHF6vb8C/okws/Ic4J64caSezYIqSZLU+zwCXAg0ESZPuiRqmt7rvcC1xe0rCMv8SHodLKiSJEm9002Ea1IBrgLOi5ilNzqN8DmuBb4HfCluHKl3sKBKkiT1Xl8H/pPwM98PgZPjxuk1pgG3AoOAnwF/HTeO1HtYUCVJknq3Swgz+jaXKZefeX1OBH5OmCl5IWGt04aoiaRexIIqSZLUuzURZvZ9ABgN3AUcHjVRz3US8AtgBGEypHMJM/dK6iIWVEmSpN5vG2FCn98BY4FfAm+MmqjnOZlQTkcCvwLOInxeJXUhC6okSVLfsBGYASwDxgP3AhNjBupB3kY4nXc4cD+h7G+JmkjqpSyokiRJfcdGwtqoTwIHAXdjSW3PmcAdwLDin+/Gcip1GwuqJElS37IaeAfwLOE0318DR0ZNlK6PEyaYGkyYtfdsPK1X6lYWVEmSpL5nFfBW4A/AgcAi4M1RE6Xns8A8oB/wA8KESDtiBpL6AguqJElS3/QKcBqwBBhFmADo1KiJ0jAQ+D7wleL4X4CP4FIyUlVYUCVJkvqu9YTTfe+m5RrLj0ZNFNf+hOVjLiQU0k8AV8QMJPU1FlRJkqS+bQtwBvATYADwXeAawqmtfckU4CFgOqG4vxuYGzWR1AdZUCVJkrQDuAC4DGgCLgFuJ6z52RfMBn4DTACeAd5COKosqcosqJIkSQIoAF8F3g/UAzOBpcDUmKG62UhgPnA9MAS4DTiRsAyPpAgsqJIkScq7mTDD7wpgEmGG388AWcRM3eEU4HeE2Xl3An9PWEbmtYiZpD7PgipJktTirwjrgn4jdpDIlhGuybyRcF3q1cBC4ICYobrIPsB/ESZDOgR4jnBK738QjiJLisiCKkmS1OJg4GTg2NhBErAROB/4OLAVOB14AvgkPfdnyNOBR4CLiuNrCUX8t9ESSWqlp765SJIkqTq+Q7gO9UFgBPAtwoRCx8QMtZfeRLi+9OeEo6bLCcvr/A2wOWIuSSUsqJIkSS2+AuwHnBM7SGKeIJwGexHhGs1phOs3rwMOjJirPfsRTk9+DJhFuNb0asIR8nsi5pJUgQVVkiSpxVbgVWBT7CAJagL+G5hMmPm2FvgYYVmWq4F940XbwxjgKsKR0s8QrqO9nXDU9x8Ia79KSpAFVZIkqcU0whqg58YOkrCXCEvRTAPuBQYRSuBKYB5xT/2dUsywHLgUGAosISyZcwbwdLxokjrCgipJktTi3cA1hGsT1bYHgdOAdxFK4BDChEqPEIrrbMI6o91tv+LrPkg47fjjxSwPAu8B3gz8ogo5JHWB2tgBJEmS1KPdVbxNB/4W+HPCGqOnEK75vAv4KXA/YUmX16sGOAo4lXCt8FuBfsX7dhDWcb2WsFyQpB7GgipJkqSu8EDxNo6wnuy5hMmI3lu8AbwCLCIc6XyueFsJrCNc45o3iDBr8KHAYYSZeKcQinD+yGwBWEpYs/V6YG3X/rUkVZMFVZIkSV3pJeBLxdvhhKJ6OmGpmrHA+4q3UrtombxoGG3/nLoZWExYNuanwAtdEVxSfBZUSZIkdZengH8t3gYSSup04EjgjYSjo+OBDOgPjCp5fBOhfD5LmC34ccIR2EeAxu6PL6naLKiSJEmqhh2Ecrmo5OP9gOHAYMJpvRlhuZ96XO5H6nMsqJIkSYqpEdhQvEnq41xmRpIkSZKUBAuqJEmSJCkJFlRJkiRJUhIsqJIkSZKkJFhQJUmSJElJsKBKkiRJkpLgMjOSJEkt5gK3A6/FDiJJfZEFVZIkqcWq4k2SFIGn+EqSJEmSkmBBlSRJkiQlwYIqSZIkSUqCBVWSJEmSlAQLqiRJkiQpCRZUSZIkSVISXGZGPcItt9wyMXaGttTV1e27YcOGQQAvvfTSpiVLlmyKnanEMGBEcXsnsCZilnIGAGNy41VAIVKWss4666xxNTU1NQDHH3/8uqOPPnpb7EyVDB48eNPpp5++vhqvdccddxywbdu2gdV4rc5YtGjR8Oeee244wJYtW7bfddddr8bOVGIgsF9u/GKsIG04MLe9FtgRK0gF+wKDitubSWz90mnTpg0bO3bsCIBhw4btPOecc1J7/23l7LPPXpllWbe//954442D+/fvv393v05n7dq1K7vhhhvGN48feuihNatWrdoZM1MZYwj/f0L4vt8cMUs5w4s3CO8bayNm2cMhhxwy8M/+7M92v//Onj07xfff3UaOHPnKqaeeuj12jmqwoCp5V155ZU1tbe3y2DnasnTpUh599NHYMdSNbr311t3bkydPprY23bfPXbt2fQu4uBqvtXPnzttqa2unVuO1OmPFihXccsstsWOoD1u8ePHu7UmTJnHeeedFTNO+u+++eyRVKPmDBg16V5ZlP+vu1+mshoYG3zt6ueXLl7N8ecuPlxdeeCHF30Mnqb6+/nTgztg5qiHdr4IkSZIkqU+xoEqSJEmSkpDuOWpSZX+ZZdni9nernu3bty+jeJ3FlClTHl62bNn7I0cqdT1wcnF7K3BMxCzlvB/499z4JCCpawUHDhz43I4d4dK7J5988tpp06ZdFTlSqWsLhcI7I2e4KsuyayNnaGXNmjU3AVMAJk6c+NqKFSuOjxyp1MeBy3aPvsEn40Wp4O/5793bb+cWzuauiGn29A9cTQNDwmDc0/BPl7X9gOo67LBbP//003dNAdi1a9e2LMuOjp0pr6mp6ZgETrVN7vPypz/9aV9gSfP4pJNOunzJkiU3RoxUzqPQ/L3PImB2xCzl/ARovgRkE8X34lRMnz79ow888MA/N483bdp0xKhRo3bFzFSqUCg8Qx88oGhBVY9TKBRenjVr1vOxc+RNmjSpqXl70KBBO4Ck8gH5i+oLpJevtIyuBFbHCNIR9fX1G88444ykPocLFizYGjsDsCG1z8v06dN3T+hTW1vbRHrf+60nszqNNdSkNUFYKyPZwjsS+7eZ0dQyqN0JF62KF2ZPAwfes/t7sFAoNKX2b+S2224bHTsDUEjt83LBBRfU58djxox5lfTeP3Lf+2wnvXz5CdWSe/8dPXr0uvz4tttuWz5//vykJsKqq6uLHSGKPtfIJUmSJElpsqBKkiRJkpJgQZUkSZIkJcGCKkmSJElKggVVkiRJkpQEC6okSZIkKQkWVEmSJElSEiyokiRJkqQkWFAlSZIkSUmojR2gl8uANwLjgUHAS8DjQOPrfN5JwFhgGLAReB5Y/TqfU5IkSZKiSv0I6kTgudxtWAceU1PymKO7K1zRJ3Ov9T+5DJ8CngSeAe4D7gAeIRTJLwFD9vJ1DgK+BbxYfM5fAwuBB4CXgd8BH6Htr+nluaz/2c7r/TOtP48faWf/m3L7fridfSVJkiRpD6kfQe0PHJobd7RQ5x8zsOvilDUq93orCCX6JuBdFfYfDXwOeAcwE9jUgde4CLiacBS2nAyYAnyPUCTPBtaX2e+xXNbzgL8DChWe8xxafx7PBL5fYd99gFnAgOL49xX2kyRJkqSKUi+oPU0N4Sjquwin8d4PPAxsBg4GzgL2Le77ZkLpnNPOc14JXJEb1wM/B5YBrwFjCEX3pOL9bwXuAt4CbC95rvuBBsLX/QBgMvDHMq/5BkLhzTsF6Ef505PfSks5XUs4UixJkiRJe8WC2rVOJnxO/wicz57l71LgZuC04vijwJeBFyo83yzgC7nxTwhHU9eV7HcF4YjoDwinDh8P/CvwmZL9NgFLgWnF8YwyGQFOpeVo9Q7CUeiRwAnAQ2X2n5HbvofKR2UlSZIkqaLUr0HtaWqBVYSCV674bQQ+AGwpjvsRimU5/QnXnGbF8Y3ABexZTpvNBz6WG3+ScHS11C9z2zPK3F/68Wv3cv97KuwjSZIkSW2yoHa9zxFOc61kNbAgN55aYb/zCBMjQSi0F9P+kckbCEdIIVyvekGZffIF8u2UP4reXDg3A18Fmko+nrcvcGxu/Msy+0iSJElSuyyoXWsX4Uhnex7ObU+ssM+Zue1babv05v00t/3WMvc/AGwrbo8ATiy5/yDgsOL2rwgzBP+hOH4LMLhk/9No+T5aQZjFV5IkSZL2mgW1az0FbO3Afmty2yMq7JMvl4v3IsOTue0jy9y/A/hNblx6VDQ/vrv4Z/NR0UHA9Db29+ipJEmSpE5zkqSuVW5pl3J25rYHlLl/IDAuN76E9tchhbDETX6t2DdU2O+XwDuL2zMIEzWRG+f3a/7z0tz9+SJ6Wpn9JUmSJGmvWVC7VkMXPc+okvGhZfdq39AKH89fhzqNMPNv85Hf5sL5CmHdVIBfE0r1AFoX2IOBScXtAnBvJ3NKkiRJkgU1Uf1KxvdQefbetmyr8PHfARsIRXggYXmcXxDWRW0+cptfLqYeeBB4G2GpmVHFx+fL6h8JpVaSJEmSOqU3FtRyp8z2NKVl9JvAz7rw+RuB+4BziuMZhILa1um6vyQU1H7AKcAteP2pJEmSpC6U+iRJ20vGgzrwmHJrf/Y022l9NPLYSju+DuXWQy03QVKl/TO8/lSSJElSF0q9oG4qGR/Qgcec1B1BIshfz3lmxb06L18opxCK/SnF8TPACyX7P0RYFxVCQT0KGFscNwD3d0NGSZIkSX1I6gX1NWB1bvyWDjzm492Updpuzm2fALy7i5//SWBVcbuGMEvvyOK43NHQXYR1UQGOAGbn7vste/4yQZIkSZL2SuoFFVqvAfoJ2r5u9iO0LJ/S091CmHio2feAQ/bi8fvRUjgryc/me3Fuu/T03mb54npxhY9LkiRJUqf0hIL6o9z20cD1wPCSfYYAXwCuA7ZUKVd3ayIU7ubrcMcCS4GPUnkiqH7AqcA8YCXwxnZeI18sB+det9JyMeX2L/24JEmSJHVKT5jF92fAIlpO7/0AcEbxY68RjhS+mbDmZxPh1NOfVj9mt1hK+Pv8kDBB1Gjgu8DXCUeWXyQU2JGEo6vHAfvsxfOXK5bLgPUV9n8UWEPriai20footyRJkiR1Sk8oqE3A+cBdwJHFjw1nz2sytwJ/A9xavWhVMR9YAfyAsE4pwAjg9HYe9zKwsZ19XgSeAg7PfazS6b0Q1kW9B/iL3McWsedsy5IkSZK013pCQYUwmc+JwCXAB2ldqDYQStx/AE8Qlj+ZX3J/d3oi93qPdfAxL+Qes7qtHYuWAscAfw6cSziNt3Q5nS2EiY9+BdxJODra2IHn/ibw9tz45ko7Fv2IcCpxR/eXJEmSpA7pKQUVoB74UvE2jHBq70b2PB21ALy/irl+yt6fUvxA8bY3moCbijcIpzSPJnwNy30eOupbxVtHLSzeJEmSJKlL9aSCmreZljU5+6r64k2SJEmSeoWeMIuvJEmSJKkPsKBKkiRJkpJgQZUkSZIkJaGnXoPaWR8CBnfB8zwN3NcFz6NOyLLsxrq6up2xc+Rdfvnlw5u3ly5dOhV4KWKcct6Q2x5CevlK/13+gTAxWDJ27mz5lpsyZcrFF1100YcjxtlDoVAYFTsDcHldXd2nYofI+/73vz+6eXv58uUjSO97f2ir0UlcFylHxyzkXO7kjNgxWmnIfw5XHQEDbo8XZk+PP56NaN4eMGDAEI2W2AAAAwZJREFUkLq6utS+B/vHDgAMTu3zsnXr1pobbrhh93jhwoVXAf8SL1FZ+fePk0nv/W10bju599+FCxe2ev/94Ac/uGL27Nmx4lTSJw8m9rWCehWwfxc8z/VYUGMa3f4u1ZVl2e7thoaGAcAB8dK0KyPtfNA1/067VKFQ2L1dU1OzD7BPvDTJGla8JaOmpuX/9sbGxhpS/97fTgq/aKhsF4PZ1SW/6O0mjbXQuG/sFHmNrRd86wnvvzEk93nJv3cANDY2jgRGxknTIQNJ7HNYIrmvcWPpP84sSypfX9YnW7kkSZIkKT197Qjq4XRNKU/q9NLe7sorr2xasGDBh2LnaMuxxx571Lhx40YCrFq16k+PPfbYC7EzlTgQOLi4vRl4JGKWcoYBx+bGDwKNFfaNYsaMGW+ura3tBzB06NAnsizr7NrD3a6pqempar1Wv379Pl8oFJI6YpU3bty4g2fOnHkgwKZNmzYuXrz4j7EzlRgBHJ0bL4oVpA1vyW0/BrwWK0gFR9FyZOtFYGXELHs48sgjx0+YMGEiwMiRIzdnWZba+28ro0aN2lqN16mtrX24qakp2f/bsyzrN3PmzDc3j5ctW/bImjVrUlvi8Big+RKjlYTv/5RMAA4qbr9GeP9IxtixY0ccd9xxu99/syx7IMuyQluPiam2tjbp946ulM2dO3f3FyLLssPnzJnzdMxAkiRJkqS+Yd68eYcVCoXdv1z3FF9JkiRJUhIsqJIkSZKkJFhQJUmSJElJsKBKkiRJkpJgQZUkSZIkJcGCKkmSJElKggVVkiRJkpQEC6okSZIkKQkWVEmSJElSEiyokiRJkqQkWFAlSZIkSUmwoEqSJEmSkmBBlSRJkiQlwYIqSZIkSUqCBVWSJEmSlAQLqiRJkiQpCRZUSZIkSVISLKiSJEmSpCRYUCVJkiRJSbCgSpIkSZKSYEGVJEmSJCWhNj8oFAoT582bFyuLJEmSJKkPKRQKE/Pj2pL77ywUCtVLI0mSJElSkaf4SpIkSZKSYEGVJEmSJCXhfwG8dZWmdb83nAAAAABJRU5ErkJggg=="
- }
- },
- "cell_type": "markdown",
- "id": "5ae8701f",
- "metadata": {},
- "source": [
- "
\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "850b1848",
- "metadata": {},
- "source": [
- "### Complexity analysis\n",
- "\n",
- "The complexity analysis is let as an exercise (see below).\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "47643bf6",
- "metadata": {},
- "source": [
- "## Exercises"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0edeee84",
- "metadata": {},
- "source": [
- "### Exercise 1\n",
- "\n",
- "Transform the following parallel implementation of the 1d Jacobi method (it is copied from above) to use latency hiding (overlap between computation of interior values and communication)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "66db180d",
- "metadata": {},
- "outputs": [],
- "source": [
- "@mpi_do manager begin\n",
- " using MPI\n",
- " comm = MPI.Comm_dup(MPI.COMM_WORLD)\n",
- " nw = MPI.Comm_size(comm)\n",
- " iw = MPI.Comm_rank(comm)+1\n",
- " function jacobi_mpi(n,niters)\n",
- " if mod(n,nw) != 0\n",
- " println(\"n must be a multiple of nw\")\n",
- " MPI.Abort(comm,1)\n",
- " end\n",
- " n_own = div(n,nw)\n",
- " u = zeros(n_own+2)\n",
- " u[1] = -1\n",
- " u[end] = 1\n",
- " u_new = copy(u)\n",
- " for t in 1:niters\n",
- " reqs = MPI.Request[]\n",
- " # Exchange cell values with neighbors\n",
- " if iw != 1\n",
- " neig_rank = (iw-1)-1\n",
- " req = MPI.Isend(view(u,2:2),comm,dest=neig_rank,tag=0)\n",
- " push!(reqs,req)\n",
- " req = MPI.Irecv!(view(u,1:1),comm,source=neig_rank,tag=0)\n",
- " push!(reqs,req)\n",
- " end\n",
- " if iw != nw\n",
- " neig_rank = (iw+1)-1\n",
- " s = n_own+1\n",
- " r = n_own+2\n",
- " req = MPI.Isend(view(u,s:s),comm,dest=neig_rank,tag=0)\n",
- " push!(reqs,req)\n",
- " req = MPI.Irecv!(view(u,r:r),comm,source=neig_rank,tag=0)\n",
- " push!(reqs,req)\n",
- " end\n",
- " MPI.Waitall(reqs)\n",
- " for i in 2:(n_own+1)\n",
- " u_new[i] = 0.5*(u[i-1]+u[i+1])\n",
- " end\n",
- " u, u_new = u_new, u\n",
- " end\n",
- " u\n",
- " @show u\n",
- " # Gather results in root process\n",
- " results = zeros(n+2)\n",
- " results[1] = -1\n",
- " results[n+2] = 1\n",
- " MPI.Gather!(view(u,2:n_own+1), view(results, 2:n+1), root=0, comm)\n",
- " if iw == 1\n",
- " @show results\n",
- " end \n",
- " end\n",
- " niters = 100\n",
- " load = 4\n",
- " n = load*nw\n",
- " jacobi_mpi(n,niters)\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ebb650d0",
- "metadata": {},
- "source": [
- "### Exercise 2\n",
- "\n",
- "Compute the complexity of the communication and computation of the three data partition strategies (1d block partition, 2d block partition, and 2d cyclic partition) when computing a single iteration of the Jacobi method in 2D. Assume that the grid is of size $N \\times N$ and the number of processes $P$ is a perfect square number, i.e. $\\sqrt{P}$ is an integer. Hint: For the complexity analysis, you can ignore the effect of the boundary conditions.\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7b3d7cb3",
- "metadata": {},
- "outputs": [],
- "source": [
- "# TODO"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6d3430ad",
- "metadata": {},
- "source": [
- "# License\n",
- "\n",
- "TODO: replace link to website\n",
- "\n",
- "This notebook is part of the course [Programming Large Scale Parallel Systems](http://localhost:8000/) at Vrije Universiteit Amsterdam and may be used under a [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) license."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3d72ff47",
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Julia 1.9.1",
- "language": "julia",
- "name": "julia-1.9"
- },
- "language_info": {
- "file_extension": ".jl",
- "mimetype": "application/julia",
- "name": "julia",
- "version": "1.9.1"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/docs/src/notebooks/jacobi_method.md b/docs/src/notebooks/jacobi_method.md
deleted file mode 100644
index d1d92fe..0000000
--- a/docs/src/notebooks/jacobi_method.md
+++ /dev/null
@@ -1,30 +0,0 @@
-```@meta
-EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/jacobi_method.ipynb"
-```
-
-```@raw html
-
- Tip
-
-
-
- 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
- .
-
-
-
-
-```
-
-```@raw html
-
-
-```
diff --git a/docs/src/notebooks/julia_async.ipynb b/docs/src/notebooks/julia_async.ipynb
deleted file mode 100644
index 9fab7f6..0000000
--- a/docs/src/notebooks/julia_async.ipynb
+++ /dev/null
@@ -1,794 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "id": "f64b009a",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "### Programming large-scale parallel systems\n",
- "\n",
- "\n",
- "# Asynchronous programming in Julia\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "bf68ad38",
- "metadata": {},
- "source": [
- "## Contents\n",
- "\n",
- "In this notebook, we will learn the basics of asynchronous programming in Julia. In particular, we will learn about:\n",
- "\n",
- "- Tasks\n",
- "- Channels\n",
- "\n",
- "Understanding these concepts is important to learn later distributed computing."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "caf64254",
- "metadata": {},
- "source": [
- "## Tasks\n",
- "\n",
- "### Creating a task\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "fe668cb1",
- "metadata": {},
- "outputs": [],
- "source": [
- "function work()\n",
- " println(\"Starting work\")\n",
- " sleep(7)\n",
- " a = rand(3,3)\n",
- " b = rand(3,3)\n",
- " r = a + b\n",
- " println(\"Finishing work\")\n",
- " r\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "67ee0328",
- "metadata": {},
- "outputs": [],
- "source": [
- "t = Task(work)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e459c5c2",
- "metadata": {},
- "source": [
- "### Scheduling a task\n",
- "\n",
- "The task has been created, but the corresponding work has not started. Note that we do not see any output from function `work` yet. To run the task we need to schedule it."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8778c199",
- "metadata": {},
- "outputs": [],
- "source": [
- "schedule(t)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f1fb9283",
- "metadata": {},
- "source": [
- "### Fetching the task result\n",
- "\n",
- "The task has been executed, but we do not see the result. To get the result we need to fetch it."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "0c7b626e",
- "metadata": {},
- "outputs": [],
- "source": [
- "fetch(t)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "fedbbd71",
- "metadata": {},
- "source": [
- "### Tasks run asynchronously\n",
- "\n",
- "It is important to note that tasks run asynchronously. To illustrate this let's create and schedule a new task."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4ccc996c",
- "metadata": {},
- "outputs": [],
- "source": [
- "t = Task(work)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "015bea27",
- "metadata": {},
- "outputs": [],
- "source": [
- "schedule(t)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "5ec0718e",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a70fcbe8",
- "metadata": {},
- "outputs": [],
- "source": [
- "sin(4π)*exp(-0.1)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "6def444b",
- "metadata": {},
- "outputs": [],
- "source": [
- "1 + 1"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d483d4d0",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "910323fd",
- "metadata": {},
- "source": [
- "### Tasks do not run in parallel\n",
- "\n",
- "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](https://en.wikipedia.org/wiki/Leibniz_formula_for_pi). The quality of the approximation increases with the value of `n`."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "b53ac640",
- "metadata": {},
- "outputs": [],
- "source": [
- "function compute_π(n)\n",
- " s = 1.0\n",
- " for i in 1:n\n",
- " s += (isodd(i) ? -1 : 1) / (i*2+1)\n",
- " end\n",
- " 4*s\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7614df94",
- "metadata": {},
- "source": [
- " Call this function with a large number. Note that it will take some time."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "63f2aec3",
- "metadata": {},
- "outputs": [],
- "source": [
- "compute_π(4_000_000_000)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "840ff590",
- "metadata": {},
- "source": [
- "Create a task that performs this computation."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f16ade6c",
- "metadata": {},
- "outputs": [],
- "source": [
- "fun = () -> compute_π(4_000_000_000)\n",
- "t = Task(fun)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "60c2567c",
- "metadata": {},
- "source": [
- "Schedule the tasks and then try to execute the 2nd cell bellow. Note that the current process will be busy running the task."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "1f28a388",
- "metadata": {},
- "outputs": [],
- "source": [
- "schedule(t)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "52e060e7",
- "metadata": {},
- "outputs": [],
- "source": [
- "1+1"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "25048665",
- "metadata": {},
- "source": [
- "### Tasks take a function with no arguments\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "87397749",
- "metadata": {},
- "outputs": [],
- "source": [
- "add(a,b) = a + b"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e99766c7",
- "metadata": {},
- "outputs": [],
- "source": [
- "t = Task(add)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "1a785bae",
- "metadata": {},
- "outputs": [],
- "source": [
- "schedule(t)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ebfce725",
- "metadata": {},
- "source": [
- "If we need, we can capture variables in the function to be run by the task as shown in the next cells."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "06f54fa0",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = rand(3,3)\n",
- "b = rand(3,3);"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4c9e586e",
- "metadata": {},
- "outputs": [],
- "source": [
- "fun = () -> a + b\n",
- "t = Task(fun)\n",
- "schedule(t)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "cd829a64",
- "metadata": {},
- "source": [
- "### Useful macro: `@async`\n",
- "\n",
- "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 `@async` macro. 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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9b5152d1",
- "metadata": {},
- "outputs": [],
- "source": [
- "@async a + b"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ac3262d8",
- "metadata": {},
- "source": [
- "### Another useful macro: `@sync`\n",
- "\n",
- "This macro is used to wait for all the tasks created with `@async` in a given block of code. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ef2e49f9",
- "metadata": {},
- "outputs": [],
- "source": [
- "@sync begin\n",
- " @async sleep(3)\n",
- " @async sleep(4)\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "98f0d685",
- "metadata": {},
- "source": [
- "## Channels\n",
- "\n",
- "### Sending data between tasks\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "b88d5308",
- "metadata": {},
- "outputs": [],
- "source": [
- "chnl = Channel{Int}()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3421c7a6",
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "@async begin\n",
- " for i in 1:5\n",
- " put!(chnl,i)\n",
- " end\n",
- " close(chnl)\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1df0508c",
- "metadata": {},
- "source": [
- "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. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "64b9436e",
- "metadata": {},
- "outputs": [],
- "source": [
- "take!(chnl)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "757a1b07",
- "metadata": {},
- "source": [
- "### Channels are iterable\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "2fc22dfa",
- "metadata": {},
- "outputs": [],
- "source": [
- "chnl = Channel{Int}()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "fa62a4df",
- "metadata": {},
- "outputs": [],
- "source": [
- "@async begin\n",
- " for i in 1:5\n",
- " put!(chnl,i)\n",
- " end\n",
- " close(chnl)\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e511e19b",
- "metadata": {},
- "outputs": [],
- "source": [
- "for i in chnl\n",
- " @show i\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b1a2a557",
- "metadata": {},
- "source": [
- "### Calls to `put!` and `take!` are blocking\n",
- "\n",
- "Note that `put!` and `take!` are blocking operations. Calling `put!` blocks the tasks until another task calls `take!` and viceversa. Thus, we need at least 2 tasks for this to work. If we call `put!` and `take!` 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 how `put!` blocks until we call `take!`. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f34373ca",
- "metadata": {},
- "outputs": [],
- "source": [
- "chnl = Channel{Int}()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "bfde2ecd",
- "metadata": {},
- "outputs": [],
- "source": [
- "@async begin\n",
- " for i in 1:5\n",
- " put!(chnl,i)\n",
- " println(\"I have put $i\")\n",
- " end\n",
- " close(chnl)\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "5830659d",
- "metadata": {},
- "outputs": [],
- "source": [
- "take!(chnl)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "006140bd",
- "metadata": {},
- "source": [
- "### Buffered channels\n",
- "\n",
- "We can be a bit more flexible and use a buffered channel. In this case, `put!` will block only if the channel is full and `take!` will block if the channel is empty. We repeat previous example, but with a buffered channel of size 2. Note that we can call `put!` until the channel is full. At this point, we need to wait to until we call `take!` which removes an item from the channel, making room for a new item."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "dfe06b5f",
- "metadata": {},
- "outputs": [],
- "source": [
- "buffer_size = 2\n",
- "chnl = Channel{Int}(buffer_size)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "6289bc2e",
- "metadata": {},
- "outputs": [],
- "source": [
- "@async begin\n",
- " for i in 1:5\n",
- " put!(chnl,i)\n",
- " println(\"I have put $i\")\n",
- " end\n",
- " close(chnl)\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "2a87cd5f",
- "metadata": {},
- "outputs": [],
- "source": [
- "take!(chnl)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9ddd66ca",
- "metadata": {},
- "source": [
- "## Questions"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "dd45ae08",
- "metadata": {},
- "outputs": [],
- "source": [
- "t = @elapsed compute_π(100_000_000)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "2a23b1c3",
- "metadata": {},
- "source": [
- "
\n",
- "Question (Q1): How long will the compute time of next cell be? \n",
- "
- 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
- .
-
-
-
-
-```
-
-```@raw html
-
-
-```
diff --git a/docs/src/notebooks/julia_basics.ipynb b/docs/src/notebooks/julia_basics.ipynb
deleted file mode 100644
index 993708d..0000000
--- a/docs/src/notebooks/julia_basics.ipynb
+++ /dev/null
@@ -1,1625 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "id": "faecad82",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "### Programming large-scale parallel systems\n",
- "\n",
- "\n",
- "# Julia basics"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0f54a3f8",
- "metadata": {},
- "source": [
- "## Contents\n",
- "\n",
- "In this notebook, we will cover the basic parts of Julia needed later to learn parallel computing. In particular, we will learn about:\n",
- "\n",
- "- Variables\n",
- "- Functions\n",
- "- Arrays\n",
- "\n",
- "For a more general introduction to Julia see the nice tutorials made available by JuliaAcademy [here](https://github.com/JuliaAcademy/JuliaTutorials/) or the official Julia educational resources [here](https://julialang.org/learning/).\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a23a48f4",
- "metadata": {},
- "source": [
- "## Using Jupyter notebooks in Julia\n",
- "\n",
- "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.\n",
- "\n",
- "
\n",
- "Tip: Did you know that Jupyter stands for Julia, Python and R?\n",
- "
\n",
- "\n",
- "### How to start a Jupyter nootebook in Julia\n",
- "\n",
- "To run a Julia Jupyther notebook, open a Julia REPL and type\n",
- "\n",
- "```julia\n",
- "julia> ]\n",
- "pkg> add IJulia\n",
- "julia> using IJulia\n",
- "julia> notebook()\n",
- "```\n",
- "A new browser window will open. Navigate to the corresponding notebook and open it.\n",
- "\n",
- "
\n",
- "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.\n",
- "
\n",
- "\n",
- "\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "bcefe994",
- "metadata": {},
- "source": [
- "### Running a cell\n",
- "To run a cell, click on a cell and press `Shift` + `Enter`. You can also use the \"Run\" button in the toolbar above."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4eae9070",
- "metadata": {},
- "outputs": [],
- "source": [
- "1+3\n",
- "4*5"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0a02fb75",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3923b7a5",
- "metadata": {},
- "outputs": [],
- "source": [
- "1+3\n",
- "4*5;"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ffd3a444",
- "metadata": {},
- "source": [
- "### Cell order is important\n",
- "\n",
- "Running the two cells below in reverse order won't work (try it). "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f7597471",
- "metadata": {},
- "outputs": [],
- "source": [
- "foo() = \"Well done!\""
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4a5d4f50",
- "metadata": {},
- "outputs": [],
- "source": [
- "foo()"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "92112bd1",
- "metadata": {},
- "source": [
- "### REPL modes\n",
- "\n",
- "This is particular to Julia notebooks. You can use package, help, and shell mode just like in the Julia REPL."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9a445020",
- "metadata": {},
- "outputs": [],
- "source": [
- "] add MPI"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "edefc54f",
- "metadata": {},
- "outputs": [],
- "source": [
- "? print"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e69fd5f1",
- "metadata": {},
- "outputs": [],
- "source": [
- "; ls"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8b94fc33",
- "metadata": {},
- "source": [
- "## Variables\n",
- "\n",
- "The usage of variables in Julia is pretty similar to Python and quite different from C/C++ and Fortran. However, there are also some differences with Python. \n",
- "\n",
- "### Creating a variable\n",
- "\n",
- "A variable is a name associated (bound) to a value. We associate variables with values with `=` as usual."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "149b7a62",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = 1"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "77b557a8",
- "metadata": {},
- "source": [
- "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)."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ba08ea11",
- "metadata": {},
- "source": [
- "### Re-assign a variable\n",
- "\n",
- "We can re-assign a variable, even with a value of another type. However, avoid changing the variable type for performance reasons."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f62e49c9",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = 2"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "86ac71ec",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = 1.0"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e67d8480",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = \"Hi!\""
- ]
- },
- {
- "cell_type": "markdown",
- "id": "767563a7",
- "metadata": {},
- "source": [
- "### Unreachable objects\n",
- "\n",
- "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"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "202e222b",
- "metadata": {},
- "outputs": [],
- "source": [
- " a = zeros(300000000);"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "35742a64",
- "metadata": {},
- "source": [
- "If we re-assign the variable to another value, the large array will be inaccessible."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e6cd1c24",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = nothing"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "2883bf04",
- "metadata": {},
- "source": [
- "### Garbage collector\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "89516ede",
- "metadata": {},
- "outputs": [],
- "source": [
- "GC.gc()"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c110aa22",
- "metadata": {},
- "source": [
- "### Type declarations are optional\n",
- "\n",
- "Julia knows the type of the object associated with a variable."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c8d93d64",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = 1\n",
- "typeof(a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "df2e48c2",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "95fd3ef4",
- "metadata": {},
- "outputs": [],
- "source": [
- "c::Int = 1\n",
- "typeof(c)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "67883016",
- "metadata": {},
- "source": [
- "If you annotate a variable with a type, then it cannot refer to objects of other types."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f3cdc27b",
- "metadata": {},
- "outputs": [],
- "source": [
- "c = \"I am a string\""
- ]
- },
- {
- "cell_type": "markdown",
- "id": "64f5aa94",
- "metadata": {},
- "source": [
- "### Variable names\n",
- "\n",
- "There is a great flexibility to choose variable names in Julia. See all the rules in the [manual](https://docs.julialang.org/en/v1/manual/variables/#man-allowed-variable-names) if you are interested."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "09eb4365",
- "metadata": {},
- "source": [
- "We can use Unicode (UTF-8 encoding) characters in variables and function names."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "94667496",
- "metadata": {},
- "outputs": [],
- "source": [
- "🐱 = \"I am a cat\"\n",
- "🐶 = \"I am a dog\"\n",
- "🐱 == 🐶"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "beaf9680",
- "metadata": {},
- "source": [
- "We can also use Greek letters and other mathematical symbols. Just write the corresponding [LaTeX](https://oeis.org/wiki/List_of_LaTeX_mathematical_symbols) command and press `Tab`. For example: `ω` is written `\\omega` + `Tab`"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "b9b3079e",
- "metadata": {},
- "outputs": [],
- "source": [
- "ω = 1.234\n",
- "sin(ω)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "4f7de553",
- "metadata": {},
- "source": [
- "In fact, some useful mathematical constants are predefined in Julia with math Greek letters."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "89d6d94e",
- "metadata": {},
- "outputs": [],
- "source": [
- "sin(π/2)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d36ba0ae",
- "metadata": {},
- "source": [
- "\n",
- "
\n",
- "Question: Which will be the value of `x` in the last line ? (Think your answer before executing next cell to find out the result) \n",
- "
\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f6f4e13e",
- "metadata": {},
- "outputs": [],
- "source": [
- "x = 1\n",
- "y = x\n",
- "y = 2\n",
- "x"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "4d2cb752",
- "metadata": {},
- "source": [
- "## Functions\n",
- "\n",
- "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.\n",
- "\n",
- "### Defining functions\n",
- "\n",
- "Functions are defined as shown in next cell. The closing `end` is necessary. Do to forget it! However, the `return` is optional. The value of last line is returned by default. Indentation is recommended, but it is also optional. That's why the closing `end` is needed."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "266b0e1b",
- "metadata": {},
- "outputs": [],
- "source": [
- "function add(a,b)\n",
- " return a + b \n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "de95c63c",
- "metadata": {},
- "source": [
- "Once defined, a function can be called using bracket notation as you would expect."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "39b452a7",
- "metadata": {},
- "outputs": [],
- "source": [
- "add(1,3)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "afceb0fc",
- "metadata": {},
- "source": [
- "### Broadcast syntax\n",
- "\n",
- "We can apply functions to arrays element by element using broadcast (dot) syntax."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e65be5e8",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = [1,2,3]\n",
- "b = [4,5,6]\n",
- "add.(a,b)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e4ca76a4",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ae0b3f6d",
- "metadata": {},
- "outputs": [],
- "source": [
- "a .* b"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b131e208",
- "metadata": {},
- "source": [
- "
\n",
- "Question: Which will be the value of `x` in the last line ?\n",
- "
\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3d637d87",
- "metadata": {},
- "outputs": [],
- "source": [
- "function q(x)\n",
- " x = 2\n",
- " x\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "882d8e59",
- "metadata": {},
- "outputs": [],
- "source": [
- "x = 1\n",
- "y = q(x)\n",
- "x"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d20bc3a4",
- "metadata": {},
- "source": [
- "### Defining functions (shorter way)\n",
- "\n",
- "For short functions, we can skip the `function` and `end` keywords as follows.\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7815aea5",
- "metadata": {},
- "outputs": [],
- "source": [
- "add_short(a,b) = a+b"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7aa48db2",
- "metadata": {},
- "outputs": [],
- "source": [
- "add_short(1,3)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f39c31c0",
- "metadata": {},
- "source": [
- "### Anonymous (lambda) functions\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "94ce2be7",
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "add_anonymous = (a,b) -> a+b"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b15b84d8",
- "metadata": {},
- "source": [
- "We can call the function by using the variable name."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "96d80a59",
- "metadata": {},
- "outputs": [],
- "source": [
- "add_anonymous(2.0,3.5)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1723e4f1",
- "metadata": {},
- "source": [
- "Note that `add_anonymous` is 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)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "5bf111a0",
- "metadata": {},
- "outputs": [],
- "source": [
- "nameof(add_anonymous)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7e23f9f8",
- "metadata": {},
- "outputs": [],
- "source": [
- "nameof(add)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e530acd3",
- "metadata": {},
- "source": [
- "### Functions are first-class objects\n",
- "\n",
- "We can work with Julia functions like with any other type of object. For instance, we can assign functions to variables."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e57fe028",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = add"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7a8f4bb0",
- "metadata": {},
- "source": [
- "Now, we can call the function using the variable name."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c7e4d4ec",
- "metadata": {},
- "outputs": [],
- "source": [
- "a(4,5)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c01be92e",
- "metadata": {},
- "source": [
- "We can also create an array of functions (this will not work in Python)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3837bf25",
- "metadata": {},
- "outputs": [],
- "source": [
- "funs = [+,-,*]"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0aac151b",
- "metadata": {},
- "source": [
- "To call a specific function in the array, we index the array and then call the returned function"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "aab45992",
- "metadata": {},
- "outputs": [],
- "source": [
- "funs[2](2,3)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f0c668c8",
- "metadata": {},
- "source": [
- "### Higher-order functions\n",
- "\n",
- "Higher order functions are functions that take and/or return other functions. And example is the `count` function in Julia.\n",
- "\n",
- "For instance, we can pass a user-defined function to count the number of even elements in an array."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "22115f82",
- "metadata": {},
- "outputs": [],
- "source": [
- "func = i->i%2==0\n",
- "a = [1,2,3,5,32,2,4]\n",
- "count(func,a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0b352cdc",
- "metadata": {},
- "source": [
- "### Do-blocks\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "533a21f3",
- "metadata": {},
- "outputs": [],
- "source": [
- "function f(i)\n",
- " m = i%2\n",
- " m != 0\n",
- "end\n",
- "count(f,[1,2,3,5,32,2,4])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "68b5165a",
- "metadata": {},
- "outputs": [],
- "source": [
- "count([1,2,3,5,32,2,4]) do i\n",
- " m = i%2\n",
- " m != 0\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "559a084d",
- "metadata": {},
- "source": [
- "### Returning multiple values\n",
- "\n",
- "Julia functions always return a single variable. To return multiple values, we can wrap them in a tuple.\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "eba30561",
- "metadata": {},
- "outputs": [],
- "source": [
- "function divrem(a,b)\n",
- " α = div(a,b)\n",
- " β = rem(a,b)\n",
- " (α,β)\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1a6e9d0d",
- "metadata": {},
- "source": [
- "The output is a tuple as expected, but we can recover the individual values by unpacking the tuple."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "cb0c2806",
- "metadata": {},
- "outputs": [],
- "source": [
- "d,r = divrem(10,3)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "31428737",
- "metadata": {},
- "outputs": [],
- "source": [
- "d"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "76424bb5",
- "metadata": {},
- "outputs": [],
- "source": [
- "r"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c73ae75b",
- "metadata": {},
- "source": [
- "### Variable number of input arguments\n",
- "\n",
- "Functions with multiple arguments are also supported. The following example iterates over the given arguments and prints them. `args` is just a tuple with all arguments."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a7f94ff0",
- "metadata": {},
- "outputs": [],
- "source": [
- "function showargs(args...)\n",
- " for (i,arg) in enumerate(args)\n",
- " println(\"args[$i] = $arg\")\n",
- " end\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7828534e",
- "metadata": {},
- "outputs": [],
- "source": [
- "showargs(1,\"Hi!\",π)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f31fc3ff",
- "metadata": {},
- "outputs": [],
- "source": [
- "showargs(6)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "dc47d83b",
- "metadata": {},
- "source": [
- "### Positional and keyword arguments\n",
- "\n",
- "Functions can combine positional and keyword arguments much like in Python, but keyword arguments start with semicolon `;` in Julia."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "37cd0314",
- "metadata": {},
- "outputs": [],
- "source": [
- "function foo(a,b;c,d)\n",
- " println(\"Positional: a=$a, b=$b. Keyword: c=$c, d=$d\")\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ee19db3d",
- "metadata": {},
- "outputs": [],
- "source": [
- "foo(3,4,d=2,c=1)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "4919ef95",
- "metadata": {},
- "source": [
- "### Optional arguments\n",
- "\n",
- "We can provide default values to arguments to make them optional."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "650e6d67",
- "metadata": {},
- "outputs": [],
- "source": [
- "function bar(a,b=0;c,d=1)\n",
- " println(\"Positional: a=$a, b=$b. Keyword: c=$c, d=$d\")\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e5d03536",
- "metadata": {},
- "outputs": [],
- "source": [
- "bar(1,c=2)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "94599fa3",
- "metadata": {},
- "source": [
- "
\n",
- "Question: Which will be the value of `x` below? \n",
- "
\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "79bd84bd",
- "metadata": {},
- "outputs": [],
- "source": [
- "function hofun(x)\n",
- " y -> x*y\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a9a7b774",
- "metadata": {},
- "outputs": [],
- "source": [
- "f2 = hofun(2)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f6606086",
- "metadata": {},
- "outputs": [],
- "source": [
- "x = f2(3)\n",
- "x"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "bc8e9bcf",
- "metadata": {},
- "source": [
- "## Arrays\n",
- "\n",
- "Julia supports multi-dimensional arrays. They are very similar to Numpy arrays in Python. Let's learn the basics of Julia arrays.\n",
- "\n",
- "\n",
- "### Array literals\n",
- "\n",
- "We can create (small) arrays from the given values using array literals.\n",
- "\n",
- "Next cell creates a vector with 3 integers. Note for Python users: there is no difference between vectors and lists in Julia."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "08652696",
- "metadata": {},
- "outputs": [],
- "source": [
- "vec = [1,2,3]"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "70f15e92",
- "metadata": {},
- "source": [
- "We can create a matrix as follows."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a3aacf71",
- "metadata": {},
- "outputs": [],
- "source": [
- "mat = [1 2 3 4\n",
- " 5 6 7 8\n",
- " 9 10 11 12]"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e096b198",
- "metadata": {},
- "source": [
- "### Array initialization\n",
- "\n",
- "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!\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9357edac",
- "metadata": {},
- "outputs": [],
- "source": [
- "zeros(4)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "dfc88af3",
- "metadata": {},
- "outputs": [],
- "source": [
- "zeros(Int,4)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "49e5a779",
- "metadata": {},
- "outputs": [],
- "source": [
- "ones(2,3)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "6bbde94a",
- "metadata": {},
- "outputs": [],
- "source": [
- "fill(5.0,3,4)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d2c44cce",
- "metadata": {},
- "outputs": [],
- "source": [
- "fill(add,3,4)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "495827bc",
- "metadata": {},
- "source": [
- "### Array comprehensions\n",
- "\n",
- "We can also create the items in the array using a loop within an array comprehension."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4b3c56b4",
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "squares = [ i^2 for i in 1:8 ]"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "cee954ce",
- "metadata": {},
- "source": [
- "### Indexing\n",
- "\n",
- "We can get and set the items of an array by indexing the array."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "df901c08",
- "metadata": {},
- "outputs": [],
- "source": [
- "squares[3]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "dc58cf4c",
- "metadata": {},
- "outputs": [],
- "source": [
- "squares[end]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "aa7936e0",
- "metadata": {},
- "outputs": [],
- "source": [
- "squares[2:4]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "6cf66ef0",
- "metadata": {},
- "outputs": [],
- "source": [
- "squares[4] = 16"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9125cbeb",
- "metadata": {},
- "outputs": [],
- "source": [
- "squares[2:3] = [4,9]"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f64021ab",
- "metadata": {},
- "source": [
- "### Immutable element type"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6650d6e7",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c79b297e",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = [10,11,12,13]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "2694f3d4",
- "metadata": {},
- "outputs": [],
- "source": [
- "a[2] = \"Hi!\""
- ]
- },
- {
- "cell_type": "markdown",
- "id": "92bae51d",
- "metadata": {},
- "source": [
- "### Arrays of any element type\n",
- "\n",
- "Arrays of fixed element type seem to be very rigid, right? Python list have not this limitation. However, we can use arrays of `Any` type, which are as flexible as Python lists, or even more since they can also contain functions. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d656b9ce",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = Any[10,11,12,13]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3d08e4ba",
- "metadata": {},
- "outputs": [],
- "source": [
- "a[3] = \"HI!\""
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "92bf2cab",
- "metadata": {},
- "outputs": [],
- "source": [
- "a"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1c72351a",
- "metadata": {},
- "source": [
- "### Loops\n",
- "\n",
- "The loop in next cell visits the elements in `a` one after the other. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "72dfc0f2",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = [10,20,30,40]\n",
- "for ai in a\n",
- " @show ai\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d4906337",
- "metadata": {},
- "source": [
- "This loop visits the integers from 1 to the length of the array and indexes the array at each of these integers."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "75f5ab89",
- "metadata": {},
- "outputs": [],
- "source": [
- "for i in 1:length(a)\n",
- " ai = a[i]\n",
- " @show (i,ai)\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e7abfdd6",
- "metadata": {},
- "source": [
- "This loop \"enumerates\" the items in the array."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "96dfab12",
- "metadata": {},
- "outputs": [],
- "source": [
- "for (i,ai) in enumerate(a)\n",
- " @show (i,ai)\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b4747037",
- "metadata": {},
- "source": [
- "### Arrays indices are 1-based by default\n",
- "\n",
- "Be aware of this if you are a C or Python user."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "694ec903",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = [10,20,30,40]\n",
- "a[0]"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a2ed5b01",
- "metadata": {},
- "source": [
- "### Slicing allocates a new array\n",
- "\n",
- "This is also different from Numpy in Python."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "961a82a3",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = [1 2 3\n",
- " 4 5 6]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "1c9f2162",
- "metadata": {},
- "outputs": [],
- "source": [
- "s = a[:,2]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4664eeb7",
- "metadata": {},
- "outputs": [],
- "source": [
- "s[2] = 0"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c356e67d",
- "metadata": {},
- "outputs": [],
- "source": [
- "a"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b75a2b83",
- "metadata": {},
- "source": [
- "### Array views\n",
- "\n",
- "If you want to modify the original array, use `view` instead."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3d8e2ab7",
- "metadata": {},
- "outputs": [],
- "source": [
- "v = view(a,:,2)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "329079bd",
- "metadata": {},
- "outputs": [],
- "source": [
- "v[1] = 0"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3b26a483",
- "metadata": {},
- "outputs": [],
- "source": [
- "a"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "fef8495e",
- "metadata": {},
- "source": [
- "## Exercises"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7d84ceaf",
- "metadata": {},
- "source": [
- "### Exercise 1\n",
- "\n",
- "Implement a function `ex1(a)` that finds the largest item in the array `a`. 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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8e6c7091",
- "metadata": {},
- "outputs": [],
- "source": [
- "# Implement here"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "370afcca",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Test\n",
- "arr = [3,4,7,3,1,7,2]\n",
- "@test ex1(arr) == (7,3)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7be19eb4",
- "metadata": {},
- "source": [
- "### Exercise 2\n",
- "\n",
- "Implement a function `ex2(f,g)` that takes two functions `f(x)` and `g(x)` and returns a new function `h(x)` representing the sum of `f` and `g`, i.e., `h(x)=f(x)+g(x)`."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "331e90b6",
- "metadata": {},
- "outputs": [],
- "source": [
- "# Implement here"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9e86c981",
- "metadata": {},
- "outputs": [],
- "source": [
- "h = ex2(sin,cos)\n",
- "xs = LinRange(0,2π,100)\n",
- "@test all(x-> h(x) == sin(x)+cos(x), xs)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "3324a991",
- "metadata": {},
- "source": [
- "### Exercise 3 (hard)\n",
- "\n",
- "Function `mandel` estimates if a given point `(x,y)` in the complex plane belongs to the [Mandelbrot set](https://en.wikipedia.org/wiki/Mandelbrot_set)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "0ac481c1",
- "metadata": {},
- "outputs": [],
- "source": [
- "function mandel(x,y,max_iters)\n",
- " z = Complex(x,y)\n",
- " c = z\n",
- " threshold=2\n",
- " for n in 1:max_iters\n",
- " if abs(z)>threshold\n",
- " return n-1\n",
- " end\n",
- " z = z^2 +c\n",
- " end\n",
- " max_iters\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "86ea2591",
- "metadata": {},
- "source": [
- "If the value of `mandel` is less than `max_iters`, the point is provably outside the Mandelbrot set. If `mandel` is equal to `max_iters`, then the point is provably inside the set. The larger `max_iters`, the better the quality of the estimate (the nicer will be your plot).\n",
- "\n",
- "Plot the value of function `mandel` for each pixel in a 2D grid of the box.\n",
- "\n",
- "$$(-1.7,0.7)\\times(-1.2,1.2).$$\n",
- "\n",
- "Use a grid resolution of at least 1000 points in each direction and `max_iters` at least 10. You can increase these values to get nicer plots. To plot the values use function `heatmap` from the Julia package `GLMakie`. Use `LinRange` to divide the horizontal and vertical axes into pixels. See the documentation of these functions for help. `GLMakie` is 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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "33137ad9",
- "metadata": {},
- "outputs": [],
- "source": [
- "# Implement here"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "357e0490",
- "metadata": {},
- "source": [
- "# License\n",
- "\n",
- "TODO: replace link to website\n",
- "\n",
- "This notebook is part of the course [Programming Large Scale Parallel Systems](http://localhost:8000/) at Vrije Universiteit Amsterdam and may be used under a [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) license."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f8d92f25",
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Julia 1.9.1",
- "language": "julia",
- "name": "julia-1.9"
- },
- "language_info": {
- "file_extension": ".jl",
- "mimetype": "application/julia",
- "name": "julia",
- "version": "1.9.1"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/docs/src/notebooks/julia_basics.md b/docs/src/notebooks/julia_basics.md
deleted file mode 100644
index baa6dea..0000000
--- a/docs/src/notebooks/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
-
-
-
- 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
- .
-
-
-
-
-```
-
-```@raw html
-
-
-```
diff --git a/docs/src/notebooks/julia_distributed.ipynb b/docs/src/notebooks/julia_distributed.ipynb
deleted file mode 100644
index 72db93d..0000000
--- a/docs/src/notebooks/julia_distributed.ipynb
+++ /dev/null
@@ -1,1335 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "id": "287c5272",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "### Programming large-scale parallel systems\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "2133c064",
- "metadata": {},
- "source": [
- "# Distributed computing in Julia\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a7b64d5a",
- "metadata": {},
- "source": [
- "## Contents\n",
- "\n",
- "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:\n",
- "\n",
- "- How to create Julia processes\n",
- "- How to execute code remotely\n",
- "- How to send and receive data\n",
- "\n",
- "With this knowledge you will be able to implement simple and complex parallel algorithms in Julia."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ec225103",
- "metadata": {},
- "source": [
- "
\n",
- "Note: Do not forget to execute the next cell before starting this notebook! \n",
- "
"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "aebf6482",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Printf\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",
- "q_1_check(answer) = answer_checker(answer,\"a\")\n",
- "q_2_check(answer) = answer_checker(answer,\"b\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "01af032c",
- "metadata": {},
- "source": [
- "## How to create Julia processes\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "036a25d7",
- "metadata": {},
- "source": [
- "### Adding processes locally\n",
- "\n",
- " 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.\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c6fed889",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Distributed"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d16faba9",
- "metadata": {},
- "outputs": [],
- "source": [
- "addprocs(3)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f07ac76c",
- "metadata": {},
- "source": [
- "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).\n",
- "\n",
- "\n",
- "
\n",
- "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.\n",
- "
\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d2e507dc",
- "metadata": {},
- "source": [
- "### Each process runs a separated Julia instance\n",
- "\n",
- "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.\n",
- " \n",
- "
\n",
- "\n",
- "
\n",
- " \n",
- " "
- ]
- },
- {
- "cell_type": "markdown",
- "id": "da2ae9aa",
- "metadata": {},
- "source": [
- "### Basic info about processes\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9214029b",
- "metadata": {},
- "outputs": [],
- "source": [
- "procs()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7a8498e6",
- "metadata": {},
- "outputs": [],
- "source": [
- "workers()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c6a8d8e4",
- "metadata": {},
- "outputs": [],
- "source": [
- "nprocs()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "6f8b2888",
- "metadata": {},
- "outputs": [],
- "source": [
- "nworkers()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8583dee0",
- "metadata": {},
- "outputs": [],
- "source": [
- "myid()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "17f1541f",
- "metadata": {},
- "outputs": [],
- "source": [
- "@everywhere println(myid())"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "44633d04",
- "metadata": {},
- "source": [
- "In previous cell, we have used the macro `@everywhere` that evaluates the given code on all processes. As a result, each process will print its own process id."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "2d104013",
- "metadata": {},
- "source": [
- "### Creating workers in other machines\n",
- "\n",
- "For large parallel computations, one typically needs to use different computers in parallel. Function `addprocs` also provides a low-level method to start workers in other machines. Next code example would create 3 workers in `server1` and 4 new workers in server `server2` (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). \n",
- "\n",
- "\n",
- "\n",
- "\n",
- "```julia\n",
- "using Distributed\n",
- "machines = [(\"user@server1\",3),(\"user@server2\",4)]\n",
- "addprocs(machines)\n",
- "```\n",
- "\n",
- "\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "3869f1f7",
- "metadata": {},
- "source": [
- "\n",
- "
\n",
- "\n",
- "
\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6c09b2e5",
- "metadata": {},
- "source": [
- "### Adding workers with ClusterManagers.jl\n",
- "\n",
- "Previous way of starting workers in other machines is very low-level. Happily, there is a Julia package called [ClusterManagers.jl](https://github.com/JuliaParallel/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).\n",
- "\n",
- "\n",
- "```julia\n",
- "using Distributed\n",
- "using ClusterManagers\n",
- "addprocs(SlurmManager(128), partition=\"debug\", t=\"00:5:00\")\n",
- "```\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7c470ed9",
- "metadata": {},
- "source": [
- "\n",
- "
\n",
- "\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8c8bc619",
- "metadata": {},
- "source": [
- "## Executing code remotely\n",
- "\n",
- "We have added new processes to our Julia session. Let's start using them!"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "76a8bf7a",
- "metadata": {},
- "source": [
- "### Function `remotecall`\n",
- "\n",
- "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 function `ones` to create a matrix locally.\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e2d0803a",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = ones(2,3)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a87fe170",
- "metadata": {},
- "source": [
- "The next cell does the same operation, but remotely on process 2. Note that `remotecall` takes the function we want to execute remotely, the process id where we want to execute it and, finally, the function arguments."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "78010eb4",
- "metadata": {},
- "outputs": [],
- "source": [
- "proc = 2\n",
- "ftr = remotecall(ones,proc,2,3)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b931f9c5",
- "metadata": {},
- "source": [
- "Note that `remotecall` does not return the result of the underlying function, but a `Future`. 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 use `fetch`."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e8163967",
- "metadata": {},
- "outputs": [],
- "source": [
- "fetch(ftr)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "613c4a20",
- "metadata": {},
- "source": [
- "### `remotecall` is asynchronous\n",
- "\n",
- "It is important to note that `remotecall` does 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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d918200b",
- "metadata": {},
- "outputs": [],
- "source": [
- "fun = (m,n) -> (sleep(10); ones(m,n))"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d5fba3bc",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "339e425b",
- "metadata": {},
- "outputs": [],
- "source": [
- "proc = 2\n",
- "ftr = remotecall(fun,proc,2,3)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8ff8dc16",
- "metadata": {},
- "outputs": [],
- "source": [
- "1+1"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "40f9d6c4",
- "metadata": {},
- "source": [
- "However, when fetching the result, the current process blocks waiting until the result is available in the remote process and arrives to its destination."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "76705f73",
- "metadata": {},
- "outputs": [],
- "source": [
- "fetch(ftr)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "75cb8b37",
- "metadata": {},
- "source": [
- "### Useful macro: `@spawnat`\n",
- "\n",
- "You have provably realized that in order to use `remotecall` we 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 `@spawnat` generates an auxiliary function from the given block of code and calls `remotecall` for us. For instance, the two following cells are equivalent."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4b550925",
- "metadata": {},
- "outputs": [],
- "source": [
- "@spawnat proc ones(2,3)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "2394cc7f",
- "metadata": {},
- "outputs": [],
- "source": [
- "fun = () -> ones(2,3)\n",
- "remotecall(fun,proc)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8631fd6f",
- "metadata": {},
- "source": [
- "### `@async` vs `@spawnat`\n",
- "\n",
- "The relation between `@async` and `@pawnat` is obvious. From the user perspective they work almost in the same way. However, `@async` generates a task that runs asynchronously in the current process, whereas `@spawnat` executes a task in a remote process in parallel. In both cases, the result is obtained using `fetch`. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7d681ed9",
- "metadata": {},
- "outputs": [],
- "source": [
- "tsk = @async begin\n",
- " sleep(3)\n",
- " zeros(2)\n",
- "end\n",
- "fetch(tsk)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c8141cff",
- "metadata": {},
- "outputs": [],
- "source": [
- "ftr = @spawnat :any begin\n",
- " sleep(3)\n",
- " zeros(2)\n",
- "end\n",
- "fetch(ftr)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "10899cd4",
- "metadata": {},
- "source": [
- "### Another usefull macro: `@fetchfrom`\n",
- "\n",
- "Macro `@fetchfrom` is the blocking version of `@spawnat`. It blocks and returns the corresponding result instead of a `Future` object. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f17da98e",
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "a = @fetchfrom proc begin\n",
- " sleep(3)\n",
- " zeros(2)\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e9d37079",
- "metadata": {},
- "source": [
- "## Data movement\n",
- "\n",
- "Data movement is a crucial part in distributed-memory computations and it is usually one of its main computational\n",
- "bottlenecks. Being aware of the data we are moving when using functions such as `remotecall` is 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."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6d6e04d9",
- "metadata": {},
- "source": [
- "### Data movement in `remotecall` / `fetch`\n",
- "\n",
- "When usig `remotecall` we send to the remote process a function and its arguments. In this example, we send function name `+` and matrices `a` and `b` to proc 4. When fetching the result we receive a copy of the matrix from proc 4."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "fd901fcc",
- "metadata": {},
- "outputs": [],
- "source": [
- "proc = 4\n",
- "a = rand(10,10)\n",
- "b = rand(10,10)\n",
- "ftr = remotecall(+,proc,a,b)\n",
- "fetch(ftr);"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "eda62267",
- "metadata": {},
- "source": [
- "### Implicit data movement\n",
- "\n",
- "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 `a` and `b` to 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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "fe9866ac",
- "metadata": {},
- "outputs": [],
- "source": [
- "proc = 4\n",
- "a = rand(10,10)\n",
- "b = rand(10,10)\n",
- "fun = () -> a+b\n",
- "ftr = remotecall(fun,proc)\n",
- "fetch(ftr);"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "3860df18",
- "metadata": {},
- "source": [
- "### Data movement with remote channels\n",
- "\n",
- "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 to `put!` runs asynchronously on process 4."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "98260fd3",
- "metadata": {},
- "outputs": [],
- "source": [
- "fun = ()->Channel{Int}()\n",
- "chnl = RemoteChannel(fun)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "0a729a00",
- "metadata": {},
- "outputs": [],
- "source": [
- "@spawnat 4 begin\n",
- " for i in 1:5\n",
- " put!(chnl,i)\n",
- " end\n",
- " close(chnl)\n",
- "end;"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c2ce01e0",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "82fbbb88",
- "metadata": {},
- "outputs": [],
- "source": [
- "take!(chnl)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1b3f8dcc",
- "metadata": {},
- "source": [
- "### Remote channels can be buffered\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "426274ac",
- "metadata": {},
- "outputs": [],
- "source": [
- "buffer_size = 2\n",
- "owner = 3\n",
- "fun = ()->Channel{Int}(buffer_size)\n",
- "chnl = RemoteChannel(fun,owner)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "db3b6fd5",
- "metadata": {},
- "outputs": [],
- "source": [
- "@spawnat 4 begin\n",
- " println(\"start\")\n",
- " for i in 1:5\n",
- " put!(chnl,i)\n",
- " println(\"I have put $i\")\n",
- " end\n",
- " close(chnl)\n",
- " println(\"stop\")\n",
- "end;"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7e85654e",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a8500649",
- "metadata": {},
- "outputs": [],
- "source": [
- "take!(chnl)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b8b411f1",
- "metadata": {},
- "source": [
- "### Remote channels are not iterable\n",
- "\n",
- "One main difference with respect to conventional channels is that remote channels cannot be iterated. Let's repeat the example above."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "18c60a68",
- "metadata": {},
- "outputs": [],
- "source": [
- "fun = ()->Channel{Int}()\n",
- "chnl = RemoteChannel(fun)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "73088996",
- "metadata": {},
- "outputs": [],
- "source": [
- "@spawnat 4 begin\n",
- " for i in 1:5\n",
- " put!(chnl,i)\n",
- " end\n",
- " close(chnl)\n",
- "end;"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "4bb35283",
- "metadata": {},
- "source": [
- "Now, try to iterate over the channel in a for loop. It will result in an error since channels are not iterable."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "1a9bcca0",
- "metadata": {},
- "outputs": [],
- "source": [
- "for j in chnl\n",
- " @show j\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8eaea2fa",
- "metadata": {},
- "source": [
- "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 the `catch` block and breaks the loop."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "72acd664",
- "metadata": {},
- "outputs": [],
- "source": [
- "while true\n",
- " try\n",
- " j = take!(chnl)\n",
- " @show j\n",
- " catch\n",
- " break\n",
- " end\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1a3986c9",
- "metadata": {},
- "source": [
- "## Questions"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "fd22b74b",
- "metadata": {},
- "source": [
- "\n",
- "
\n",
- "Question (Q1): How many integers are transferred between master and worker? Including both directions. \n",
- "
\n",
- "Question (Q4): Which value will be the value of `x` ? \n",
- "
\n",
- "\n",
- "Which value will be the value of `x` ?"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "b8678fd1",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = zeros(Int,3)\n",
- "proc = myid()\n",
- "@sync @spawnat proc a[2] = 2\n",
- "x = a[2]\n",
- "x"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9e985c61",
- "metadata": {},
- "source": [
- "## Remember: each process runs in a separated Julia instance\n",
- "\n",
- "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.\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "cdc07cba",
- "metadata": {},
- "source": [
- "### Functions are defined in a single process\n",
- "\n",
- "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_ones` is only being defined in the local process. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4544ca4c",
- "metadata": {},
- "outputs": [],
- "source": [
- "function sleep_ones(m,n)\n",
- " sleep(4)\n",
- " ones(m,n)\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "b54a1a84",
- "metadata": {},
- "outputs": [],
- "source": [
- "proc = 3\n",
- "remotecall_fetch(sleep_ones,proc,3,4)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9ec5d8fe",
- "metadata": {},
- "source": [
- "To fix this, we can define the function on all processes with the `@everywhere` macro."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7eb1e210",
- "metadata": {},
- "outputs": [],
- "source": [
- "@everywhere function sleep_ones(m,n)\n",
- " sleep(4)\n",
- " ones(m,n)\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "b55cbedf",
- "metadata": {},
- "outputs": [],
- "source": [
- "proc = 3\n",
- "remotecall_fetch(sleep_ones,proc,3,4)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9f71e7fd",
- "metadata": {},
- "source": [
- "### Anonymous functions are available everywhere\n",
- "\n",
- "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 `@everywhere` to work in a remote call.\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d9734913",
- "metadata": {},
- "outputs": [],
- "source": [
- "fun = (m,n) -> (sleep(4);ones(m,n))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f08131c3",
- "metadata": {},
- "outputs": [],
- "source": [
- "proc = 3\n",
- "remotecall_fetch(fun,proc,3,4)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ebec4f5e",
- "metadata": {},
- "source": [
- "### Each proc uses packages independently\n",
- "\n",
- "When using a package in a process, it is not available in the other ones. For instance, if we load the `LinearAlgebra` package in the current process and use one of its exported functions in another process, we will get an error."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e0cf3545",
- "metadata": {},
- "outputs": [],
- "source": [
- "using LinearAlgebra"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "0498748c",
- "metadata": {},
- "outputs": [],
- "source": [
- "@fetchfrom 3 norm([1,2,3])"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "97462cbd",
- "metadata": {},
- "source": [
- "To fix this, we can load the package on all processors with the `@everywhere` macro."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "14ee3498",
- "metadata": {},
- "outputs": [],
- "source": [
- "@everywhere using LinearAlgebra"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f8b3ba73",
- "metadata": {},
- "outputs": [],
- "source": [
- "@fetchfrom 3 norm([1,2,3])"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "3743ae18",
- "metadata": {},
- "source": [
- "### Each process has its own active package environment\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "abf3afe5",
- "metadata": {},
- "outputs": [],
- "source": [
- "] activate ."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "4d8bfa90",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3c6f29fb",
- "metadata": {},
- "outputs": [],
- "source": [
- "@everywhere using Pkg"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ee04b3cc",
- "metadata": {},
- "outputs": [],
- "source": [
- "@spawnat 2 Pkg.status();"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8244ae89",
- "metadata": {},
- "source": [
- "To fix this, you need to activate the current directory on all processes."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "99991a93",
- "metadata": {},
- "outputs": [],
- "source": [
- "@everywhere Pkg.activate(\".\")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "dfc8ea20",
- "metadata": {},
- "outputs": [],
- "source": [
- "@spawnat 2 Pkg.status();"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "4bfbe073",
- "metadata": {},
- "source": [
- "## Easy ways of parallelizing code\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "89e4b22b",
- "metadata": {},
- "source": [
- "### Useful macro: @distributed\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f5e27d24",
- "metadata": {},
- "outputs": [],
- "source": [
- "function compute_π(n)\n",
- " s = 1.0\n",
- " for i in 1:n\n",
- " s += (isodd(i) ? -1 : 1) / (i*2+1)\n",
- " end\n",
- " 4*s\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d0921cb9",
- "metadata": {},
- "source": [
- "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)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d4efb990",
- "metadata": {},
- "outputs": [],
- "source": [
- "function compute_π_dist(n)\n",
- " s = 1.0\n",
- " r = @distributed (+) for i in 1:n\n",
- " (isodd(i) ? -1 : 1) / (i*2+1)\n",
- " end\n",
- " 4*(s+r)\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8dbc2240",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ec604251",
- "metadata": {},
- "outputs": [],
- "source": [
- "@time compute_π(4_000_000_000)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c21b12b5",
- "metadata": {},
- "source": [
- "Run next cell to measure the performance of the parallel function."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "60fca9ae",
- "metadata": {},
- "outputs": [],
- "source": [
- "@time compute_π_dist(4_000_000_000)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f996ec0b",
- "metadata": {},
- "source": [
- "### Useful function: `pmap`\n",
- "\n",
- "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 `pmap` consider 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 for `pmap`."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3d137f35",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = [ rand(300,300) for i in 1:60];"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b95d3558",
- "metadata": {},
- "source": [
- "First, lets measure the serial performance"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9f6e302b",
- "metadata": {},
- "outputs": [],
- "source": [
- "using LinearAlgebra"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "74f16257",
- "metadata": {},
- "outputs": [],
- "source": [
- "@time svd.(a);"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "25de96e4",
- "metadata": {},
- "source": [
- "If we use `pmap` instead of broadcast, the different calls to `svd` will be distributed over the available processes."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4f9c5c2d",
- "metadata": {},
- "outputs": [],
- "source": [
- "@time pmap(svd,a);"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ad00cf8f",
- "metadata": {},
- "source": [
- "## Summary\n",
- "\n",
- "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 `remotecall` and `RemoteChannel`, but there are other functions and macros like `pmap` and `@distributed` that simplify the implementation of parallel algorithms."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9a49ad48",
- "metadata": {},
- "source": [
- "# License\n",
- "\n",
- "TODO: replace link to website\n",
- "\n",
- "This notebook is part of the course [Programming Large Scale Parallel Systems](http://localhost:8000/) at Vrije Universiteit Amsterdam and may be used under a [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) license."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8e36ae43",
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Julia 1.9.1",
- "language": "julia",
- "name": "julia-1.9"
- },
- "language_info": {
- "file_extension": ".jl",
- "mimetype": "application/julia",
- "name": "julia",
- "version": "1.9.1"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/docs/src/notebooks/julia_distributed.md b/docs/src/notebooks/julia_distributed.md
deleted file mode 100644
index d9bdc3e..0000000
--- a/docs/src/notebooks/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
-
-
-
- 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
- .
-
-
-
-
-```
-
-```@raw html
-
-
-```
diff --git a/docs/src/notebooks/julia_intro.ipynb b/docs/src/notebooks/julia_intro.ipynb
deleted file mode 100644
index 4034f11..0000000
--- a/docs/src/notebooks/julia_intro.ipynb
+++ /dev/null
@@ -1,746 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "id": "ae2a0512",
- "metadata": {},
- "source": [
- "\n",
- "# Is Julia fast?"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e040636a",
- "metadata": {},
- "source": [
- "## Contents\n",
- "\n",
- "With this notebook, you will learn\n",
- "\n",
- "- Two basic julia concepts related performance:\n",
- " - type-inference\n",
- " - JIT compilation\n",
- "- Some Julia syntax\n",
- "- Some useful Julia packages"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "beef4d5e",
- "metadata": {},
- "source": [
- "## Using Jupyter notebooks in Julia\n",
- "\n",
- "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.\n",
- "\n",
- "
\n",
- "Tip: Did you know that Jupyter stands for Julia, Python and R?\n",
- "
\n",
- "\n",
- "### How to start a Jupyter nootebook in Julia\n",
- "\n",
- "To run a Julia Jupyther notebook, open a Julia REPL and type\n",
- "\n",
- "```julia\n",
- "julia> ]\n",
- "pkg> add IJulia\n",
- "julia> using IJulia\n",
- "julia> notebook()\n",
- "```\n",
- "A new browser window will open. Navigate to the corresponding notebook and open it.\n",
- "\n",
- "
\n",
- "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.\n",
- "
\n",
- "\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7cdd5195",
- "metadata": {},
- "source": [
- "### Running a cell\n",
- "\n",
- "To run a cell, click on a cell and press `Shift` + `Enter`. You can also use the \"Run\" button in the toolbar above."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ec0f9041",
- "metadata": {},
- "outputs": [],
- "source": [
- "1+3\n",
- "4*5"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a8474251",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "fb3eaec4",
- "metadata": {},
- "outputs": [],
- "source": [
- "1+3\n",
- "4*5;"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0823876a",
- "metadata": {},
- "source": [
- "### Cell order is important\n",
- "\n",
- "Running the two cells below in reverse order won't work (try it). "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "bc1da6cb",
- "metadata": {},
- "outputs": [],
- "source": [
- "foo() = \"Well done!\""
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ddea57d8",
- "metadata": {},
- "outputs": [],
- "source": [
- "foo()"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "629189d9",
- "metadata": {},
- "source": [
- "### REPL modes\n",
- "\n",
- "This is particular to Julia notebooks. You can use package, help, and shell mode just like in the Julia REPL."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "defa851f",
- "metadata": {},
- "outputs": [],
- "source": [
- "] add MPI"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "1e6fbe66",
- "metadata": {},
- "outputs": [],
- "source": [
- "? print"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "123c0f6a",
- "metadata": {},
- "outputs": [],
- "source": [
- "; ls"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "44c03fde",
- "metadata": {},
- "source": [
- "## How fast is Julia code?"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f1e716db",
- "metadata": {},
- "source": [
- "NB. Most of the examples below are taken from the lecture by S.G. Johnson at MIT. See here:\n",
- "https://github.com/mitmath/18S096/blob/master/lectures/lecture1/Boxes-and-registers.ipynb"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6c088d9e",
- "metadata": {},
- "source": [
- "### Example\n",
- "\n",
- "Sum entries of a given array $a = [a_1,a_2,...,a_n]$\n",
- "\n",
- " $$s = \\sum_{i=1}^n a_i$$\n",
- "\n",
- "### "
- ]
- },
- {
- "cell_type": "markdown",
- "id": "34e3c7dd",
- "metadata": {},
- "source": [
- "### Hand-written sum function"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a672edac",
- "metadata": {
- "code_folding": []
- },
- "outputs": [],
- "source": [
- "function sum_hand(a)\n",
- " s = zero(eltype(a))\n",
- " for ai in a\n",
- " s += ai\n",
- " end\n",
- " s\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0494a1b3",
- "metadata": {},
- "source": [
- "### Test it\n",
- "\n",
- "The Julia macro `@test` which is provided in the `Test` package is useful to write (unit) tests in Julia."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ebac7b9e",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Test"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "aa5505d4",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = rand(5)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "fc11923c",
- "metadata": {},
- "outputs": [],
- "source": [
- "@test sum_hand(a) ≈ sum(a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e23f3ddf",
- "metadata": {},
- "source": [
- "## Benchmarking\n",
- "\n",
- "In Julia, the most straight-forward way of measuring the computation time of a piece of code is with the macro `@time`. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "caff0307",
- "metadata": {},
- "outputs": [],
- "source": [
- "a = rand(10^7);"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "6d073236",
- "metadata": {},
- "outputs": [],
- "source": [
- "@time sum_hand(a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "443ac8da",
- "metadata": {},
- "source": [
- "Note that `@time` also measures the compile time of a function if it's the first call to that function. So make sure to run `@time` twice on a freshly compiled function in order to get a more meaningful result."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c664522e",
- "metadata": {},
- "source": [
- "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 `@time` macro, only `@elapsed` discards the result of the computation and returns the elapsed time instead."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "bba864c1",
- "metadata": {},
- "outputs": [],
- "source": [
- "@elapsed sum_hand(a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "70317d33",
- "metadata": {},
- "source": [
- "## BenchmarkTools\n",
- "\n",
- "The `BenchmarkTools` extension package provides useful macros for sampling runtimes automatically. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "424ffe03",
- "metadata": {},
- "outputs": [],
- "source": [
- "using BenchmarkTools"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "079ef7d0",
- "metadata": {},
- "source": [
- "First of all, the `@benchmark` macro 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. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "b5bb9a12",
- "metadata": {},
- "outputs": [],
- "source": [
- "bch_sum_hand = @benchmark sum_hand($a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c77d454b",
- "metadata": {},
- "source": [
- "For quick sanity checks, one can use the `@btime` macro, which is a convenience wrapper around `@benchmark`. It returns only the minimum execution time and memory allocations. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "fee5493d",
- "metadata": {},
- "outputs": [],
- "source": [
- "@btime sum_hand($a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "3a1a774c",
- "metadata": {},
- "source": [
- "Similar to the `@elapsed` macro, `BenchmarkTool`'s `@belapsed` discards the return value of the function and instead returns the minimum runtime in seconds. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8e0cdb8c",
- "metadata": {},
- "outputs": [],
- "source": [
- "@belapsed sum_hand($a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "fb530b87",
- "metadata": {},
- "source": [
- "As opposed to `@time` and `@elapsed`, `@btime` and `@belapsed` run the code several times and return the minimum runtime, thus eliminating possible compilation times from the measurement. "
- ]
- },
- {
- "cell_type": "markdown",
- "id": "782be14f",
- "metadata": {},
- "source": [
- "### Built-in sum function"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "07df4939",
- "metadata": {},
- "outputs": [],
- "source": [
- "bch_sum = @benchmark sum($a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "18972fa7",
- "metadata": {},
- "source": [
- "### Hand-written sum in Python\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "038702c0",
- "metadata": {},
- "outputs": [],
- "source": [
- "using PyCall"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8e1dcbfa",
- "metadata": {},
- "outputs": [],
- "source": [
- "py\"\"\"\n",
- "def sum_py_hand(A):\n",
- " s = 0.0\n",
- " for a in A:\n",
- " s += a\n",
- " return s\n",
- "\"\"\"\n",
- "sum_py_hand = py\"sum_py_hand\""
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "60ca4517",
- "metadata": {},
- "outputs": [],
- "source": [
- "@test sum(a) ≈ sum_py_hand(a)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "2f009076",
- "metadata": {},
- "outputs": [],
- "source": [
- "bch_sum_py_hand = @benchmark sum_py_hand($a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "22d29afa",
- "metadata": {},
- "source": [
- "### Numpy sum "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "0f7ffb8e",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Conda"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e667fc80",
- "metadata": {},
- "outputs": [],
- "source": [
- "numpy = pyimport(\"numpy\")\n",
- "sum_numpy = numpy[\"sum\"]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ee0361df",
- "metadata": {},
- "outputs": [],
- "source": [
- "@test sum_numpy(a) ≈ sum(a)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d2f4e9ea",
- "metadata": {},
- "outputs": [],
- "source": [
- "bch_sum_numpy = @benchmark sum_numpy($a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d7c3cbd6",
- "metadata": {},
- "source": [
- "### Sumary of the results\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "5ce4c97a",
- "metadata": {},
- "outputs": [],
- "source": [
- "timings = [bch_sum_hand,bch_sum,bch_sum_py_hand,bch_sum_numpy]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9d5caabd",
- "metadata": {},
- "outputs": [],
- "source": [
- "methods = [\"sum_hand\",\"sum\",\"sum_py_hand\",\"sum_numpy\"]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8b7a3661",
- "metadata": {},
- "outputs": [],
- "source": [
- "using DataFrames"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "82b11610",
- "metadata": {},
- "outputs": [],
- "source": [
- "df = DataFrame(method=methods,time=timings)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f03c4281",
- "metadata": {},
- "source": [
- "### Improving the hand-written sum in Julia\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7696719c",
- "metadata": {},
- "outputs": [],
- "source": [
- "# ✍️ Exercise 3\n",
- "function sum_hand_fast(a)\n",
- " s = 0.0\n",
- " @simd for ai in a\n",
- " s += ai\n",
- " end\n",
- " s\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "0b461036",
- "metadata": {},
- "outputs": [],
- "source": [
- "@test sum_hand_fast(a) ≈ sum(a)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "fc228cba",
- "metadata": {},
- "outputs": [],
- "source": [
- "@benchmark sum_hand_fast($a)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "328f9128",
- "metadata": {},
- "source": [
- "## Conlcusions so far\n",
- "\n",
- "- Julia code (for loops) are much faster than in Python\n",
- "- Julia code can be as fast as optimized C code"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f9ce5464",
- "metadata": {},
- "source": [
- "## Why Julia is fast?\n",
- "\n",
- "- Julia is a compiled language (like C, C++, Fortran)\n",
- "- Julia is JIT compiled (C, C++, Fortran are AOT compiled)\n",
- "- Type declarations are optional in Julia\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "aa5254e4",
- "metadata": {},
- "source": [
- "# Conclusion: Why we use Julia in this course\n",
- "\n",
- "- Julia code is fast (it can be as fast as C)\n",
- "- Julia is a high-level language with simpler syntax than C \n",
- "- Julia supports different parallel programming models\n",
- "\n",
- "We will look into the third point in a later section of this course. \n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d4efafe2",
- "metadata": {},
- "source": [
- "## Solution to the exercises\n",
- "\n",
- "### Solution to Exercise 1"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "17f73c0f",
- "metadata": {},
- "outputs": [],
- "source": [
- "function sum_hand(a)\n",
- " s = 0.0\n",
- " for ai in a\n",
- " s += ai\n",
- " end\n",
- " s\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d12d6400",
- "metadata": {},
- "source": [
- "### Solution to Exercise 2"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e39094da",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Statistics\n",
- "\n",
- "a = rand(10^7)\n",
- "num_it = 15\n",
- "runtimes = zeros(num_it)\n",
- "for i in 1:num_it\n",
- " runtimes[i] = @elapsed sum_hand(a)\n",
- "end\n",
- "@show mean(runtimes) \n",
- "@show std(runtimes)\n",
- "@show minimum(runtimes)\n",
- "@show maximum(runtimes);"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7943f1cc",
- "metadata": {},
- "outputs": [],
- "source": [
- "# ✍️ Exercise 3\n",
- "function sum_hand_fast(a)\n",
- " s = 0.0\n",
- " @simd for ai in a\n",
- " s += ai\n",
- " end\n",
- " s\n",
- "end"
- ]
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Julia 1.9.0",
- "language": "julia",
- "name": "julia-1.9"
- },
- "language_info": {
- "file_extension": ".jl",
- "mimetype": "application/julia",
- "name": "julia",
- "version": "1.9.0"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/docs/src/notebooks/julia_intro.md b/docs/src/notebooks/julia_intro.md
deleted file mode 100644
index 2770451..0000000
--- a/docs/src/notebooks/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
-
-
-
- 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
- .
-
- 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
- .
-
-
-
-
-```
-
-```@raw html
-
-
-```
diff --git a/docs/src/notebooks/julia_tutorial.ipynb b/docs/src/notebooks/julia_tutorial.ipynb
deleted file mode 100644
index 1c98c21..0000000
--- a/docs/src/notebooks/julia_tutorial.ipynb
+++ /dev/null
@@ -1,523 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "### Programming large-scale parallel systems\n",
- "### Lectures on Julia for HPC\n",
- "\n",
- "\n",
- "# Tutorial: Getting started with Julia\n",
- "\n",
- "by Francesc Verdugo (VU Amsterdam)\n",
- "\n",
- "Version fall 2022"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Contents\n",
- "\n",
- "- Using the Julia REPL\n",
- "- Running serial and parallel code\n",
- "- Installing and managing packages"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Introduction\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Preliminaries\n",
- "\n",
- "This is a tutorial. To follow it:\n",
- "\n",
- "- Download and install Julia from https://julialang.org/\n",
- "- Download and install VSCode and its Julia extension. How to: https://www.julia-vscode.org/docs/dev/gettingstarted/\n",
- "- Continue reading this document \n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## The Julia REPL\n",
- "\n",
- "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`+`P` to open the command bar, and execute `Julia: 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:\n",
- "\n",
- "```\n",
- "julia> \n",
- "```\n",
- "\n",
- "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*.\n",
- "\n",
- "\n",
- "The usage is as follows:\n",
- "- You write some input\n",
- "- press enter\n",
- "- you get the output\n",
- "\n",
- "For instance, try this\n",
- "\n",
- "```julia\n",
- "julia> 1 + 1\n",
- "```"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Hello world!\n",
- "\n",
- "A \"Hello world\" example looks like this in Julia\n",
- "\n",
- "```julia\n",
- "julia> println(\"Hello, world!\")\n",
- "```\n",
- "\n",
- "Try to run it in the REPL."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Help mode\n",
- "\n",
- "\n",
- "\n",
- "Curious about what function `println` does? Enter into *help* mode to look into the documentation. This is done by typing a question mark (`?`) into the inut field:\n",
- "\n",
- "```julia\n",
- "julia> ?\n",
- "```\n",
- "\n",
- "After typing `?`, the command prompt changes to `help?>`. It means we are in help mode. Now, we can type a function name to see its documentation.\n",
- "\n",
- "```julia\n",
- "help?> println\n",
- "```"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Package and shell modes\n",
- "\n",
- "The REPL comes with two more modes, namely *package* and *shell* modes. To enter package mode type\n",
- "\n",
- "```julia\n",
- "julia> ]\n",
- "```\n",
- "\n",
- "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.\n",
- "\n",
- "To enter shell mode type semicolon (`;`)\n",
- "```julia\n",
- "julia> ;\n",
- "```\n",
- "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, \n",
- "\n",
- "```julia\n",
- "shell> ls\n",
- "```\n",
- "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.\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Running more complex code\n",
- "\n",
- "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 name `hello.jl`, execute it as follows\n",
- "\n",
- "```julia\n",
- "julia> include(\"hello.jl\")\n",
- "```\n",
- "\n",
- "
\n",
- "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 function pwd(). You can change to another directory with function cd() if needed. Also, make sure that the file extension is .jl.\n",
- "
\n",
- "\n",
- "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.\n",
- "\n",
- "```\n",
- "$ julia hello.jl\n",
- "```\n",
- "\n",
- "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 `PATH` environment variable. To check that Julia is properly installed, you can use\n",
- "\n",
- "```\n",
- "$ julia --version\n",
- "```\n",
- "\n",
- "If this runs without error and you see a version number, you are good to go!\n",
- "\n",
- "\n",
- "
\n",
- "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.\n",
- "
\n",
- "\n",
- "
\n",
- "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.\n",
- "
\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Running parallel code\n",
- "\n",
- "\n",
- "Since we are in a parallel computing course, let's run a parallel \"hello world\" example in Julia. Open a Julia REPL and write \n",
- "\n",
- "```julia\n",
- "julia> using Distributed\n",
- "julia> @everywhere println(\"Hello, world! I am proc $(myid()) from $(nprocs())\")\n",
- "```\n",
- "\n",
- "Here, we are using the `Distributed` package, 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.\n",
- "\n",
- "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 `addprocs` function.\n",
- "\n",
- "```julia\n",
- "julia> addprocs(3)\n",
- "```\n",
- "We have added 3 new processes, plus the old one, we have 4 processes. Run the code again.\n",
- "\n",
- "```julia\n",
- "julia> @everywhere println(\"Hello, world! I am proc $(myid()) from $(nprocs())\")\n",
- "```\n",
- "\n",
- "Now, you should see output from 4 processes.\n",
- "\n",
- "It is possible to specify the number of processes when starting Julia from the terminal with the `-p` argument (useful, e.g., when running in a cluster). If you launch Julia from the terminal as\n",
- "\n",
- "```\n",
- "$ julia -p 3\n",
- "```\n",
- "and then run \n",
- "\n",
- "```julia\n",
- "julia> @everywhere println(\"Hello, world! I am proc $(myid()) from $(nprocs())\")\n",
- "```\n",
- "\n",
- "You should get output from 4 processes as before.\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Installing packages\n",
- "\n",
- "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.\n",
- "\n",
- "\n",
- "Copy the following block of code into a new file named `\"hello_mpi.jl\"`\n",
- "\n",
- "```julia\n",
- "# file hello_mpi.jl\n",
- "using MPI\n",
- "MPI.Init()\n",
- "comm = MPI.COMM_WORLD\n",
- "rank = MPI.Comm_rank(comm)\n",
- "nranks = MPI.Comm_size(comm)\n",
- "println(\"Hello world, I am rank $rank of $nranks\")\n",
- "```\n",
- "\n",
- "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.\n",
- "\n",
- "Now, run the file from the REPL\n",
- "```julia\n",
- "julia> incude(\"hello_mpi.jl\")\n",
- "```\n",
- "\n",
- "It provably didn't work, right? Read the error message and note that the MPI package needs to be installed to run this code.\n",
- "\n",
- "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 `]`\n",
- "```julia\n",
- "julia> ]\n",
- "```\n",
- "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).\n",
- "\n",
- "To install the MPI package, type\n",
- "```julia\n",
- "(@v1.8) pkg> add MPI\n",
- "```\n",
- "Congrats, you have installed MPI!\n",
- "\n",
- "
\n",
- " 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 .jl needs to be ommited. In this case, we have isntalled the MPI.jl package even though we have only typed MPI in the REPL.\n",
- "
\n",
- "\n",
- "
\n",
- " Tip:\n",
- "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 of MPI.jl for further details.\n",
- "
\n",
- "\n",
- "To check that the package was installed properly, exit package mode by pressing the backspace key several times, and run it again\n",
- "\n",
- "```julia\n",
- "julia> incude(\"hello_mpi.jl\")\n",
- "```\n",
- "\n",
- "Now, it should work, but you provably get output from a single MPI rank only.\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Running MPI code\n",
- "\n",
- "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\n",
- "```\n",
- "$ mpiexec -np 4 julia hello_mpi.jl\n",
- "```\n",
- "But it will provably don't work since the version of `mpiexec` needs to match with the MPI version we are using from Julia. You can find the path to the `mpiexec` binary you need to use with these commands\n",
- "\n",
- "```julia\n",
- "julia> using MPI\n",
- "julia> MPI.mpiexec_path\n",
- "```\n",
- "\n",
- "and then try again\n",
- "```\n",
- "$ /path/to/my/mpiexec -np 4 julia hello_mpi.jl\n",
- "```\n",
- "with your particular path.\n",
- "\n",
- "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:\n",
- "```julia\n",
- "julia> using MPI\n",
- "julia> mpiexec(cmd->run(`$cmd -np 4 julia hello_mpi.jl`))\n",
- "```\n",
- "\n",
- "Now, you should see output from 4 ranks."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Installing packages locally\n",
- "\n",
- "We have installed the `MPI` package 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.\n",
- "\n",
- "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 `activate` command followed by the path to the folder you want to activate.\n",
- "```julia\n",
- "(@v1.8) pkg> activate .\n",
- "```\n",
- " Previous command will activate the current working directory. Note that the dot `.` is indeed the path to the current folder.\n",
- " \n",
- "The prompt has changed to `(lessons) pkg>` indicating that we are in the project within the `lessons` folder. The particular folder name can be different in your case.\n",
- "\n",
- "
\n",
- " Tip: You can activate a project directly when opening Julia from the terminal using the --project flag. 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 variable JULIA_PROJECT with the path of the folder you want to activate.\n",
- "
\n",
- " \n",
- " \n",
- "
\n",
- " Tip: The active project folder and the current working directory are two independent concepts! For instance, (@v1.8) pkg> activate folderB and then julia> cd(\"folderA\"), will activate the project in folderB and change the current working directory to folderA.\n",
- " \n",
- "
\n",
- "\n",
- "At this point all package-related operations will be local to the new project. For instance, install the `DataFrames` package.\n",
- "\n",
- "```julia\n",
- "(lessons) pkg> add DataFrames\n",
- "```\n",
- "Use the package to check that it is installed\n",
- "\n",
- "```julia\n",
- "julia> using DataFrames\n",
- "julia> DataFrame(a=[1,2],b=[3,4])\n",
- "```\n",
- "Now, we can return to the global project to check that `DataFrames` has not been installed there. To return to the global environment, use `activate` without a folder name.\n",
- "\n",
- "```julia\n",
- "(lessons) pkg> activate\n",
- "```\n",
- "The prompt is again `(@v1.8) pkg>`\n",
- "\n",
- "Now, try to use `DataFrames`.\n",
- "\n",
- "```julia\n",
- "julia> using DataFrames\n",
- "julia> DataFrame(a=[1,2],b=[3,4])\n",
- "```\n",
- "You should get an error or a warning unless you already had `DataFrames` installed globally.\n",
- "\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Project and Manifest files\n",
- "\n",
- "The information about a project is stored in two files `Project.toml` and `Manifest.toml`.\n",
- "\n",
- "- `Project.toml` contains the packages explicitly installed (the direct dependencies)\n",
- "\n",
- "- `Manifest.toml` contains direct and indirect dependencies along with the concrete version of each package.\n",
- "\n",
- "\n",
- "In other words, `Project.toml` contains the packages relevant for the user, whereas `Manifest.toml` is the detailed snapshot of all dependencies. The `Manifest.toml` can be used to reproduce the same envinonment in another machine.\n",
- "\n",
- "You can see the path to the current `Project.toml` file by using the `status` operator (or `st` in its short form) while in package mode\n",
- "\n",
- "```julia\n",
- "(@v1.8) pkg> status\n",
- "```\n",
- "\n",
- "The information about the `Manifest.toml` can be inspected by passing the `-m` flag.\n",
- "\n",
- "```julia\n",
- "(@v1.8) pkg> status -m\n",
- "```"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Installing packages from a project file\n",
- "\n",
- "Project files can be used to install lists of packages defined by others. E.g., to install all the dependencies of a Julia application.\n",
- "\n",
- "Assume that a colleague has sent to you a `Project.toml` file with this content:\n",
- "\n",
- "```\n",
- "[deps]\n",
- "BenchmarkTools = \"6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf\"\n",
- "DataFrames = \"a93c6f00-e57d-5684-b7b6-d8193f3e46c0\"\n",
- "MPI = \"da04e1cc-30fd-572f-bb4f-1f8673147195\"\n",
- "```\n",
- "\n",
- "Copy the contents of previous code block into a file called `Project.toml` and place it in an empty folder named `newproject`. It is important that the file is named `Project.toml`. You can create a new folder from the REPL with \n",
- "\n",
- "```julia\n",
- "julia> mkdir(\"newproject\")\n",
- "```\n",
- "\n",
- "To install all the packages registered in this file you need to activate the folder containing your `Project.toml` file\n",
- "```julia\n",
- "(@v1.8) pkg> activate newproject\n",
- "```\n",
- "and then *instantiating* it\n",
- "```julia\n",
- "(newproject) pkg> instantiate\n",
- "```\n",
- "\n",
- "The instantiate command will download and install all listed packages and their dependencies in just one click."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Getting help in package mode\n",
- "\n",
- "You can get help about a particular package operator by writing `help` in front of it\n",
- "\n",
- "```julia\n",
- "(@v1.8) pkg> help activate\n",
- "```\n",
- "\n",
- "You can get an overview of all package commands by typing `help` alone\n",
- "```julia\n",
- "(@v1.8) pkg> help\n",
- "```"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Package operations in Julia code\n",
- "\n",
- "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 `Pkg` package. For instance\n",
- "\n",
- "```julia\n",
- "julia> using Pkg\n",
- "julia> Pkg.status()\n",
- "```\n",
- "is equivalent to call `status` in package mode.\n",
- "```julia\n",
- "(@v1.8) pkg> status\n",
- "```"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Conclusion\n",
- "\n",
- "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\n",
- "\n",
- "- Julia Manual https://docs.julialang.org/en/v1/manual/getting-started/\n",
- "- Package manager https://pkgdocs.julialang.org/v1/getting-started/\n",
- "\n",
- "If you want to interact with the Julia community on discourse, sign in at https://discourse.julialang.org/"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# License\n",
- "\n",
- "TODO: replace link to website\n",
- "\n",
- "This notebook is part of the course [Programming Large Scale Parallel Systems](http://localhost:8000/) at Vrije Universiteit Amsterdam and may be used under a [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) license."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Julia 1.8.5",
- "language": "julia",
- "name": "julia-1.8"
- },
- "language_info": {
- "file_extension": ".jl",
- "mimetype": "application/julia",
- "name": "julia",
- "version": "1.8.5"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/docs/src/notebooks/julia_tutorial.md b/docs/src/notebooks/julia_tutorial.md
deleted file mode 100644
index 0e7850f..0000000
--- a/docs/src/notebooks/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
-
-
-
- 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
- .
-
\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a358ee60",
- "metadata": {},
- "source": [
- "### Goals\n",
- "\n",
- "We want to\n",
- "\n",
- "- compute the product in parallel using more than one process (distributed implementation)\n",
- "- study the performance of different parallelization alternatives\n",
- "- implement the algorithms using Julia\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "495ef679",
- "metadata": {},
- "source": [
- "### Assumptions\n",
- "\n",
- "- All matrices `A`,`B`, and `C` are initially stored in the master process\n",
- "- The result will be overwritten in `C`"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "5828e243",
- "metadata": {},
- "source": [
- "
\n",
- "\n",
- "
\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c57f340c",
- "metadata": {},
- "source": [
- "### Steps\n",
- "\n",
- "To develop and study the parallel implementation, we will follow these steps:\n",
- "\n",
- "- Identify the parts of the sequential algorithm that can be parallelized\n",
- "- Consider different parallelization strategies\n",
- "- Discuss the (theoretical) performance of these implementations\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ca56a7fe",
- "metadata": {},
- "source": [
- "## Serial implementation\n",
- "\n",
- "We start by considering the (naive) sequential algorithm:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "af8dfb37",
- "metadata": {},
- "outputs": [],
- "source": [
- "@everywhere function matmul_seq!(C,A,B)\n",
- " m = size(C,1)\n",
- " n = size(C,2)\n",
- " l = size(A,2)\n",
- " @assert size(A,1) == m\n",
- " @assert size(B,2) == n\n",
- " @assert size(B,1) == l\n",
- " z = zero(eltype(C))\n",
- " for j in 1:n\n",
- " for i in 1:m\n",
- " Cij = z\n",
- " for k in 1:l\n",
- " @inbounds Cij += A[i,k]*B[k,j]\n",
- " end\n",
- " C[i,j] = Cij\n",
- " end\n",
- " end\n",
- " C\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f967d2ea",
- "metadata": {},
- "source": [
- "
\n",
- "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.\n",
- "
\n",
- "\n",
- "Run the following cell to compare the performance of our hand-written function with respect to the built in function `mul!`\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "725387f6",
- "metadata": {},
- "outputs": [],
- "source": [
- "using LinearAlgebra\n",
- "N = 1000\n",
- "A = rand(N,N)\n",
- "B = rand(N,N)\n",
- "C = rand(N,N)\n",
- "@btime matmul_seq!(C,A,B)\n",
- "@btime mul!(C,A,B);"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0eedd28a",
- "metadata": {},
- "source": [
- "### Where do we can exploit parallelism?\n",
- "\n",
- "Look at the three nested loops in the sequential implementation:\n",
- "\n",
- "```julia\n",
- "for j in 1:n\n",
- " for i in 1:m\n",
- " Cij = z\n",
- " for k in 1:l\n",
- " @inbounds Cij += A[i,k]*B[k,j]\n",
- " end\n",
- " C[i,j] = Cij\n",
- " end\n",
- "end\n",
- "```\n",
- "- Loops over `i` and `j` are trivially parallelizable.\n",
- "- The loop over `k` can be parallelized but it requires a reduction."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b50aecff",
- "metadata": {},
- "source": [
- "### Parallel algorithms\n",
- "\n",
- "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:\n",
- "\n",
- "- Algorithm 1: each worker computes a single entry of C\n",
- "- Algorithm 2: each worker computes a single row of C\n",
- "- Algorithm 3: each worker computes a block rows of C"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6a706283",
- "metadata": {},
- "source": [
- "
\n",
- "\n",
- "
\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "bd759cb2",
- "metadata": {},
- "source": [
- "## Parallel algorithm 1"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "bdae0a02",
- "metadata": {},
- "source": [
- "### Data dependencies\n",
- "\n",
- "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.\n",
- "\n",
- "In algorithm 1, each worker computes only an entry of the result matrix C."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "acfb354b",
- "metadata": {},
- "source": [
- "
\n",
- "Question: Which are the data dependencies of the computations done by the worker in charge of computing entry C[i,j] ? \n",
- "
\n",
- "\n",
- " a) column A[:,i] and row B[j,:]\n",
- " b) row A[i,:] and column B[:,j]\n",
- " c) the whole matrices A and B\n",
- " d) row A[i,:] and the whole matrix B"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a8b7d1e1",
- "metadata": {},
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, c, or d \n",
- "alg_1_deps_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "06e1977a",
- "metadata": {},
- "source": [
- "### Implementation\n",
- "\n",
- "Taking into account the data dependencies, the parallel algorithm 1 can be efficiently implemented following these steps from the worker perspective:\n",
- "\n",
- "1. The worker receives the corresponding row A[i,:] and column B[:,j] from the master process\n",
- "2. The worker computes the dot product of A[i,:] and B[:,j]\n",
- "3. The worker sends back the result of C[i,j] to the master process"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "70087bce",
- "metadata": {},
- "source": [
- "
\n",
- "\n",
- "
\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9d22ccea",
- "metadata": {},
- "source": [
- "A possible implementation of this algorithm in Julia is as follows:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e4697fda",
- "metadata": {},
- "outputs": [],
- "source": [
- "function matmul_dist_1!(C, A, B)\n",
- " m = size(C,1)\n",
- " n = size(C,2)\n",
- " l = size(A,2)\n",
- " @assert size(A,1) == m\n",
- " @assert size(B,2) == n\n",
- " @assert size(B,1) == l\n",
- " z = zero(eltype(C))\n",
- " @assert nworkers() == m*n\n",
- " iw = 0 \n",
- " @sync for j in 1:n\n",
- " for i in 1:m\n",
- " Ai = A[i,:]\n",
- " Bj = B[:,j]\n",
- " iw += 1\n",
- " w = workers()[iw]\n",
- " ftr = @spawnat w begin\n",
- " Cij = z\n",
- " for k in 1:l\n",
- " @inbounds Cij += Ai[k]*Bj[k]\n",
- " end\n",
- " Cij\n",
- " end\n",
- " @async C[i,j] = fetch(ftr)\n",
- " end\n",
- " end\n",
- " C\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f0e5a38b",
- "metadata": {},
- "source": [
- "You can execute the following cells to test this implementation."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "13920a31",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Test\n",
- "N = 2\n",
- "A = rand(N,N)\n",
- "B = rand(N,N)\n",
- "C = similar(A)\n",
- "@test matmul_dist_1!(C,A,B) ≈ A*B"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f69d3333",
- "metadata": {},
- "source": [
- "### Performance\n",
- "\n",
- "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 \n",
- "\n",
- "$$\n",
- "S_P = \\frac{T_1}{T_P},\n",
- "$$\n",
- "\n",
- "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$:\n",
- "\n",
- "$$\n",
- "S^{*}_p = P.\n",
- "$$\n",
- "\n",
- "The ratio of the actual speedup over the optimal one is called the parallel efficiency\n",
- "\n",
- "$$\n",
- "E_p = \\frac{S_p}{S^{*}_p} = \\frac{T_1/T_P}{P}.\n",
- "$$\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b8eb224d",
- "metadata": {},
- "source": [
- "### Experimental speedup\n",
- "\n",
- "The following cell measures the speedup of parallel algorithm 1. Do we achieve the optimal speedup?"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "cc698aa8",
- "metadata": {},
- "outputs": [],
- "source": [
- "N = 2\n",
- "A = rand(N,N)\n",
- "B = rand(N,N)\n",
- "C = similar(A)\n",
- "T1 = @belapsed matmul_seq!(C,A,B)\n",
- "C = similar(A)\n",
- "TP = @belapsed matmul_dist_1!(C,A,B)\n",
- "P = nworkers()\n",
- "println(\"Speedup = \", T1/TP)\n",
- "println(\"Optimal speedup = \", P)\n",
- "println(\"Efficiency = \", 100*(T1/TP)/P, \"%\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "dac6a50b",
- "metadata": {},
- "source": [
- "### Communication overhead\n",
- "\n",
- "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.\n",
- "\n",
- "Remember, algorithm 1 consisted of these main steps:\n",
- "\n",
- "1. The worker receives the corresponding row A[i,:] and column B[:,j] from the master process\n",
- "2. The worker computes the dot product of A[i,:] and B[:,j]\n",
- "3. The worker sends back the result of C[i,j] to the master process\n",
- "\n",
- "
\n",
- "Question: How many scalars are communicated from and to a worker? Assume that matrices A, B, and C are N by N matrices.\n",
- "
\n",
- "Question: How many operations are done in a worker? \n",
- "
\n",
- "\n",
- " a) O(N)\n",
- " b) O(N²)\n",
- " c) O(N³)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "fcc9b903",
- "metadata": {},
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, or c\n",
- "alg_1_comp_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d4c301de",
- "metadata": {},
- "source": [
- "From these results we can conclude:\n",
- "\n",
- "- The communication complexity is O(N)\n",
- "- The computation complexity is O(N)\n",
- "- The ratio communication over computation (the communication overhead) is O(1)\n",
- "\n",
- "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.\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b15cbaf4",
- "metadata": {},
- "source": [
- "## Parallel algorithm 2\n",
- "\n",
- "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.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "62e5c637",
- "metadata": {},
- "source": [
- "
\n",
- "Question: Which are the data dependencies of the computations done by the worker in charge of computing row C[i,:] ? \n",
- "
\n",
- "\n",
- " a) column A[:,i] and row B[j,:]\n",
- " b) row A[i,:] and column B[:,j]\n",
- " c) the whole matrices A and B\n",
- " d) row A[i,:] and the whole matrix B"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "cdb46cd8",
- "metadata": {},
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, c, or d \n",
- "alg_2_deps_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a9d84ac2",
- "metadata": {},
- "source": [
- "### Implementation\n",
- "\n",
- "These are the main steps of the implementation of algorithm 2:\n",
- "\n",
- "1. The worker receives the corresponding row A[i,:] and matrix B from the master process\n",
- "2. The worker computes the product of row A[i,:] times B\n",
- "3. The worker sends back the result of row C[i,:] to the master process\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "fb6b572b",
- "metadata": {},
- "source": [
- "\n",
- "
\n",
- "\n",
- "
\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a662afa4",
- "metadata": {},
- "source": [
- "A possible implementation of this algorithm in Julia is as follows:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "365dc58e",
- "metadata": {},
- "outputs": [],
- "source": [
- "function matmul_dist_2!(C, A, B)\n",
- " m = size(C,1)\n",
- " n = size(C,2)\n",
- " l = size(A,2)\n",
- " @assert size(A,1) == m\n",
- " @assert size(B,2) == n\n",
- " @assert size(B,1) == l\n",
- " z = zero(eltype(C))\n",
- " @assert nworkers() == m\n",
- " iw = 0\n",
- " @sync for i in 1:m\n",
- " Ai = A[i,:]\n",
- " iw += 1\n",
- " w = workers()[iw]\n",
- " ftr = @spawnat w begin\n",
- " Ci = fill(z,l)\n",
- " for j in 1:n\n",
- " for k in 1:l\n",
- " @inbounds Ci[j] += Ai[k]*B[k,j]\n",
- " end\n",
- " end\n",
- " Ci\n",
- " end\n",
- " @async C[i,:] = fetch(ftr)\n",
- " end\n",
- " C\n",
- " end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8de835b9",
- "metadata": {},
- "source": [
- "Test it using next cell"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "267ac8b2",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Test\n",
- "N = 4\n",
- "A = rand(N,N)\n",
- "B = rand(N,N)\n",
- "C = similar(A)\n",
- "@test matmul_dist_2!(C,A,B) ≈ A*B"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f1f30faf",
- "metadata": {},
- "source": [
- "### Experimental speedup"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "fe42f069",
- "metadata": {},
- "outputs": [],
- "source": [
- "N = 4\n",
- "A = rand(N,N)\n",
- "B = rand(N,N)\n",
- "C = similar(A)\n",
- "T1 = @belapsed matmul_seq!(C,A,B)\n",
- "C = similar(A)\n",
- "TP = @belapsed matmul_dist_2!(C,A,B)\n",
- "P = nworkers()\n",
- "println(\"Speedup = \", T1/TP)\n",
- "println(\"Optimal speedup = \", P)\n",
- "println(\"Efficiency = \", 100*(T1/TP)/P, \"%\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e2c6f60a",
- "metadata": {},
- "source": [
- "### Complexity\n",
- "\n",
- "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:\n",
- "\n",
- "1. The worker receives the corresponding row A[i,:] and matrix B from the master process\n",
- "2. The worker computes the product of row A[i,:] times B\n",
- "3. The worker sends back the result of row C[i,:] to the master process\n",
- "\n",
- "
\n",
- "Question: Which is the complexity of the communication and computations done by a worker in algorithm 2?\n",
- "
\n",
- "\n",
- " a) O(N) communication and O(N^2) computation\n",
- " b) O(N^2) communication and O(N^2) computation\n",
- " c) O(N^3) communication and O(N^3) computation\n",
- " d) O(N) communication and O(N) computation"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "1bf8feff",
- "metadata": {},
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, c, or d \n",
- "alg_2_complex_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a2038e04",
- "metadata": {},
- "source": [
- "The communication and computation cost are still of the same order of magnitude even though we have increased the grain size. "
- ]
- },
- {
- "cell_type": "markdown",
- "id": "71088fb9",
- "metadata": {},
- "source": [
- "## Parallel algorithm 3\n",
- "\n",
- "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.\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f1b8c712",
- "metadata": {},
- "source": [
- "
\n",
- "Question: Which are the data dependencies of the computations done by the worker in charge of computing the range of rows C[rows,:] ? \n",
- "
\n",
- "\n",
- "\n",
- " a) A[rows,:] and B[:,rows] \n",
- " b) the whole matrix A and B[:,rows]\n",
- " c) A[rows,:] and the whole matrix B\n",
- " d) the whole matrices A and B"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "2825385a",
- "metadata": {},
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, c, or d \n",
- "alg_3_deps_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "429faa32",
- "metadata": {},
- "source": [
- "### Implementation\n",
- "\n",
- "These are the main steps of the implementation of algorithm 3:\n",
- "\n",
- "1. The worker receives the corresponding rows A[rows,:] and matrix B from the master process\n",
- "2. The worker computes the product of A[rows,:] times B\n",
- "3. The worker sends back the result of C[rows,:] to the master process"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c14ebcb3",
- "metadata": {},
- "source": [
- "
\n",
- "\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7d7952c1",
- "metadata": {},
- "source": [
- "The implementation of this variant is let as an exercise (see below)."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0323d6d8",
- "metadata": {},
- "source": [
- "### Communication overhead\n",
- "\n",
- "Let us analyze the (theoretical) communication overhead for algorithm 3.\n",
- "\n",
- "
\n",
- "Question: Which is the complexity of the communication and computations done by a worker in algorithm 3?\n",
- "
\n",
- "\n",
- " a) O(N²) communication and O(N³) computation\n",
- " b) O(N²) communication and O(N³/P) computation\n",
- " c) O(N²) communication and O(N²) computation\n",
- " d) O(N²/P) communication and O(N³/P) computation"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "50b8bf53",
- "metadata": {},
- "outputs": [],
- "source": [
- "answer = \"x\" # replace x with a, b, c, or d \n",
- "alg_3_complex_check(answer)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "3f0d99e6",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7bb65f2e",
- "metadata": {},
- "source": [
- "## Summary\n",
- "\n",
- "The table below compares the three parallel algorithms. \n",
- "\n",
- "\n",
- "\n",
- "| Algorithm | Parallelism (#workers) | Communication per worker | Computation per worker | Ratio communication/ computation |\n",
- "|---|---|---|---|---|\n",
- "| 1 | N² | 2N + 1 | N | O(1) |\n",
- "| 2 | N | 2N + N² | N² | O(1) |\n",
- "| 3 | P | N² + 2N²/P | N³/P | O(P/N) |\n",
- "\n",
- "\n",
- "- Matrix-matrix multiplication is trivially parallelizable (all entries in the result matrix can be computed in parallel, at least in theory)\n",
- "- However, we cannot exploit all the potential parallelism in a distributed system due to communication overhead\n",
- "- We need a sufficiently large grain size to obtain a near optimal speedup\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8b83e744",
- "metadata": {},
- "source": [
- "## Exercises"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a628a1df",
- "metadata": {},
- "source": [
- "### Exercise 1\n",
- "\n",
- "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.\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8e50b923",
- "metadata": {},
- "outputs": [],
- "source": [
- "function matmul_dist_3!(C,A,B)\n",
- " m = size(C,1)\n",
- " n = size(C,2)\n",
- " l = size(A,2)\n",
- " @assert size(A,1) == m\n",
- " @assert size(B,2) == n\n",
- " @assert size(B,1) == l\n",
- " @assert mod(m,nworkers()) == 0\n",
- " # Implement here\n",
- " \n",
- " C\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "4506dcfb",
- "metadata": {},
- "source": [
- "Use test-driven development to implement the algorithm. Use this test:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "28cde36a",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Test\n",
- "P = nworkers()\n",
- "load = 100\n",
- "N = load*P\n",
- "A = rand(N,N)\n",
- "B = rand(N,N)\n",
- "C = similar(A)\n",
- "@test matmul_dist_3!(C,A,B) ≈ A*B"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "03952b0b",
- "metadata": {},
- "source": [
- "Measure the performance of your implementation by running next cell. Do you get close to the optimal speedup?"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "b3aa2b7c",
- "metadata": {},
- "outputs": [],
- "source": [
- "P = nworkers()\n",
- "load = 100\n",
- "N = load*P\n",
- "A = rand(N,N)\n",
- "B = rand(N,N)\n",
- "C = similar(A)\n",
- "T1 = @belapsed matmul_seq!(C,A,B)\n",
- "C = similar(A)\n",
- "TP = @belapsed matmul_dist_3!(C,A,B)\n",
- "println(\"Speedup = \", T1/TP)\n",
- "println(\"Optimal speedup = \", P)\n",
- "println(\"Efficiency = \", 100*(T1/TP)/P, \"%\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "fa8d7f40",
- "metadata": {},
- "source": [
- "### Exercise 2"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0e7c607e",
- "metadata": {},
- "source": [
- "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.\n",
- "See the following code:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "023b20d1",
- "metadata": {},
- "outputs": [],
- "source": [
- "function matmul_dist_1_v2!(C, A, B)\n",
- " m = size(C,1)\n",
- " n = size(C,2)\n",
- " l = size(A,2)\n",
- " @assert size(A,1) == m\n",
- " @assert size(B,2) == n\n",
- " @assert size(B,1) == l\n",
- " z = zero(eltype(C))\n",
- " @sync for j in 1:n\n",
- " for i in 1:m\n",
- " Ai = A[i,:]\n",
- " Bj = B[:,j]\n",
- " ftr = @spawnat :any begin\n",
- " Cij = z\n",
- " for k in 1:l\n",
- " @inbounds Cij += Ai[k]*Bj[k]\n",
- " end\n",
- " Cij\n",
- " end\n",
- " @async C[i,j] = fetch(ftr)\n",
- " end\n",
- " end\n",
- " C\n",
- "end"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "52005ca1",
- "metadata": {},
- "source": [
- "With this new implementation, we can multiply matrices of arbitrary size with a fixed number of workers. Test it:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c1d3595b",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Test\n",
- "N = 50\n",
- "A = rand(N,N)\n",
- "B = rand(N,N)\n",
- "C = similar(A)\n",
- "@test matmul_dist_1_v2!(C,A,B) ≈ A*B"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ab609c18",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d7d31710",
- "metadata": {},
- "outputs": [],
- "source": [
- "N = 100\n",
- "A = rand(N,N)\n",
- "B = rand(N,N)\n",
- "C = similar(A)\n",
- "P = nworkers()\n",
- "T1 = @belapsed matmul_seq!(C,A,B)\n",
- "C = similar(A)\n",
- "TP = @belapsed matmul_dist_1_v2!(C,A,B)\n",
- "println(\"Speedup = \", T1/TP)\n",
- "println(\"Optimal speedup = \", P)\n",
- "println(\"Efficiency = \", 100*(T1/TP)/P, \"%\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8e171362",
- "metadata": {},
- "source": [
- "# License\n",
- "\n",
- "TODO: replace link to website\n",
- "\n",
- "This notebook is part of the course [Programming Large Scale Parallel Systems](http://localhost:8000/) at Vrije Universiteit Amsterdam and may be used under a [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) license."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "86b7b044",
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Julia 1.9.1",
- "language": "julia",
- "name": "julia-1.9"
- },
- "language_info": {
- "file_extension": ".jl",
- "mimetype": "application/julia",
- "name": "julia",
- "version": "1.9.1"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/docs/src/notebooks/matrix_matrix.md b/docs/src/notebooks/matrix_matrix.md
deleted file mode 100644
index a1d4690..0000000
--- a/docs/src/notebooks/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
-
-
-
- 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
- .
-
-
-
-
-```
-
-```@raw html
-
-
-```
diff --git a/docs/src/notebooks/mpi_tutorial.ipynb b/docs/src/notebooks/mpi_tutorial.ipynb
deleted file mode 100644
index e142a42..0000000
--- a/docs/src/notebooks/mpi_tutorial.ipynb
+++ /dev/null
@@ -1,413 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "id": "8d800917",
- "metadata": {},
- "source": [
- "# Tutorial: Using MPI in Julia\n",
- "Message Passing Interface (MPI) is a standardized and portable library specification for communication between parallel processes in distributed memory systems. Julia offers a convenient way to work with MPI for creating efficient parallel and distributed applications. In this tutorial, you will learn how to use MPI from Julia to perform parallel computing tasks.\n",
- "\n",
- "## MPI launches separate Julia instances\n",
- "When you run an MPI-enabled Julia script, MPI takes care of spawning multiple instances of the Julia executable, each acting as a separate process. These workers can communicate with each other using MPI communication functions. This enables parallel processing and distributed computation. Here's a summary of how it works:\n",
- "\n",
- "-- TODO: insert picture here --\n",
- "\n",
- "- **MPI Spawns Processes**: The `mpiexec` command launches multiple instances of the Julia executable, creating separate worker processes. In this example, 4 Julia workers are spawned.\n",
- "\n",
- "- **Worker Communication**: These workers can communicate with each other using MPI communication functions, allowing them to exchange data and coordinate actions.\n",
- "\n",
- "- **Parallel Tasks**: The workers execute parallel tasks simultaneously, working on different parts of the computation to potentially speed up the process.\n",
- "\n",
- "\n",
- " \n",
- "\n",
- "\n",
- "## Installing MPI.jl and MPIClusterManagers Packages\n",
- "To use MPI in Julia, you'll need the MPI.jl package, and if you intend to run MPI programs in a Jupyter Notebook, you'll also need the MPIClusterManagers package. These packages provide the necessary bindings to the MPI library and cluster management capabilities. To install the packages, open a terminal and run the following commands:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3cb5f151",
- "metadata": {},
- "outputs": [],
- "source": [
- "using Pkg\n",
- "Pkg.add(\"MPI\")\n",
- "Pkg.add(\"MPIClusterManagers\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "ed45a4b2",
- "metadata": {},
- "source": [
- "
\n",
- " Tip:\n",
- "The package MPI.jl is the Julia interface to MPI. 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 of MPI.jl for further details.\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7a36916e",
- "metadata": {},
- "source": [
- "## Writing a HelloWorld MPI Program in Julia\n",
- "Let's start by creating a simple MPI program that prints a message along with the rank of each worker. \n",
- "\n",
- "Create a new Julia script, for example, `mpi_hello_world.jl`:\n",
- "\n",
- "```julia\n",
- "using MPI\n",
- "\n",
- "# Initialize MPI\n",
- "MPI.Init()\n",
- "\n",
- "# Get the default communicator (MPI_COMM_WORLD) for all processes\n",
- "comm = MPI.COMM_WORLD\n",
- "\n",
- "# Get the number of processes in this communicator\n",
- "nranks = MPI.Comm_size(comm)\n",
- "\n",
- "# Get the rank of the current process within the communicator\n",
- "rank = MPI.Comm_rank(comm)\n",
- "\n",
- "# Print a message with the rank of the current process\n",
- "println(\"Hello, I am process $rank of $nranks processes!\")\n",
- "\n",
- "# Finalize MPI\n",
- "MPI.Finalize()\n",
- "```"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6caa8d74",
- "metadata": {},
- "source": [
- "### MPI Communicators\n",
- "In MPI, a **communicator** is a context in which a group of processes can communicate with each other. `MPI_COMM_WORLD` is one of the MPI standard communicators, it represents all processes in the MPI program. Custom communicators can also be created to group processes based on specific requirements or logical divisions. \n",
- "\n",
- "The **rank** of a processor is a unique identifier assigned to each process within a communicator. It allows processes to distinguish and address each other in communication operations. "
- ]
- },
- {
- "cell_type": "markdown",
- "id": "19f41e38",
- "metadata": {},
- "source": [
- "## Running the HelloWorld MPI Program\n",
- "\n",
- "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\n",
- "```\n",
- "$ mpiexec -np 4 mpi_hello_world.jl\n",
- "```\n",
- "In this command, `-np 4` specifies the desired number of processes. \n",
- "But it will probably not work since the version of `mpiexec` needs to match with the MPI version we are using from Julia. You can find the path to the `mpiexec` binary you need to use with these commands\n",
- "\n",
- "```julia\n",
- "julia> using MPI\n",
- "julia> MPI.mpiexec_path\n",
- "```\n",
- "\n",
- "and then try again\n",
- "```\n",
- "$ /path/to/my/mpiexec -np 4 julia mpi_hello_world.jl\n",
- "```\n",
- "with your particular path.\n",
- "\n",
- "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:\n",
- "```julia\n",
- "julia> using MPI\n",
- "julia> mpiexec(cmd->run(`$cmd -np 4 julia mpi_hello_world.jl`))\n",
- "```\n",
- "\n",
- "Now, you should see output from 4 ranks.\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0592e58c",
- "metadata": {},
- "source": [
- "## Running MPI Programs in Jupyter Notebook with MPIClusterManagers\n",
- "If you want to run your MPI code from a Jupyter Notebook, you can do so using the `MPIClusterManagers` package.\n",
- "\n",
- "1. Load the packages and start an MPI cluster with the desired number of workers:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "cf66dd39",
- "metadata": {},
- "outputs": [],
- "source": [
- "using MPIClusterManagers\n",
- "# Distributed package is needed for addprocs()\n",
- "using Distributed\n",
- "\n",
- "manager = MPIWorkerManager(4)\n",
- "addprocs(manager)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d40fe3ee",
- "metadata": {},
- "source": [
- "2. Run your MPI code inside a `@mpi_do` block to execute it on the cluster workers:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "0a51d1f2",
- "metadata": {},
- "outputs": [],
- "source": [
- "@mpi_do manager begin\n",
- " using MPI\n",
- " comm = MPI.COMM_WORLD\n",
- " rank = MPI.Comm_rank(comm)\n",
- " println(\"Hello from process $rank\")\n",
- "end\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "38ed88c1",
- "metadata": {},
- "source": [
- "MPI is automatically initialized and finalized within the `@mpi_do` block.\n",
- "\n",
- "3. Remove processes when done:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e0b53cc1",
- "metadata": {},
- "outputs": [],
- "source": [
- "rmprocs(manager)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "5466a650",
- "metadata": {},
- "source": [
- "## Point-to-Point Communication with MPI\n",
- "MPI provides point-to-point communication using blocking send and receiving functions `MPI.send`, `MPI.recv`; or their non-blocking versions `MPI.Isend`, and `MPI.Irecv!`. These functions allow individual processes to send and receive data between each other.\n",
- "\n",
- "### Blocking communication\n",
- "\n",
- "Let's demonstrate how to send and receive with an example:\n",
- "\n",
- "```julia\n",
- "using MPI\n",
- "\n",
- "MPI.Init()\n",
- "\n",
- "comm = MPI.COMM_WORLD\n",
- "rank = MPI.Comm_rank(comm)\n",
- "\n",
- "# Send and receive messages using blocking MPI.send and MPI.recv\n",
- "if rank == 0\n",
- " data = \"Hello from process $rank !\"\n",
- " MPI.send(data, comm, dest=1)\n",
- "elseif rank == 1\n",
- " received_data = MPI.recv(comm, source=0)\n",
- " println(\"Process $rank received: $received_data\")\n",
- "end\n",
- "\n",
- "MPI.Finalize()\n",
- "```"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d4dfe654",
- "metadata": {},
- "source": [
- "In this example, process 0 sends a message using `MPI.send`, and process 1 receives it using `MPI.recv`.\n",
- "\n",
- "### Non-blocking communication\n",
- "\n",
- "To demonstrate asynchronous communication, let's modify the example using `MPI.Isend` and `MPI.Irecv!`:\n",
- "\n",
- "```julia\n",
- "using MPI\n",
- "\n",
- "MPI.Init()\n",
- "\n",
- "comm = MPI.COMM_WORLD\n",
- "rank = MPI.Comm_rank(comm)\n",
- "\n",
- "# Asynchronous communication using MPI.Isend and MPI.Irecv!\n",
- "if rank == 0\n",
- " data = \"Hello from process $rank !\"\n",
- " request = MPI.Isend(data, comm, dest=1)\n",
- " # Other computation can happen here\n",
- " MPI.Wait(request)\n",
- "elseif rank == 1\n",
- " received_data = Array{UInt8}(undef, 50) # Preallocate buffer\n",
- " request = MPI.Irecv!(received_data, comm, source=0)\n",
- " # Other computation can happen here\n",
- " MPI.Wait(request)\n",
- " println(\"Process $rank received: $(String(received_data))\")\n",
- "end\n",
- "\n",
- "MPI.Finalize()\n",
- "```"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "024db538",
- "metadata": {},
- "source": [
- "In this example, process 0 uses `MPI.Isend` to send the message asynchronously. This function returns immediately, allowing the sender process to continue its execution. However, the actual sending of data is done asynchronously in the background. Similar to `MPI.Isend`, `MPI.Irecv!` returns immediately, allowing the receiver process to continue executing. \n",
- "\n",
- "
\n",
- "Important: In asynchronous communication, always use MPI.Wait() to ensure the communication is finished before accessing the send or receive buffer.\n",
- "
\n",
- "\n",
- "\n",
- "## Collective Communication with MPI\n",
- "MPI provides collective communication functions for communication involving multiple processes. Let's explore some of these functions:\n",
- "\n",
- "- MPI.Gather: Gathers data from all processes to a single process.\n",
- "- MPI.Scatter: Distributes data from one process to all processes.\n",
- "- MPI.Bcast: Broadcasts data from one process to all processes.\n",
- "- MPI.Barrier: Synchronizes all processes.\n",
- "\n",
- "\n",
- "Let's illustrate the use of `MPI.Gather` and `MPI.Scatter` with an example:\n",
- "\n",
- "```julia\n",
- "# TODO: check if this runs correctly\n",
- "using MPI\n",
- "using Random\n",
- "\n",
- "MPI.Init()\n",
- "\n",
- "comm = MPI.COMM_WORLD\n",
- "rank = MPI.Comm_rank(comm)\n",
- "size = MPI.Comm_size(comm)\n",
- "\n",
- "# Root processor generates random data\n",
- "data = rand(rank == 0 ? size * 2 : 0)\n",
- "\n",
- "# Scatter data to all processes\n",
- "local_data = Vector{Float64}(undef, 2)\n",
- "MPI.Scatter!(data, local_data, comm, root=0)\n",
- "\n",
- "# Compute local average\n",
- "local_average = sum(local_data) / length(local_data)\n",
- "\n",
- "# Gather local averages at the root processor\n",
- "gathered_averages = Vector{Float64}(undef, size)\n",
- "MPI.Gather!(local_average, gathered_averages, comm, root=0)\n",
- "\n",
- "if rank == 0\n",
- " # Compute global average of sub-averages\n",
- " global_average = sum(gathered_averages) / size\n",
- " println(\"Global average: $global_average\")\n",
- "end\n",
- "\n",
- "MPI.Finalize()\n",
- "```"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e65cb53f",
- "metadata": {},
- "outputs": [],
- "source": [
- "using MPI\n",
- "using Random\n",
- "\n",
- "# TODO: check if this runs correctly\n",
- "\n",
- "MPI.Init()\n",
- "\n",
- "comm = MPI.COMM_WORLD\n",
- "rank = MPI.Comm_rank(comm)\n",
- "size = MPI.Comm_size(comm)\n",
- "\n",
- "# Root processor generates random data\n",
- "data = rand(rank == 0 ? size * 2 : 0)\n",
- "\n",
- "# Scatter data to all processes\n",
- "local_data = Vector{Float64}(undef, 2)\n",
- "MPI.Scatter!(data, local_data, comm, root=0)\n",
- "\n",
- "# Compute local average\n",
- "local_average = sum(local_data) / length(local_data)\n",
- "\n",
- "# Gather local averages at the root processor\n",
- "gathered_averages = Vector{Float64}(undef, size)\n",
- "MPI.Gather!(local_average, gathered_averages, comm, root=0)\n",
- "\n",
- "if rank == 0\n",
- " # Compute global average of sub-averages\n",
- " global_average = sum(gathered_averages) / size\n",
- " println(\"Global average: $global_average\")\n",
- "end\n",
- "\n",
- "MPI.Finalize()"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "dfd5da9e",
- "metadata": {},
- "source": [
- "In this example, the root processor generates random data and then scatters it to all processes using MPI.Scatter. Each process calculates the average of its local data, and then the local averages are gathered using MPI.Gather. The root processor computes the global average of all sub-averages and prints it."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "5e8f6e6a",
- "metadata": {},
- "source": [
- "# License\n",
- "\n",
- "TODO: replace link to website\n",
- "\n",
- "This notebook is part of the course [Programming Large Scale Parallel Systems](http://localhost:8000/) at Vrije Universiteit Amsterdam and may be used under a [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) license."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c9364808",
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Julia 1.9.1",
- "language": "julia",
- "name": "julia-1.9"
- },
- "language_info": {
- "file_extension": ".jl",
- "mimetype": "application/julia",
- "name": "julia",
- "version": "1.9.1"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/docs/src/notebooks/mpi_tutorial.md b/docs/src/notebooks/mpi_tutorial.md
deleted file mode 100644
index a57c6a2..0000000
--- a/docs/src/notebooks/mpi_tutorial.md
+++ /dev/null
@@ -1,30 +0,0 @@
-```@meta
-EditURL = "https://github.com/fverdugo/XM_40017/blob/main/docs/src/notebooks/mpi_tutorial.ipynb"
-```
-
-```@raw html
-
- Tip
-
-
-
- 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
- .
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-![](\"https://upload.wikimedia.org/wikipedia/commons/thumb/3/39/VU_logo.png/800px-VU_logo.png?20161029201021\"){kind=logo}
\n",
- "\n",
- "### Programming large-scale parallel systems"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e8549215",
- "metadata": {},
- "source": [
- "# Jacobi method"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6e0ef563",
- "metadata": {},
- "source": [
- "## Contents\n",
- "\n",
- "In this notebook, we will learn\n",
- "\n",
- "- How to paralleize a Jacobi method\n",
- "- How the data partition can impact the performance of a distributed algorithm\n",
- "- How to use latency hiding\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 2,
- "id": "1dc78750",
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "jacobi_4_check (generic function with 1 method)"
- ]
- },
- "execution_count": 2,
- "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\")\n",
- "jacobi_3_check(answer) = answer_checker(answer, \"c\")\n",
- "jacobi_4_check(anwswer) = answer_checker(answer, \"d\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d4cb59d5",
- "metadata": {},
- "source": [
- "## The Jacobi method for the Laplace equation\n",
- "\n",
- "\n",
- "The [Jacobi method](https://en.wikipedia.org/wiki/Jacobi_method) is a numerical tool to solve systems of linear algebraic equations. One of the main applications of the method is to solve boundary value problems (BVPs). I.e., given the values at the boundary (of a grid), the Jacobi method will find the interior values that fulfill a certain equation.\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6a2bdbc6",
- "metadata": {},
- "source": [
- "![](\"https://raw.githubusercontent.com/fverdugo/W2bYx4/main/notebooks/figures/logo.png\"){kind=logo}
\n",
- "