# mpi_tutorial

### Programming large-scale parallel systems¶

# Distributed computing with MPI¶

## Contents¶

In this notebook, we will learn the basics of parallel computing using the Message Passing Interface (MPI) from Julia. In particular, we will learn:

- How to run parallel MPI code in Julia
- How to use basic collective communication directives
- How to use basic point-to-point communication directives

For further information on how to use MPI from Julia see https://github.com/JuliaParallel/MPI.jl

## What is MPI ?¶

- MPI stands for the "Message Passing Interface"
- It is a standardized library specification for communication between parallel processes in distributed-memory systems.
- It is the gold-standard for distributed computing in HPC systems since the 90s
- It is huge: the MPI standard has more than 1k pages (see https://www.mpi-forum.org/docs/mpi-4.0/mpi40-report.pdf)
- There are several implementations of this standard (OpenMPI, MPICH, IntelMPI)
- The interface is in C and FORTRAN (C++ was deprecated)
- There are Julia bindings via the package MPI.jl https://github.com/JuliaParallel/MPI.jl

### What is MPI.jl ?¶

- It is not a Julia implementation of the MPI standard
- It is a wrapper to the C interface of MPI
- You need a C MPI installation in your system

MPI.jl provides a convenient Julia API to access MPI. For instance, this is how you get the id (rank) of the current process.

    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)

Internally, MPI.jl uses `ccall` which is a mechanism that allows you to call C functions from Julia. In this, example we are calling the C function `MPI_Comm_rank` from the underlying MPI installation.

    comm = MPI.COMM_WORLD
    rank_ref = Ref{Cint}()
    ccall((:MPI_Comm_rank, MPI.API.libmpi), Cint, (MPI.API.MPI_Comm, Ptr{Cint}), comm, rank_ref)
    rank = Int(rank_ref[])

### Installing MPI in Julia¶

The Jupyter Julia kernel installed by IJulia activates the folder where the notebook is located as the default environment, which causes the main process and the worker processes to not share the same environment. Therefore, we need to set the environment as the global environment.

In [ ]:

    ] activate

MPI can be installed as any other Julia package using the package manager.

In [ ]:

    ] add MPI

Note: 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. See more information in https://github.com/JuliaParallel/MPI.jl

### Basic information about MPI processes¶

The following cells give information about MPI processes, such as the rank id, the total number of processes and the name of the host running the code respectively. Before calling this functions one needs to initialize MPI.

In [ ]:

    using MPI
    MPI.Init()

In [ ]:

    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)

In [ ]:

    nranks = MPI.Comm_size(comm)

In [ ]:

    MPI.Get_processor_name()

Note that this note notebook is not running with different MPI processes (yet). So using MPI will only make sense later when we add more processes.

### Hello-world example¶

In [ ]:

    using MPI
    MPI.Init()
    comm = MPI.COMM_WORLD
    nranks = MPI.Comm_size(comm)
    rank = MPI.Comm_rank(comm)
    println("Hello, I am process $rank of $nranks processes!")

### Hello world in C¶

    #include <mpi.h>
    #include <stdio.h>
    int main(int argc, char** argv) {
        MPI_Init(NULL, NULL);
        int world_size;
        MPI_Comm_size(MPI_COMM_WORLD, &world_size);
        int world_rank;
        MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
        char processor_name[MPI_MAX_PROCESSOR_NAME];
        int name_len;
        MPI_Get_processor_name(processor_name, &name_len);
        printf("Hello from %s, I am rank %d of %d ranks!\n",
               processor_name, world_rank, world_size);
        MPI_Finalize();
    }

## Running MPI code¶

### Creating MPI processes (aka ranks)¶

- MPI processes are created with the driver program `mpiexec`
- `mpiexec` takes an application and runs it on different ranks
- The application calls MPI directives to communicate between these ranks
- The application can be Julia running your script in particular.

### Execution model¶

- MPI programs are typically run with a Single Program Multiple Data (SPMD) model
- But the standard supports Multiple Program Multiple Data (MPMD)

