{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7606d30a",
   "metadata": {},
   "source": [
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/3/39/VU_logo.png/800px-VU_logo.png?20161029201021\" width=\"350\">\n",
    "\n",
    "### Programming large-scale parallel systems"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ac1e5d9",
   "metadata": {},
   "source": [
    "# Intro to MPI (collectives)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a341be2e",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "\n",
    "In this notebook, we will continue learning the basics of parallel computing using the Message Passing Interface (MPI) from Julia. In particular, we will learn:\n",
    "\n",
    "- How to use basic collective communication directives\n",
    "- How to use MPI communicators\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c87b3c82",
   "metadata": {},
   "source": [
    "## Collective communication\n",
    "\n",
    "MPI provides a set of routines for communication involving multiple processes. These are called *collective communication* operations. Some usual collective operations are:\n",
    "\n",
    "\n",
    "- `MPI_Barrier`: Synchronize all processes\n",
    "- `MPI_Bcase`: Send the same data from one to all processes\n",
    "- `MPI_Gather`: Gather data from all processes to one\n",
    "- `MPI_Scatter`: Scatter data from one to all all processes\n",
    "- `MPI_Reduce`: Reduce data from all processes to a single one\n",
    "- `MPI_Scan`: Scan (or prefix) reduction\n",
    "- `MPI_Allgather`: Like a `MPI_Gather` but all processes receive the result\n",
    "- `MPI_Allreduce`: Like `MPI_Reduce` but all processes receive the result\n",
    "- `MPI_Alltoall`:  Exchange data from all to all processes\n",
    "\n",
    "\n",
    "We will discuss some of them in greater detail in this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ffa5e56",
   "metadata": {},
   "source": [
    "## Why collective operations?\n",
    "\n",
    "Point-to-point communication functions provide all the building blocks needed in parallel programs and could be used to implement the collective functions described above. Then, why does MPI provide collective communication functions? There are several reasons:\n",
    "\n",
    "- Ease of use: It is handy for users to have these functions readily available instead of having to implement them.\n",
    "- Performance: Library implementations typically use efficient algorithms (such as reduction trees).\n",
    "- Hardware support: Hardware vendors can optimize the MPI library for their machine and even develop hardware specially design to perform MPI operations efficiently. For instance, [IMB Blue Gene](https://en.wikipedia.org/wiki/IBM_Blue_Gene) was equipped with  auxiliary networks for MPI collectives.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53a0f5fa",
   "metadata": {},
   "source": [
    "## Semantics of collective operations\n",
    "\n",
    "These are key properties of collective operations:\n",
    "\n",
    "\n",
    "- Completeness: All the collective communication directives above are *complete* operations. Thus, it is safe to use and reset the buffers once the function returns.\n",
    "- Standard mode: Collective directives are in standard mode only, like `MPI_Send`. Assuming that they block is erroneous, assuming that they do not block is also erroneous.\n",
    "- Synchronization:  Completion of a call does not guarantee that other processes have completed the operation. A collective operation may or may not have the effect of synchronizing all processes, the only exception is `MPI_Barrier` of course.\n",
    "\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b> Recent versions of the MPI standard also include non-blocking (i.e., incomplete) versions of collective operations (not covered in this notebook). A particularly funny one is the non-blocking barrier `MPI_Ibarrier`.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9e2a219",
   "metadata": {},
   "source": [
    "## MPI_Barrier\n",
    "\n",
    "This function is used to synchronizes a group of processes. All processes block until all have reached the barrier. It is often invoked at the end of for loops to make sure all processes have finished the current loop iteration to move to the next one. We will see a practical example later in another notebook when studying the traveling sales person problem (TSP).\n",
    "\n",
    "In Julia:\n",
    "```julia\n",
    "MPI.Barrier(comm)\n",
    "```\n",
    "\n",
    "In C:\n",
    "```c\n",
    "int MPI_Barrier(MPI_Comm comm)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9753e4b",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "In this example the ranks sleep for a random amount of time and then they call barrier. It is guaranteed that the message \"Done!\" will be printed after all processes printed \"I woke up\" since we used a barrier. Try also to comment out the call to `MPI.Barrier`. You will see that the message can be printed in any order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6e78389",
   "metadata": {},
   "outputs": [],
   "source": [
    "using MPI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a66c717",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    sleep(rand(1:3))\n",
    "    println(\"[rank $rank] I woke up\")\n",
    "    MPI.Barrier(comm)\n",
    "    println(\"[rank $rank] Done!\")\n",
    "end\n",
    "run(`$(mpiexec()) -np 4 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df0223ff",
   "metadata": {},
   "source": [
    "## MPI_Reduce\n",
    "\n",
    "This function combines values provided by different processors according to a given reduction operation. The result is received in a single process (called the root process).  The root process can be any process and it is rank 0 by default in Julia.\n",
    "\n",
    "In Julia:\n",
    "```julia\n",
    "MPI.Reduce!(sendbuf, recvbuf, op, comm; root)\n",
    "```\n",
    "\n",
    "In C:\n",
    "```\n",
    "int MPI_Reduce(const void *sendbuf, void *recvbuf, int count,\n",
    "               MPI_Datatype datatype, MPI_Op op, int root,\n",
    "               MPI_Comm comm)\n",
    "```\n",
    "\n",
    "\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b> `MPI_Reduce` (and related functions like `MPI_Allreduce`) use operations that are *commutative* and *associative*, such as addition, multiplication, maximum, minimum, logical and, logical or. Otherwise the result will be undefined as it will depend in which order the operations are done.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fc1aa28",
   "metadata": {},
   "source": [
    "### Example\n",
    "Each process generates a random integer. Compute the sum of all these integers. Receive the result on rank 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "341b8438",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    value = rand(1:4)\n",
    "    root = 0\n",
    "    sndbuf = [value]\n",
    "    if rank == root\n",
    "        rcvbuf = [0]\n",
    "    else\n",
    "        # This can be anything as only the root process\n",
    "        # Will receive data\n",
    "        rcvbuf = nothing\n",
    "    end\n",
    "    println(\"[rank $rank] I am contributing with $sndbuf\")\n",
    "    MPI.Reduce!(sndbuf, rcvbuf, +, comm; root)\n",
    "    println(\"[rank $rank] The sum is $rcvbuf\")\n",
    "end\n",
    "run(`$(mpiexec()) -np 3 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ea77a8f",
   "metadata": {},
   "source": [
    "## Reducing multiple values\n",
    "\n",
    "If more than one item is provided in the send buffers, `MPI_Reduce` does the reduction element by element.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b26fd0b2",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Each process generates 2 random integers. Take a look to the result of reduce. Note that the values have ben reduced \"element by element\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a145948",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    sndbuf = rand(1:4,2)\n",
    "    root = 0\n",
    "    if rank == root\n",
    "        rcvbuf = [0,0]\n",
    "    else\n",
    "        rcvbuf = nothing\n",
    "    end\n",
    "    println(\"[rank $rank] I am contributing with $sndbuf\")\n",
    "    MPI.Reduce!(sndbuf, rcvbuf, +, comm; root)\n",
    "    println(\"[rank $rank] The sum is $rcvbuf\")\n",
    "end\n",
    "run(`$(mpiexec()) -np 3 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "539df34f",
   "metadata": {},
   "source": [
    "## MPI_Allreduce\n",
    "\n",
    "This is a variant of `MPI_Reduce` in which all processes get the result."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be0381bb",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "This is the same example as above, but all processes get the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1595f4f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    sndbuf = rand(1:4,2)\n",
    "    # All ranks allocate a receive buffer\n",
    "    rcvbuf = [0,0]\n",
    "    println(\"[rank $rank] I am contributing with $sndbuf\")\n",
    "    MPI.Allreduce!(sndbuf, rcvbuf, +, comm)\n",
    "    println(\"[rank $rank] The sum is $rcvbuf\")\n",
    "end\n",
    "run(`$(mpiexec()) -np 3 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97dc2886",
   "metadata": {},
   "source": [
    "## MPI_Gather\n",
    "\n",
    "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).\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b>  This function assumes that the amount of data sent from each rank is the same. See `MPI_Gatherv` below for more general cases.\n",
    "</div>\n",
    "\n",
    "\n",
    "In Julia:\n",
    "```julia\n",
    "MPI.Gather!(sendbuf, recvbuf, comm; root=0)\n",
    "```\n",
    "\n",
    "In C:\n",
    "```c\n",
    "int MPI_Gather(const void *sendbuf, int sendcount, MPI_Datatype sendtype,\n",
    "    void *recvbuf, int recvcount, MPI_Datatype recvtype, int root,\n",
    "    MPI_Comm comm)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d0a1e9d",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Each process sends an integer to rank 0.\n"
   ]
  },
  {
   "attachments": {
    "g13884.png": {
     "image/png": "[encoded data, 111604 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "4bdf9c02",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 111604 characters]\"  align=\"left\" width=\"350\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f8a70c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    msg = 10*(rank+2)\n",
    "    root = 0\n",
    "    sndbuf = [msg]\n",
    "    if rank == root\n",
    "        nranks = MPI.Comm_size(comm)\n",
    "        rcvbuf = zeros(Int,nranks)\n",
    "    else\n",
    "        # This can be anything as only the root process\n",
    "        # Will receive data\n",
    "        rcvbuf = nothing\n",
    "    end\n",
    "    println(\"[rank $rank] I am sending $sndbuf\")\n",
    "    MPI.Gather!(sndbuf, rcvbuf, comm; root)\n",
    "    println(\"[rank $rank] I received $rcvbuf\")\n",
    "end\n",
    "run(`$(mpiexec()) -np 3 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f85f08d",
   "metadata": {},
   "source": [
    "## MPI_Allgather\n",
    "\n",
    "This is a variant of `MPI_Gather` in which all processes get the result."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bf2ef34",
   "metadata": {},
   "source": [
    "### Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "070e1fc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    msg = 10*(rank+2)\n",
    "    sndbuf = [msg]\n",
    "    nranks = MPI.Comm_size(comm)\n",
    "    rcvbuf = zeros(Int,nranks)\n",
    "    println(\"[rank $rank] I am sending $sndbuf\")\n",
    "    MPI.Allgather!(sndbuf, rcvbuf, comm)\n",
    "    # All ranks receive the result\n",
    "    println(\"[rank $rank] I received $rcvbuf\")\n",
    "end\n",
    "run(`$(mpiexec()) -np 3 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb6a9034",
   "metadata": {},
   "source": [
    "## MPI_Gatherv\n",
    "\n",
    "This function is similar to `MPI_Gather`, but it is used when there is a different amount of data sent from each process. There is also `MPI_Allgatherv` for which all processes will receive the result."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0d2998c",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Each processes sends `rank+1` integers to rank 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23814278",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    sndbuf = fill(10*rank,rank+1)\n",
    "    root = 0\n",
    "    if rank == root\n",
    "        nranks = MPI.Comm_size(comm)\n",
    "        # Define the amount of data to be received\n",
    "        # from each rank\n",
    "        counts = [ rank+1 for rank in 0:(nranks-1)]\n",
    "        # Total amount of receive data\n",
    "        ndata = sum(counts)\n",
    "        data = zeros(Int,ndata)\n",
    "        rcvbuf = MPI.VBuffer(data,counts)\n",
    "    else\n",
    "        # This can be anything as only the root process\n",
    "        # Will receive data\n",
    "        rcvbuf = nothing\n",
    "    end\n",
    "    println(\"[rank $rank] I am sending $sndbuf\")\n",
    "    MPI.Gatherv!(sndbuf, rcvbuf, comm; root)\n",
    "    if rank == root\n",
    "        println(\"[rank $rank] I received $data\")\n",
    "    end\n",
    "end\n",
    "run(`$(mpiexec()) -np 4 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25c2fb1e",
   "metadata": {},
   "source": [
    "## How to get the amount of data to be received?\n",
    "\n",
    "Remember that `MPI_Probe` allows you to query the amount of data to be received with `MPI_Recv`. Unfortunately, this mechanism does not work for MPI collectives. A workaround is to first communicate the amount of data to be gathered, allocate the receive buffer accordingly, and then communicate the actual message.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12f9e2ca",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Each process sends a random amount of integers to rank 0. The root process will not know the amount of data to be gathered from each rank in advance. We need an auxiliary gather to inform about the message size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b10430a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    # Get number of ints to send\n",
    "    count = rand(1:4)\n",
    "    # Generate data to send\n",
    "    sndbuf = fill(rank,count)\n",
    "    # Communicate counts\n",
    "    root = 0\n",
    "    if rank == root\n",
    "        nranks = MPI.Comm_size(comm)\n",
    "        counts = zeros(Int,nranks)\n",
    "    else\n",
    "        counts = nothing\n",
    "    end\n",
    "    MPI.Gather!([count], counts, comm; root)\n",
    "    # Allocate receive buffer\n",
    "    if rank == root\n",
    "        ndata = sum(counts)\n",
    "        data = zeros(Int,ndata)\n",
    "        rcvbuf = MPI.VBuffer(data,counts)\n",
    "    else\n",
    "        rcvbuf = nothing\n",
    "    end\n",
    "    # Communicate data\n",
    "    println(\"[rank $rank] I am sending $sndbuf\")\n",
    "    MPI.Gatherv!(sndbuf, rcvbuf, comm; root)\n",
    "    if rank == root\n",
    "        println(\"[rank $rank] I received $data\")\n",
    "    end\n",
    "end\n",
    "run(`$(mpiexec()) -np 4 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b4254d1",
   "metadata": {},
   "source": [
    "## MPI_Scatter\n",
    "\n",
    "This is the inverse operation of `MPI_Gather`. The root rank sends a distinct message to each rank.\n",
    "\n",
    "In Julia:\n",
    "```julia\n",
    "MPI.Scatter!(sendbuf, recvbuf, comm; root=0)\n",
    "```\n",
    "\n",
    "In C:\n",
    "```c\n",
    "int MPI_Scatter(const void *sendbuf, int sendcount, MPI_Datatype sendtype,\n",
    "    void *recvbuf, int recvcount, MPI_Datatype recvtype, int root,\n",
    "    MPI_Comm comm)\n",
    "```\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b> `MPI_Scatter` sends the same amount of data to each rank. To send a variable amount of data there is `MPI_Scatterv` (`MPI.Scatterv!` in julia), which works similar to `MPI_Gatherv`. \n",
    "</div>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e55aa353",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Rank 0 sends a unique integer to each rank. "
   ]
  },
  {
   "attachments": {
    "g13389.png": {
     "image/png": "[encoded data, 111664 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "74d5b606",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 111664 characters]\"  align=\"left\" width=\"350\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ea413fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    nranks = MPI.Comm_size(comm)\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    root = 0\n",
    "    rcv = Ref(0) \n",
    "    if rank == root\n",
    "        snd = [10*(i+1) for i in 1:nranks]\n",
    "        println(\"[rank $rank] I am sending: $snd\")\n",
    "    else\n",
    "        # Only the root rank sends data\n",
    "        # so this can be anything\n",
    "        snd = nothing\n",
    "    end    \n",
    "    MPI.Scatter!(snd,rcv,comm;root)\n",
    "    println(\"[rank $rank] I have received: $(rcv[])\")\n",
    "end\n",
    "run(`$(mpiexec()) -np 4 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed8da7f9",
   "metadata": {},
   "source": [
    "## MPI_Bcast\n",
    "\n",
    "Similar to `MPI_Scatter`, but we send the same message to all processes.\n",
    "\n",
    "In Julia:\n",
    "```julia\n",
    "MPI.Bcast!(buf, comm; root)\n",
    "```\n",
    "\n",
    "In C:\n",
    "```c\n",
    "int MPI_Bcast(void *buffer, int count, MPI_Datatype datatype,\n",
    "    int root, MPI_Comm comm)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97369713",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Rank 0 sends the same integer to all processes."
   ]
  },
  {
   "attachments": {
    "g1657.png": {
     "image/png": "[encoded data, 96496 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "26b56c81",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 96496 characters]\"  align=\"left\" width=\"350\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4de15781",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    nranks = MPI.Comm_size(comm)\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    root = 0\n",
    "    buffer = Ref(0)\n",
    "    if rank == root\n",
    "        buffer[] = 20\n",
    "        println(\"[rank $rank] I am sending: $(buffer[])\")\n",
    "    end    \n",
    "    MPI.Bcast!(buffer,comm;root)\n",
    "    println(\"[rank $rank] I have received: $(buffer[])\")\n",
    "end\n",
    "run(`$(mpiexec()) -np 4 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e83be96d",
   "metadata": {},
   "source": [
    "## Communicators\n",
    "\n",
    "A key component in MPI programming are so-called communicators. A communicator object has to main purposes:\n",
    "\n",
    "\n",
    "1. To provide an isolated communication context\n",
    "2. To define a group of processes\n",
    "\n",
    "\n",
    "### Communication context\n",
    "\n",
    "When you write an MPI program it is very likely that you are going to use libraries that also use MPI to send messages. Ideally, these libraries should not interfere with application messages. Using tags to isolate the messages send by your application does not solve the problem. MPI communicators fix this problem as they provided an isolated communication context. For instance, `MPI_SEND` and `MPI_RECV` specify a communicator. `MPI_RECV` can only receive messages sent to same communicator. The same is also true for collective communication directives. If two libraries use different communicators, their message will never interfere. In particular it is recommended to never use the default communicator, `MPI_COMM_WORLD`, directly when working with other libraries. A new isolated communicator can be created with `MPI_Comm_dup`.\n",
    "\n",
    "\n",
    "### Groups of processes\n",
    "\n",
    "On the other hand,  imagine that we want to use an MPI communication directive like `MPI_Gather`, but we only want a subset of the processes to participate in the operation. So far, we have used always the default communication `MPI_COMM_WORLD`, which represents all processes. Thus, by using this communicator, we are including all processes in the operation. We can create other communicators that contain only a subset of processes. To this end, we can use function `MPI_Comm_split`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "556ca223",
   "metadata": {},
   "source": [
    "## MPI_Comm_dup\n",
    "\n",
    "Duplicate a given communicator, creating a new one with an new communication context.\n",
    "\n",
    "\n",
    "In Julia:\n",
    "```julia\n",
    "newcomm = MPI.Comm_dup(comm)\n",
    "```\n",
    "\n",
    "In C:\n",
    "```c\n",
    "int MPI_Comm_dup(MPI_Comm comm, MPI_Comm *newcomm)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9acfb792",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Send value 2 from rank 2 to rank 3. Then send value 1 also from rank 2 to rank 3. Note that we are using 2 different communicators for each send receive pair in such a way that we first receive value 1 and then value 2, even if this not was the other in which the data was sent. Try also to define `comm2 = comm` and see how the result will change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e50d44f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    # We create a new communicator\n",
    "    comm2 = MPI.Comm_dup(comm)\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    if rank == 2\n",
    "        sndbuf = [2]\n",
    "        MPI.Send(sndbuf, comm2; dest=3, tag=0)\n",
    "        sndbuf = [1]\n",
    "        MPI.Send(sndbuf, comm; dest=3, tag=0)\n",
    "    end\n",
    "    if rank == 3\n",
    "        rcvbuf = zeros(Int,1)\n",
    "        MPI.Recv!(rcvbuf, comm; source=2, tag=0)\n",
    "        @show rcvbuf\n",
    "        MPI.Recv!(rcvbuf, comm2; source=2, tag=0)\n",
    "        @show rcvbuf\n",
    "    end\n",
    "end\n",
    "run(`$(mpiexec()) -np 4 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e1b13eb",
   "metadata": {},
   "source": [
    "## MPI_Comm_split\n",
    "\n",
    "Split a given communicator.\n",
    "\n",
    "\n",
    "In Julia:\n",
    "\n",
    "```julia\n",
    "newcomm = MPI.Comm_split(comm, color, key)\n",
    "```\n",
    "\n",
    "In C:\n",
    "```c\n",
    "int MPI_Comm_split(MPI_Comm comm, int color, int key,\n",
    "    MPI_Comm *newcomm)\n",
    "```\n",
    "\n",
    "There are two key parameters:\n",
    "\n",
    "- `color`: all processes with the same color will be grouped in the same new communicator.\n",
    "- `key`: The processes will be ranked in the new communicator according to key, breaking ties with the rank in the old communicator. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b128840",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Create two groups of processes, one for processes with rank id even, and another for the rest. Note that the rank id will be function of the given communicator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b18cb74d",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    color = mod(rank,2)\n",
    "    key = 0\n",
    "    newcomm = MPI.Comm_split(comm, color, key)\n",
    "    newrank = MPI.Comm_rank(newcomm)\n",
    "    println(\"Hello, I am rank $rank in world and rank $newrank in group $color.\")\n",
    "end\n",
    "run(`$(mpiexec()) -np 4 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b22e4bf4",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Processes with odd rank id generate a random integer. Sum all these integers and send them to the first process with odd rank id. Note that we are using a new communicator to perform the reduction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c0e504b",
   "metadata": {},
   "outputs": [],
   "source": [
    "code = quote\n",
    "    using MPI\n",
    "    MPI.Init()\n",
    "    comm = MPI.COMM_WORLD\n",
    "    rank = MPI.Comm_rank(comm)\n",
    "    color = mod(rank,2)\n",
    "    key = 0\n",
    "    newcomm = MPI.Comm_split(comm, color, key)\n",
    "    # Only odd rank ids will execute this\n",
    "    if color == 1\n",
    "        newrank = MPI.Comm_rank(newcomm)\n",
    "        newroot = 0\n",
    "        sndbuf = [rand(1:4)]\n",
    "        if newrank == newroot\n",
    "            rcvbuf = [0]\n",
    "        else\n",
    "            rcvbuf = nothing\n",
    "        end\n",
    "        println(\"[rank $rank, newrank $newrank] I am contributing with $sndbuf\")\n",
    "        MPI.Reduce!(sndbuf, rcvbuf, +, newcomm;root=newroot)\n",
    "        println(\"[rank $rank, newrank $newrank] The sum is $rcvbuf\")\n",
    "    end\n",
    "end\n",
    "run(`$(mpiexec()) -np 4 julia --project=. -e $code`);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d465ebce",
   "metadata": {},
   "source": [
    "Try to run the code without splitting the communicator. I.e., replace `newcomm = MPI.Comm_split(comm, color, key)` with `newcomm = comm`. Try to figure out what will happen before executing the code."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d334aea1",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "- MPI also defines operations involving several processes called, collective operations.\n",
    "- These are provided both for convenience and performance.\n",
    "- The semantics are equivalent to \"standard mode\" `MPI_Send`, but there are also non-blocking versions (not discussed in this notebook).\n",
    "- Discovering message sizes is often done by communicating the message size, instead of using `MPI_Probe`.\n",
    "- Finally, we discussed MPI communicators. They provide two key features: isolated communication context and creating groups of processes. They are useful, for instance, to combine different libraries using MPI in the same application, and to use collective operations in a subset of the processes.\n",
    "\n",
    "After learning this material and the previous MPI notebook, you have a solid basis to start implementing sophisticated parallel algorithms using MPI."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "843b40cd",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c2045d9",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Implement the parallel matrix-matrix multiplication (Algorithm 3) using MPI collectives instead of point-to-point communication. I.e., this is the same exercise as in previous notebook, but using different functions for communication."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e8f6e6a",
   "metadata": {},
   "source": [
    "# License\n",
    "\n",
    "This notebook is part of the course [Programming Large Scale Parallel Systems](https://www.francescverdugo.com/XM_40017) at Vrije Universiteit Amsterdam and may be used under a [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) license."
   ]
  }
 ],
 "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
}
