Programming large-scale parallel systems

Jacobi method

Contents

In this notebook, we will learn

• How to paralleize a Jacobi method
• How the data partition can impact the performance of a distributed algorithm
• How to use latency hiding

Note: Do not forget to run the next cell before starting studying this notebook.

using Printf

function answer_checker(answer,solution)
    if answer == solution
        "🥳 Well done! "
    else
        "It's not correct. Keep trying! 💪"
    end |> println
end
gauss_seidel_1_check(answer) = answer_checker(answer,"c")
jacobi_1_check(answer) = answer_checker(answer, "d")
jacobi_2_check(answer) = answer_checker(answer, "b")
jacobi_3_check(answer) = answer_checker(answer, "c")
jacobi_4_check(anwswer) = answer_checker(answer, "d")

The Jacobi method for the Laplace equation

The Jacobi method is a numerical tool to solve systems of linear algebraic equations. One of the main applications of the Jacobi method is to solve the equations resulting from boundary value problems (BVPs). I.e., given the values at the boundary (of a grid), we are interested in finding the interior values that fulfill a certain equation.

When solving a Laplace equation in 1D, the Jacobi method leads to the following iterative scheme: The entry $i$ of vector $u$ at iteration $t+1$ is computed as:

$u^{t+1}_i = \dfrac{u^t_{i-1}+u^t_{i+1}}{2}$

Serial implementation

The following code implements the iterative scheme above for boundary conditions -1 and 1 on a grid with $n$ interior points.

function jacobi(n,niters)
    u = zeros(n+2)
    u[1] = -1
    u[end] = 1
    u_new = copy(u)
    for t in 1:niters
        for i in 2:(n+1)
            u_new[i] = 0.5*(u[i-1]+u[i+1])
        end
        u, u_new = u_new, u
    end
    u
end

jacobi(5,0)

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.

Where can we exploit parallelism?

Look at the two nested loops in the sequential implementation:

for t in 1:nsteps
    for i in 2:(n+1)
        u_new[i] = 0.5*(u[i-1]+u[i+1])
    end
    u, u_new = u_new, u
end

• The outer loop cannot be parallelized. The value of u at step t+1 depends on the value at the previous step t.
• The inner loop can be parallelized.

The Gauss-Seidel method

The usage of u_new seems a bit unnecessary at first sight, right?. If we remove it, we get another method called Gauss-Seidel.

function gauss_seidel(n,niters)
    u = zeros(n+2)
    u[1] = -1
    u[end] = 1
    for t in 1:niters
        for i in 2:(n+1)
            u[i] = 0.5*(u[i-1]+u[i+1])
        end
    end
    u
end

Note that the final solution is nearly the same for a large enough number of iterations.

gauss_seidel(5,1000)

Question: Which of the two loops in the Gauss-Seidel method are trivially parallelizable?

for t in 1:niters
    for i in 2:(n+1)
        u[i] = 0.5*(u[i-1]+u[i+1])
    end
end

a) Both of them
b) The outer, but not the inner
c) None of them
d) The inner, but not the outer

answer = "x" # replace x with a, b, c or d
gauss_seidel_1_check(answer)

Parallelization of the Jacobi method

Now, let us parallelize the Jacobi method.

Parallelization strategy

Remember that a sufficiently large grain size is needed to achieve performance in a distributed algorithm. For Jacobi, one could update each entry of vector u_new in a different worker, but this would not be efficient. Instead, we use a parallelization strategy with a larger grain size that is analogous to the algorithm 3 we studied for the matrix-matrix multiplication:

• Each worker updates a consecutive section of the array u_new

The following figure displays the data distribution over 3 processes.

Data dependencies

Recall the Jacobi update:

u_new[i] = 0.5*(u[i-1]+u[i+1])

Thus, in order to update the local entries in u_new, we also need remote entries of vector u located in neighboring processes. Figure below shows the entries of u needed to update the local entries of u_new in a particular process (CPU 2).

Efficiency

• We update $N/P$ entries in each process at each iteration, where $N$ is the total length of the vector and $P$ the number of processes
• We need to get remote entries from 2 neighbors (2 messages per iteration)
• We need to communicate 1 entry per message
• Communication/computation ration is $O(P/N)$ (potentially scalable if $P<<N$)