### Hello world¶

Julia code typically needs to be in a file to run it in with MPI. Let's us write our hello world in a file:

In [ ]:

    code = raw"""
    using MPI
    MPI.Init()
    comm = MPI.COMM_WORLD
    nranks = MPI.Comm_size(comm)
    rank = MPI.Comm_rank(comm)
    println("Hello, I am process $rank of $nranks processes!")
    """
    filename = tempname()*".jl"
    write(filename,code);

Now, we can run it

In [ ]:

    using MPI
    mpiexec(cmd->run(`$cmd -np 4 julia --project=. $filename`));

Note that function mpiexec provided by MPI.jl is a convenient way of accessing the mpiexec program that matches the MPI installation used my Julia.

### MPIClusterManagers¶

- This package allows you to create Julia workers that can call MPI functions
- This is useful to combine Distributed.jl and MPI.jl
- E.g., we can run MPI code interactively (from a notebook)
- Link: https://github.com/JuliaParallel/MPIClusterManagers.jl

In [ ]:

    ] add MPIClusterManagers

In [ ]:

    using MPIClusterManagers
    using Distributed
    if procs() == workers()
        nranks = 3
        manager = MPIWorkerManager(nranks)
        addprocs(manager)
    end

In [ ]:

    @everywhere workers() begin
        using MPI
        MPI.Init()
        comm = MPI.COMM_WORLD
        nranks = MPI.Comm_size(comm)
        rank = MPI.Comm_rank(comm)
        println("Hello, I am process $rank of $nranks processes!")
    end

Note: Note that the rank ids start with 0.

## MPI Communicators¶

In MPI, a communicator represents a group of processes that can communicate with each other. `MPI_COMM_WORLD` (`MPI.COMM_WORLD` from Julia) is a built-in communicator that represents all processes available 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 (integer) identifier assigned to each process within a communicator. It allows processes to distinguish and address each other in communication operations.

### Duplicating a communicator¶

It is a good practice to not using the built-in communicators directly, and use a copy instead with `MPI.Comm_dup`. Different libraries using the same communicator can lead to unexpected interferences.

## Collective communication¶

MPI provides collective communication functions for communication involving multiple processes. Some usual collective directives are:

- `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.

See more collective directives available from Julia here: https://juliaparallel.org/MPI.jl/stable/reference/collective/

### Gather¶

Each rank sends a message to the root rank (the root rank also sends a message to itself). The root rank receives all these values in a buffer (e.g. a vector).

In [ ]:

    @everywhere workers() begin
        comm = MPI.Comm_dup(MPI.COMM_WORLD)
        nranks = MPI.Comm_size(comm)
        rank = MPI.Comm_rank(comm)
        root = 0
        snd = 10*(rank+2)
        println("I am sending $snd")
        rcv = MPI.Gather(snd,comm;root)
        if rank == root
            println("I have received: $rcv")
        end
    end

### Scatter¶

The root rank contains a buffer (e.g., a vector) of values (one value for each rank in a communicator). Scatter sends one value to each rank (the root rank also receives a value). The root rank can be any process in a communicator.

In [ ]:

    @everywhere workers() begin
        comm = MPI.Comm_dup(MPI.COMM_WORLD)
        nranks = MPI.Comm_size(comm)
        rank = MPI.Comm_rank(comm)
        root = 0
        rcv = Ref(0)
        if rank == root
            snd = [10*(i+1) for i in 1:nranks]
            println("I am sending: $snd")
        else
            snd = nothing
        end
        MPI.Scatter!(snd,rcv,comm;root)
        println("I have received: $(rcv[])")
    end

### Bcast (broadcast)¶

Similar to scatter, but we send the same message to all processes.

In [ ]:

    @everywhere workers() begin
        comm = MPI.Comm_dup(MPI.COMM_WORLD)
        nranks = MPI.Comm_size(comm)
        rank = MPI.Comm_rank(comm)
        root = 0
        buffer = Ref(0)
        if rank == root
            buffer[] = 20
            println("I am sending: $(buffer[])")
        end
        MPI.Bcast!(buffer,comm;root)
        println("I have received: $(buffer[])")
    end

### Barrier¶

In [ ]:

    @everywhere workers() begin
        comm = MPI.Comm_dup(MPI.COMM_WORLD)
        sleep(rand(1:3))
        MPI.Barrier(comm)
        println("Done!")
    end

## Point-to-Point communication¶

MPI also provides point-to-point communication directives for arbitrary communication between processes. Point-to-point communications are two-sided: there is a sender and a receiver. Here, we will discuss these basic directives:

- `MPI.Isend`, and `MPI.Irecv!` (non-blocking directives)
- `MPI.Send`, and `MPI.Recv!` (blocking directives)

Non-blocking directives return immediately and return an `MPI.Request` object. This request object can be queried with functions like `MPI.Wait`. It is mandatory to wait on the request object before reading the receive buffer, or before writing again on the send buffer.

For blocking directives, it is save to read/write from/to the receive/send buffer once the function has returned. By default, blocking directives might wait (or might not wait) for a matching send/receive. For fine control, MPI offers advanced blocking directives with different blocking behaviors (called communication modes, see section 3.9 of the MPI standard 4.0). Blocking communication will be discussed later in the course.

### Blocking communication¶

If we start a receive before a matching send, we will block in the call to `MPI.Recv!`. Run the next cell and note that the message is not printed since the process is blocked at `MPI.Recv!` waiting for a matching send.

In [ ]:

    @spawnat 4 begin
        buffer = Ref(0)
        comm = MPI.COMM_WORLD
        MPI.Recv!(buffer,comm;source=2-2,tag=0)
        println("I have received $(buffer[]).")
    end;

If you run the next cell containing the corresponding send, the communication will take place.

In [ ]:

    @spawnat 2 begin
        buffer = Ref(2)
        comm = MPI.COMM_WORLD
        MPI.Send(buffer,comm;dest=4-2,tag=0)
        println("I have send $(buffer[]). It is now safe to overwite the buffer.")
    end;

### MPI does not integrate well with Julia Tasks¶

MPI blocks without yielding (we cannot switch to other Julia tasks). Run next cell:

In [ ]:

    @spawnat 4 begin
        buffer = Ref(0)
        comm = MPI.COMM_WORLD
        MPI.Recv!(buffer,comm;source=2-2,tag=0)
        println("I have received $(buffer[]).")
    end;

Now try to spawn other tasks on process 4 by running next cell. This task will not be served yet.

In [ ]:

    @spawnat 4 println("Hello!");

We first need to unlock the receive with a matching send. Then the task printing "Hello!" will be finally served.

In [ ]:

    @spawnat 2 begin
        buffer = Ref(2)
        comm = MPI.COMM_WORLD
        MPI.Send(buffer,comm;dest=4-2,tag=0)
        println("I have send $(buffer[]). It is now safe to overwite the buffer.")
    end;

### Non-blocking communication¶

In [ ]:

    @spawnat 4 begin
        buffer = Ref(0)
        comm = MPI.COMM_WORLD
        req = MPI.Irecv!(buffer,comm;source=2-2,tag=0)
        println("Not yet safe to read the buffer")
        MPI.Wait(req)
        println("I have received $(buffer[]).")
    end;

In [ ]:

    @spawnat 2 begin
        buffer = Ref(2)
        comm = MPI.COMM_WORLD
        req = MPI.Isend(buffer,comm;dest=4-2,tag=0)
        println("Not yet safe to write the buffer")
        MPI.Wait(req)
        println("I have send $(buffer[]). It is now safe to overwite the buffer.")
    end;

### Combining Julia Tasks and non-blocking communication¶

We can implement our own blocking mechanism that combines well with the Julia task scheduler. This can be done by testing on the request in a while loop. The `yield` in the loop provides the opportunity of other tasks to run.

In [ ]:

    @spawnat 4 begin
        buffer = Ref(0)
        comm = MPI.COMM_WORLD
        req = MPI.Irecv!(buffer,comm;source=2-2,tag=0)
        while MPI.Test(req) == false
            yield()
        end
        println("I have received $(buffer[]).")
    end;

In [ ]:

    @spawnat 4 println("Hello!");

In [ ]:

    @spawnat 2 begin
        buffer = Ref(2)
        comm = MPI.COMM_WORLD
        req = MPI.Isend(buffer,comm;dest=4-2,tag=0)
        while MPI.Test(req) == false
            yield()
        end
        println("I have send $(buffer[]). It is now safe to overwite the buffer.")
    end;

### Example¶

The first rank generates a message and sends it to the last rank. The last rank receives the message and multiplies it by a coefficient. The last rank sends the result back to the first rank.

In [ ]:

    @everywhere workers() begin
        comm = MPI.Comm_dup(MPI.COMM_WORLD)
        rank = MPI.Comm_rank(comm)
        nranks = MPI.Comm_size(comm)
        snder = 0
        rcver = nranks-1
        buffer = Ref(0)
        if rank == snder
            msg = 10*(rank+2)
            println("I am sending: $msg")
            buffer[] = msg
            req = MPI.Isend(buffer,comm;dest=rcver,tag=0)
            MPI.Wait(req)
            req = MPI.Irecv!(buffer,comm,source=rcver,tag=0)
            MPI.Wait(req)
            msg = buffer[]
            println("I have received: $msg")
        end
        if rank == rcver
            req = MPI.Irecv!(buffer,comm,source=snder,tag=0)
            MPI.Wait(req)
            msg = buffer[]
            println("I have received: $msg")
            coef = (rank+2)
            msg = msg*coef
            println("I am sending: $msg")
            buffer[] = msg
            req = MPI.Isend(buffer,comm;dest=snder,tag=0)
            MPI.Wait(req)
        end
    end

Important: In non-blocking communication, use `MPI.Wait()` before modifying the send buffer or using the receive buffer.

### Example (with blocking directives)¶

In [ ]:

    @everywhere workers() begin
        comm = MPI.Comm_dup(MPI.COMM_WORLD)
        rank = MPI.Comm_rank(comm)
        nranks = MPI.Comm_size(comm)
        snder = 0
        rcver = nranks-1
        buffer = Ref(0)
        if rank == snder
            msg = 10*(rank+2)
            println("I am sending: $msg")
            buffer[] = msg
            MPI.Send(buffer,comm;dest=rcver,tag=0)
            MPI.Recv!(buffer,comm,source=rcver,tag=0)
            msg = buffer[]
            println("I have received: $msg")
        end
        if rank == rcver
            MPI.Recv!(buffer,comm,source=snder,tag=0)
            msg = buffer[]
            println("I have received: $msg")
            coef = (rank+2)
            msg = msg*coef
            println("I am sending: $msg")
            buffer[] = msg
            MPI.Send(buffer,comm;dest=snder,tag=0)
        end
    end

Important: Blocking directives might look simpler to use, but they can lead to dead locks if the sends and receives are not issued in the right order. Non-blocking directives can also lead to dead locks, but when waiting for the request, not when calling the send/receive functions.

## Exercises¶

### Exercise 1¶

Implement this simple algorithm: Rank 0 generates a message (an integer). Rank 0 sends the message to rank 1. Rank 1 receives the message, increments the message by 1, and sends the result to rank 2. Rank 2 receives the message, increments the message by 1, and sends the result to rank 3. Etc. The last rank sends back the message to rank 0 closing the ring. See the next figure. Implement the communications using MPI.

### Exercise 2¶

Implement the same algorithm as in Exercise 1, but now without using MPI. Implement the communications using the native `Distributed` module provided by Julia. In this case, start using process 1 instead of rank 0.

# License¶

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.