Extension to 2D

Now, let us study the method for a 2D example.

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:

$u^{t+1}_{(i,j)} = \dfrac{u^t_{(i-1,j)}+u^t_{(i+1,j)}+u^t_{(i,j-1)}+u^t_{(i,j+1)}}{4}$

Note that each entry is updated as the average of the four neighbors (top,bottom,left,right) of that entry in the previous iteration.

Serial implementation

The next code implements a simple example, where the boundary values are equal to 1.

function jacobi_2d(n,niters)
    u = zeros(n+2,n+2)
    u[1,:] = u[end,:] = u[:,1] = u[:,end] .= 1
    u_new = copy(u)
    for t in 1:niters
        for j in 2:(n+1)
            for i in 2:(n+1)
                north = u[i,j+1]
                south = u[i,j-1]
                east = u[i+1,j]
                west = u[i-1,j]
                u_new[i,j] = 0.25*(north+south+east+west)
            end
        end
        u, u_new = u_new, u
    end
    u
end

u = jacobi_2d(10,0)

Where can we exploit parallelism?

for t in 1:niters
    for j in 2:(n+1)
        for i in 2:(n+1)
            north = u[i,j+1]
            south = u[i,j-1]
            east = u[i+1,j]
            west = u[i-1,j]
            u_new[i,j] = 0.25*(north+south+east+west)
        end
    end
    u, u_new = u_new, u
end

• The outer loop cannot be parallelized (like in the 1d case).
• The two inner loops can be parallelized

Parallelization strategies

In 2d one has more flexibility in order to distribute the data over the processes. We consider these three alternatives:

• 1D block partition (each worker handles a subset of consecutive rows and all columns)
• 2D block partition (each worker handles a subset of consecutive rows and columns)
• 2D cyclic partition (each workers handles a subset of alternating rows ans columns)

The three partition types are depicted in the following figure for 4 processes.

Which of the thee alternatives is more efficient? To answer this question we need to quantify how much data is processed and communicated in each case. The following analysis assumes that the grid is of $N$ by $N$ cells and that the number of processes is $P$.

1D block partition

The following figure shows the portion of vector u_new is updated at each iteration by a particular process (CPU 3) left picture, and which entries of u are needed to update this data, right picture. We use analogous figures for the other partitions below.

• We update $N^2/P$ items per iteration
• We need data from 2 neighbors (2 messages per iteration)
• We communicate $N$ items per message
• Communication/computation ratio is $O(P/N)$

2D block partition

• We update $N^2/P$ items per iteration
• We need data from 4 neighbors (4 messages per iteration)
• We communicate $N/\sqrt{P}$ items per message
• Communication/computation ratio is $O(\sqrt{P}/N)$

2D cyclic partition

• We update $N^2/P$ items
• We need data from 4 neighbors (4 messages per iteration)
• We communicate $N^2/P$ items per message (the full data owned by the neighbor)
• Communication/computation ratio is $O(1)$

Which partition is the best one?

• Both 1d and 2d block partitions are potentially scalable if $P<<N$
• The 2d block partition is with the lowest communication complexity
• The 1d block partition requires to send less messages (It can be useful if the fixed cost of sending a message is high)
• The best strategy for a given problem size will thus depend on the machine, but the 2d block partition is the one used in practice since it has the lowest communication complexity.
• Cyclic partitions are impractical for this application (but they are useful in others)

Implementation

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.

Ghost (aka halo) cells

A usual way of implementing the Jacobi method and related algorithms 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.

Thus, the algorithm is usually implemented following two main phases at each iteration Jacobi:

1. Fill the ghost entries with communications
2. Do the Jacobi update sequentially at each process

Code

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.

] add MPI MPIClusterManagers

using MPIClusterManagers 
using Distributed

if procs() == workers()
    nw = 3
    manager = MPIWorkerManager(nw)
    addprocs(manager)
end

# Test cell, remove me
u = [-1, 0, 0, 0, 0, 1]
view(u, 6:6)

@mpi_do manager begin
    using MPI
    comm = MPI.Comm_dup(MPI.COMM_WORLD)
    nw = MPI.Comm_size(comm)
    iw = MPI.Comm_rank(comm)+1
    function jacobi_mpi(n,niters)
        if mod(n,nw) != 0
            println("n must be a multiple of nw")
            MPI.Abort(comm,1)
        end
        n_own = div(n,nw)
        u = zeros(n_own+2)
        u[1] = -1
        u[end] = 1
        u_new = copy(u)
        for t in 1:niters
            reqs = MPI.Request[]
            # Exchange cell values with neighbors
            if iw != 1
                neig_rank = (iw-1)-1
                req = MPI.Isend(view(u,2:2),comm,dest=neig_rank,tag=0)
                push!(reqs,req)
                req = MPI.Irecv!(view(u,1:1),comm,source=neig_rank,tag=0)
                push!(reqs,req)
            end
            if iw != nw
                neig_rank = (iw+1)-1
                s = n_own+1
                r = n_own+2
                req = MPI.Isend(view(u,s:s),comm,dest=neig_rank,tag=0)
                push!(reqs,req)
                req = MPI.Irecv!(view(u,r:r),comm,source=neig_rank,tag=0)
                push!(reqs,req)
            end
            MPI.Waitall(reqs)
            for i in 2:(n_own+1)
                u_new[i] = 0.5*(u[i-1]+u[i+1])
            end
            u, u_new = u_new, u
        end
        u
        @show u
        # Gather results in root process
        results = zeros(n+2)
        results[1] = -1
        results[n+2] = 1
        MPI.Gather!(view(u,2:n_own+1), view(results, 2:n+1), root=0, comm)
        if iw == 1
            @show results
        end            
    end
    niters = 100
    load = 4
    n = load*nw
    jacobi_mpi(n,niters)
end

Question: How many messages per iteration are sent from a process away from the boundary?

a) 1
b) 2
c) 3
d) 4

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. 

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.

answer = "x" # replace x with a, b, c or d.
jacobi_4_check(answer)

Latency hiding

Can our implementation above be improved? 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).

Exercises

Exercise 1

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)

@mpi_do manager begin
    using MPI
    comm = MPI.Comm_dup(MPI.COMM_WORLD)
    nw = MPI.Comm_size(comm)
    iw = MPI.Comm_rank(comm)+1
    function jacobi_mpi(n,niters)
        if mod(n,nw) != 0
            println("n must be a multiple of nw")
            MPI.Abort(comm,1)
        end
        n_own = div(n,nw)
        u = zeros(n_own+2)
        u[1] = -1
        u[end] = 1
        u_new = copy(u)
        for t in 1:niters
            reqs = MPI.Request[]
            # Exchange cell values with neighbors
            if iw != 1
                neig_rank = (iw-1)-1
                req = MPI.Isend(view(u,2:2),comm,dest=neig_rank,tag=0)
                push!(reqs,req)
                req = MPI.Irecv!(view(u,1:1),comm,source=neig_rank,tag=0)
                push!(reqs,req)
            end
            if iw != nw
                neig_rank = (iw+1)-1
                s = n_own+1
                r = n_own+2
                req = MPI.Isend(view(u,s:s),comm,dest=neig_rank,tag=0)
                push!(reqs,req)
                req = MPI.Irecv!(view(u,r:r),comm,source=neig_rank,tag=0)
                push!(reqs,req)
            end
            MPI.Waitall(reqs)
            for i in 2:(n_own+1)
                u_new[i] = 0.5*(u[i-1]+u[i+1])
            end
            u, u_new = u_new, u
        end
        u
        @show u
        # Gather results in root process
        results = zeros(n+2)
        results[1] = -1
        results[n+2] = 1
        MPI.Gather!(view(u,2:n_own+1), view(results, 2:n+1), root=0, comm)
        if iw == 1
            @show results
        end            
    end
    niters = 100
    load = 4
    n = load*nw
    jacobi_mpi(n,niters)
end

License

TODO: replace link to website

This notebook is part of the course Programming Large Scale Parallel Systems at Vrije Universiteit Amsterdam and may be used under a CC BY 4.0 license.