gpu_intro.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "gpu_intro.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/evaneschneider/parallel-programming/blob/master/gpu_intro.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "metadata": {
        "id": "QicF_avFwaKK",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Introduction to GPU programming with Numba"
      ]
    },
    {
      "metadata": {
        "id": "ZfMrCWC-N_7S",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "This notebook borrows *heavily* from [seibert's 2018 gtc numba tutorial](https://github.com/ContinuumIO/gtc2018-numba). I highly recommend that tutorial in its entireity if you want more practice with Numba and GPUs."
      ]
    },
    {
      "metadata": {
        "id": "ePRO1HsDwkp2",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "Yesterday we discussed the principles of parallel programming, and explored the key aspects of using Numba - the `@jit` decorator, benchmarking, and the `@vectorize` decorator for Numpy UFuncs. Today we are going to expand on that basis and use Numba to do parallel calculations in python by taking advantage of Numba's GPU interface (and Google's free GPUs - thanks Colaboratory!)."
      ]
    },
    {
      "metadata": {
        "id": "x5CRqd4awPqx",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import math\n",
        "from numba import vectorize, cuda\n",
        "from matplotlib import pyplot as plt\n",
        "%matplotlib inline"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "drQnGa0vfCha",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### Problem 0 - Accessing the GPU\n",
        "\n",
        "**0a)** In order to run Numba functions using the GPU, we have to do a couple of things. First, go to the Runtime menu, click on 'Change Runtime Type', and in the pop-up box, under 'Hardware Accelerator', select 'GPU'. Save the Runtime.\n",
        "\n",
        "**0b)** Ideally, that's all we should have to do. But in practice, even though the CUDA libararies are installed, for some reason Colab usually can't find them. So, we'll figure out where they are, and then point Colab to them."
      ]
    },
    {
      "metadata": {
        "id": "qYzz2fEegbpG",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "!find / -iname 'libdevice'\n",
        "!find / -iname 'libnvvm.so'"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "PqW3vHXPgjqk",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "Paste the location of the libraries into the following code box (if it's different, otherwise you can just run the code):"
      ]
    },
    {
      "metadata": {
        "id": "na0O3h9hgoVa",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import os\n",
        "os.environ['NUMBAPRO_LIBDEVICE'] = \"/usr/local/cuda-10.0/nvvm/libdevice\"\n",
        "os.environ['NUMBAPRO_NVVM'] = \"/usr/local/cuda-10.0/nvvm/lib64/libnvvm.so\""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "g007BPp8fHtO",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "And that should do it! Okay, now that we've pointed Numba to the correct libraries, let's get going. To start, we are going to return to the first function we created yesterday - the vector add."
      ]
    },
    {
      "metadata": {
        "id": "JV9ELUHwxnZh",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### Problem 1 - Vector Addition on GPUs"
      ]
    },
    {
      "metadata": {
        "id": "1veva5o3sYB4",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "The simplest way to access the GPU through Numba is to return to our vectorized ufunc from yesterday. As you may recall, Numpy Universal Functions operate on vectors, or arrays. If we specify the `cuda` target, Numba will automatically write a CUDA kernel for us, and run the function on the GPU! Let's try it out:"
      ]
    },
    {
      "metadata": {
        "id": "Klw3GHXtwWyW",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "@vectorize(['int64(int64, int64)'], target='cuda')\n",
        "def add_ufunc(x, y):\n",
        "    return x + y"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "b9UKqsW7gOda",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "x = np.arange(10)\n",
        "y = 2 * x"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "NGQhcF7xtT-B",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "add_ufunc(x, y)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "noAuRtx_tpil",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "Cool, it worked! But what actually just happened? Well, a lot of things. Numba automatically:\n",
        "+ Compiled a CUDA kernel to execute the ufunc operation in parallel over all the input elements.\n",
        "+ Allocated GPU memory for the inputs and the output.\n",
        "+ Copied the input data to the GPU.\n",
        "+ Executed the CUDA kernel with the correct kernel dimensions given the input sizes.\n",
        "+ Copied the result back from the GPU to the CPU.\n",
        "+ Returned the result as a NumPy array on the host."
      ]
    },
    {
      "metadata": {
        "id": "Kw7gZF1muGpy",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**1a)** Determine how fast the CUDA addition function is. Compare that to a function compiled for the CPU. How does the GPU do?\n",
        "\n",
        "*You'll probably want to write two functions with separate names to compare them.*"
      ]
    },
    {
      "metadata": {
        "id": "ZKCwRcCEtXcd",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "C0YFAGN4ROCr",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here\n",
        "def add_ufunc_cpu(x, y):\n",
        "    return x + y\n",
        "\n",
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "AkMt5Fy0ugRm",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**1b)** Wow, the GPU is a LOT slower! Why might that be?\n",
        "\n",
        "*Try to think of several reasons.*"
      ]
    },
    {
      "metadata": {
        "id": "g8BsBv1_uzqM",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "Add text here.\n"
      ]
    },
    {
      "metadata": {
        "id": "FJVE9XZSQSMj",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "###Problem 2 - Memory Management\n",
        "\n",
        "As we saw in the last problem, Numba can automatically handle transferring data to and from the GPU for us. However, that's not always what we want. Sometimes we will want to perform several functions in a row on the GPU without transferring the data back to the CPU in between."
      ]
    },
    {
      "metadata": {
        "id": "ZltkamG0KFYv",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**2a)** Remake the addition ufunc to operate on and return 32 bit floats, and target the GPU."
      ]
    },
    {
      "metadata": {
        "id": "OqWl14DCKKso",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here\n",
        "def add_ufunc(x, y):\n",
        "    return x + y"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "8-ZpffA3KoV1",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "Now, let's give it a bit more work to do:"
      ]
    },
    {
      "metadata": {
        "id": "dS1dCw0PKTdj",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "n = 100000\n",
        "x = np.arange(n).astype(np.float32)\n",
        "y = 2 * x"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "ABw7nZWFKs7b",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "As we saw in the last problem, copying the data to and from the GPU for every function is not necessarily the most efficient way to use the GPU. To address this, Numba provides the `to_device` function in the `cuda` module to allocate and copy arrays to the GPU:"
      ]
    },
    {
      "metadata": {
        "id": "p3eu33UFuVfM",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "x_device = cuda.to_device(x)\n",
        "y_device = cuda.to_device(y)\n",
        "\n",
        "print(x_device)\n",
        "print(x_device.shape)\n",
        "print(x_device.dtype)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "f5U0-goyLG7t",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "`x_device` and `y_device` are now Numba \"device arrays\" that are in many ways equivalent to Numpy ndarrays except that they live in the GPU's global memory, rather than on the CPU. These device arrays can be passed to Numba cuda functions just the way Numpy arrays can, but without the memory copying overhead."
      ]
    },
    {
      "metadata": {
        "id": "QCF1pxsmLwJF",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**2b)** Try out your function using host vs device arrays. How does the time compare?"
      ]
    },
    {
      "metadata": {
        "id": "JNJHfuSSMAf_",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "R9Y59-0zIUno",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "nhqU-TCKMQZB",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "You should see a big performance improvement already, but we are still allocating a device array for the output of the ufunc and copying it back to the host. We can create an output buffer on the GPU with the numba.cuda.device_array() function:"
      ]
    },
    {
      "metadata": {
        "id": "ShhIJX80L_Qh",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "out_device = cuda.device_array(shape=(n,), dtype=np.float32)  # does not initialize the contents, much like np.empty()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "T4B1aBILMc_E",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "And then we can use a special out keyword argument to the ufunc to specify the output buffer:"
      ]
    },
    {
      "metadata": {
        "id": "umP9Ab7YMY8H",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "%timeit add_ufunc(x_device, y_device, out=out_device)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "hhF9NpXONd5E",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "You should see an even bigger improvement. Once we've finished all of our calculations on the GPU, we can copy the array back from the device using the `copy_to_host` method:"
      ]
    },
    {
      "metadata": {
        "id": "8QAMR4qKMgoS",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "out_host = out_device.copy_to_host()\n",
        "print(out_host[:10])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "h7mdK5DUVNgi",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**2c)** Remake a new version of the addition ufunc with 32bit floats that targets the cpu. Compare the resulting time to execute with the gpu version you just timed."
      ]
    },
    {
      "metadata": {
        "id": "gTG333oXUumf",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here\n",
        "def add_ufunc_cpu(x, y):\n",
        "    return x + y\n",
        "%timeit add_ufunc_cpu(x, y)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "Uu0I4ck3VzA-",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**2d)** Now go back and try the two functions (gpu and cpu) with even larger arrays. When does the GPU start to win? Does the execution time on the GPU scale with the number of array elements the same way that the CPU version does?"
      ]
    },
    {
      "metadata": {
        "id": "M5s3oC-nnYtA",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "If your result is like mine, you may have seen a slight timing advantage in the GPU version with a million array elements, but it was close (and that wasn't even counting the data transfer time). That's because we're still not giving the GPU enough work to keep all those cores busy all the time! By the time we hit 10 million, the GPU was clearly winning. The time it took the CPU function continued to increase linearly with array size, but the GPU function time increased much more slowly."
      ]
    },
    {
      "metadata": {
        "id": "eeve4KxkO2OZ",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**2e)** Let's practice some more memory management. Given the following ufuncs:"
      ]
    },
    {
      "metadata": {
        "id": "5VeKgSFXNxxa",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "@vectorize(['float32(float32, float32)'], target='cuda')\n",
        "def add_ufunc(x, y):\n",
        "    return x + y\n",
        "\n",
        "@vectorize(['float32(float32, float32, float32)'], target='cuda')\n",
        "def make_pulses(i, period, amplitude):\n",
        "    return max(math.sin(i / period) - 0.3, 0.0) * amplitude"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "PaAgMJf3PITJ",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "Convert the following code to use device allocations so that there are only host<->device copies at the beginning and end. Then benchmark the performance change.\n",
        "\n",
        "*Hint #1: You may want to keep the following code in place for reference, and create the version with GPU arrays in the additional cells below.*\n",
        "\n",
        "*Hint #2: Determine how many arrays you will need on the device.*"
      ]
    },
    {
      "metadata": {
        "id": "Jlk_PrSXoH_u",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "n = 100000\n",
        "noise = (np.random.normal(size=n) * 3).astype(np.float32)\n",
        "t = np.arange(n, dtype=np.float32)\n",
        "period = n / 23"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "_Bbpugz2ZpqK",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "pulses = make_pulses(t, period, 100.0)\n",
        "waveform = add_ufunc(pulses, noise)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "AbMU8BVfPUeO",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "plt.plot(waveform)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "EpDq3V6eO_CY",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "RNcWlpFWPhoz",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "fjXXoHB1gG9o",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### Problem 3 - Writing Cuda Kernels\n",
        "\n",
        "While targeting ufuncs with the `cuda` syntax is the most straightforward way to access the GPU with Numba, it may not be flexible enough for your needs. If you want to write a more detailed GPU program, at some point you are probably going to need to write CUDA kernels.\n",
        "\n",
        "As discussed in the lecture, the CUDA programming model allows you to abstract the GPU hardware into a software model composed of a **grid** containing **blocks** of **threads**. These threads are the smallest individual unit in the programming model, and they execute together in groups (traditionally called **warps**, consisting of 32 thread each). Determiming the best size for your grid of thread blocks is a complicated problem that often depends on the specific algorithm and hardware you're using, but here a few good rules of thumb:\n",
        "+ the size of a block should be a multiple of 32 threads, with typical block sizes between 128 and 512 threads per block.\n",
        "+ the size of the grid should ensure the full GPU is utilized where possible. Launching a grid where the number of blocks is 2x-4x the number of \"multiprocessors\" on the GPU is a good starting place. Something in the range of 20 - 100 blocks is usually a good starting point.\n",
        "+ The CUDA kernel launch overhead does depend on the number of blocks, so it may not be best to launch a grid where the number of threads equals the number of input elements when the input size is very big. We'll show a pattern for dealing with large inputs below.\n",
        "\n",
        "As a first example, let's return to our vector addition function, but this time, we'll target it with the `cuda.jit` decorator:"
      ]
    },
    {
      "metadata": {
        "id": "aHogP9LBdCSf",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "@cuda.jit\n",
        "def add_kernel(x, y, out):\n",
        "    tidx = cuda.threadIdx.x # this is the unique thread ID within a 1D block\n",
        "    bidx = cuda.blockIdx.x  # Similarly, this is the unique block ID within the 1D grid\n",
        "\n",
        "    block_dimx = cuda.blockDim.x  # number of threads per block\n",
        "    grid_dimx = cuda.gridDim.x    # number of blocks in the grid\n",
        "    \n",
        "    start = tidx + bidx * block_dimx\n",
        "    stride = block_dimx * grid_dimx\n",
        "\n",
        "    # assuming x and y inputs are same length\n",
        "    for i in range(start, x.shape[0], stride):\n",
        "        out[i] = x[i] + y[i]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "xyxOba6DkD0G",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "\n",
        "That's a lot more typing than our ufunc example, and it is much more limited: it only works on 1D arrays, it doesn't verify input sizes match, etc. Most of the function is spent figuring out how to turn the block and grid indices and dimensions into unique offsets in the input arrays. The pattern of computing a starting index and a stride is a common way to ensure that your grid size is independent of the input size. The striding will maximize bandwidth by ensuring that threads with consecuitive indices are accessing consecutive memory locations as much as possible. Thread indices beyond the length of the input (x.shape[0], since x is a NumPy array) automatically skip over the for loop.\n",
        "\n",
        "Let's call the function now on some data:"
      ]
    },
    {
      "metadata": {
        "id": "I4e4bd0fhoR9",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "n = 100000\n",
        "x = np.arange(n).astype(np.float32)\n",
        "y = 2 * x\n",
        "out = np.empty_like(x)\n",
        "\n",
        "threads_per_block = 128\n",
        "blocks_per_grid = 30\n",
        "\n",
        "add_kernel[blocks_per_grid, threads_per_block](x, y, out)\n",
        "print(out[:10])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "phlsKQvgkqR1",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "The calling syntax is designed to mimic the way CUDA kernels are launched in C, where the number of blocks per grid and threads per block are specified in the square brackets, and the arguments to the function are specified afterwards in parentheses.\n",
        "\n",
        "Note that, unlike the ufunc, the arguments are passed to the kernel as full NumPy arrays. A thread within the kernel can access any element in the array it wants, regardless of its position in the thread grid. This is why CUDA kernels are significantly more powerful than ufuncs. (But with great power, comes a greater amount of typing...)\n",
        "\n",
        "Numba has created some [helper functions](http://numba.pydata.org/numba-doc/dev/cuda/kernels.html#absolute-positions) to cut down on the typing. We can write the previous kernel much more simply as:"
      ]
    },
    {
      "metadata": {
        "id": "W8idvpu1kkX6",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "@cuda.jit\n",
        "def add_kernel(x, y, out):\n",
        "    start = cuda.grid(1)      # the 1 argument means a one dimensional thread grid, this returns a single value\n",
        "    stride = cuda.gridsize(1) # ditto\n",
        "\n",
        "    # assuming x and y inputs are same length\n",
        "    for i in range(start, x.shape[0], stride):\n",
        "        out[i] = x[i] + y[i]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "dWZ81VlKmeZT",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "As before, using NumPy arrays forces Numba to allocate GPU memory, copy the arguments to the GPU, run the kernel, then copy the argument arrays back to the host. This not very efficient, so you will often want to allocate device arrays.\n",
        "\n",
        "**3a)** Allocate device arrays for x, y, and the output, then try out your new Cuda kernel using the pre-copied device arrays. Compare the time to a version without moving the data first."
      ]
    },
    {
      "metadata": {
        "id": "SYFkeEw4mqvM",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "ykopd4uNmvZK",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "OL5K9JFom0Xn",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "0xCY7lSNuOhj",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "#### Atomic Operations and avoiding Race Conditions\n",
        "\n",
        "CUDA, like many general purpose parallel execution frameworks, makes it possible to have race conditions in your code. A race condition in CUDA arises when threads read or write a memory location that might be modified by another independent thread. Generally speaking, you need to worry about:\n",
        "+ read-after-write hazards: One thread is reading a memory location at the same time another thread might be writing to it.\n",
        "+ write-after-write hazards: Two threads are writing to the same memory location, and only one write will be visible when the kernel is complete.\n",
        "\n",
        "A common strategy to avoid both of these hazards is to organize your CUDA kernel algorithm such that each thread has exclusive responsibility for unique subsets of output array elements, and/or to never use the same array for both input and output in a single kernel call. (Iterative algorithms can use a double-buffering strategy if needed, and switch input and output arrays on each iteration.)\n",
        "However, there are many cases where different threads need to combine results. Consider something very simple, like: \"every thread increments a global counter.\" Implementing this in your kernel requires each thread to:\n",
        "1. Read the current value of a global counter.\n",
        "2. Compute counter + 1.\n",
        "3. Write that value back to global memory.\n",
        "However, there is no guarantee that another thread has not changed the global counter between steps 1 and 3. To resolve this problem, CUDA provides \"atomic operations\" which will read, modify and update a memory location in one, indivisible step. Numba supports several of these functions, [described here](http://numba.pydata.org/numba-doc/dev/cuda/intrinsics.html#supported-atomic-operations).\n",
        "\n",
        "As an example, let's make a thread counter kernel:"
      ]
    },
    {
      "metadata": {
        "id": "6rHpfkqhvC7i",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "@cuda.jit\n",
        "def thread_counter_race_condition(global_counter):\n",
        "    global_counter[0] += 1  # This is bad\n",
        "    \n",
        "@cuda.jit\n",
        "def thread_counter_safe(global_counter):\n",
        "    cuda.atomic.add(global_counter, 0, 1)  # Safely add 1 to offset 0 in global_counter array"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "1Q6NbcwFvGzf",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# This gets the wrong answer\n",
        "global_counter = cuda.to_device(np.array([0], dtype=np.int32))\n",
        "thread_counter_race_condition[64, 64](global_counter)\n",
        "\n",
        "print('Should be %d:' % (64*64), global_counter.copy_to_host())"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "SB7ZeUGCvrGw",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# This works correctly\n",
        "global_counter = cuda.to_device(np.array([0], dtype=np.int32))\n",
        "thread_counter_safe[64, 64](global_counter)\n",
        "\n",
        "print('Should be %d:' % (64*64), global_counter.copy_to_host())"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "39f5mwaMv0C9",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**3b)** Let's practice writing a function that requires an atomic operation - a histogramming kernel. This will take an array of input data, a range and a number of bins, and count how many of the input data elements land in each bin. Below is an example CPU implementation of histogramming:"
      ]
    },
    {
      "metadata": {
        "id": "mFBD1af_m-n5",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def cpu_histogram(x, xmin, xmax, histogram_out):\n",
        "    '''Increment bin counts in histogram_out, given histogram range [xmin, xmax).'''\n",
        "    # Note that we don't have to pass in nbins explicitly, because the size of histogram_out determines it\n",
        "    nbins = histogram_out.shape[0]\n",
        "    bin_width = (xmax - xmin) / nbins\n",
        "    \n",
        "    # This is a very slow way to do this with NumPy, but looks similar to what you will do on the GPU\n",
        "    for element in x:\n",
        "        bin_number = np.int32((element - xmin)/bin_width)\n",
        "        if bin_number >= 0 and bin_number < histogram_out.shape[0]:\n",
        "            # only increment if in range\n",
        "            histogram_out[bin_number] += 1"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "eMScRD4Orc-3",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "x = np.random.normal(size=10000, loc=0, scale=1).astype(np.float32)\n",
        "xmin = np.float32(-4.0)\n",
        "xmax = np.float32(4.0)\n",
        "histogram_out = np.zeros(shape=10, dtype=np.int32)\n",
        "\n",
        "cpu_histogram(x, xmin, xmax, histogram_out)\n",
        "\n",
        "histogram_out"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "Lr158nPNwLlH",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "In the space below, create a cuda version of this kernel, then run it to check that you get the same answer as the CPU version.\n",
        "\n",
        "*Hint: You can use much of the same syntax that we used in the CUDA addition kernel.*"
      ]
    },
    {
      "metadata": {
        "id": "Ij3ZUgvBrdcq",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "@cuda.jit\n",
        "def cuda_histogram(x, xmin, xmax, histogram_out):\n",
        "    '''Increment bin counts in histogram_out, given histogram range [xmin, xmax).'''\n",
        "    \n",
        "    # add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "Tm17T1Z3s06f",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "7_5MBLU6yUZC",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### Problem 4 - Return to the Fractals!\n",
        "\n",
        "Yesterday we defined two functions to create an instance of the Julia set:"
      ]
    },
    {
      "metadata": {
        "id": "EMa5G0pus5ZI",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def julia(x, y, max_iters):\n",
        "    \"\"\"\n",
        "    Given the real and imaginary parts of a complex number,\n",
        "    determine if it is a candidate for membership in the Julia\n",
        "    set given a fixed number of iterations.\n",
        "    \"\"\"\n",
        "    i = 0\n",
        "    c = complex(-0.8, 0.156)\n",
        "    a = complex(x,y)\n",
        "    for i in range(max_iters):\n",
        "        a = a*a + c\n",
        "        if (a.real*a.real + a.imag*a.imag) > 1000:\n",
        "            return 0\n",
        "    return 255"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "PjV_qG1oyiPH",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def create_fractal(min_x, max_x, min_y, max_y, image, iters):\n",
        "    height = image.shape[0]\n",
        "    width = image.shape[1]\n",
        "\n",
        "    pixel_size_x = (max_x - min_x) / width\n",
        "    pixel_size_y = (max_y - min_y) / height\n",
        "    for x in range(width):\n",
        "        real = min_x + x * pixel_size_x\n",
        "        for y in range(height):\n",
        "            imag = min_y + y * pixel_size_y\n",
        "            color = julia(real, imag, iters)\n",
        "            image[y, x] = color\n",
        "\n",
        "    return image"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "VUokQRuayu1m",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "image = np.zeros((500, 750), dtype=np.uint8)\n",
        "create_fractal(-2.0, 2.0, -1.0, 1.0, image, 200)\n",
        "plt.imshow(image)\n",
        "plt.viridis()\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "SnyP0npyzWhF",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "In order to turn this into a GPU implementation, we'd like to have a kernel function (create_fractal) call another function (julia) on the device. Numba has a way of specifying functions that will be called from within a kernel by passing the `cuda.jit` decorator an argument:"
      ]
    },
    {
      "metadata": {
        "id": "5GqaXTVay11s",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "@cuda.jit(device=True)\n",
        "def julia(x, y, max_iters):\n",
        "    \"\"\"\n",
        "    Given the real and imaginary parts of a complex number,\n",
        "    determine if it is a candidate for membership in the Julia\n",
        "    set given a fixed number of iterations.\n",
        "    \"\"\"\n",
        "    i = 0\n",
        "    c = complex(-0.8, 0.156)\n",
        "    a = complex(x,y)\n",
        "    for i in range(max_iters):\n",
        "        a = a*a + c\n",
        "        if (a.real*a.real + a.imag*a.imag) > 1000:\n",
        "            return 0\n",
        "    return 255"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "bKbyLkba5Haw",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**Multi-dimensional grids**\n",
        "\n",
        "For some problems, it makes sense to define a two- or three-dimensional grid of thread blocks. That way, when you're indexing a single thread, you can map it to, say, the pixel position in an image. Multi-dimensional grids are created by passing tuples to the kernel function. You can ensure that you launch a big enough grid by calculating the size of each dimension as a function of the array size and number of threads per block:"
      ]
    },
    {
      "metadata": {
        "id": "DCkBC9P-56Pc",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "threadsperblock = 16\n",
        "xblocks = (image.shape[1] + (threadsperblock - 1)) // threadsperblock\n",
        "yblocks = (image.shape[0] + (threadsperblock - 1)) // threadsperblock"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "Pr73XJ-D6BAP",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "Then, within a kernel, you can determine the absolute thread position by calling the `grid` helper function, as in `x, y = cuda.grid(2)`."
      ]
    },
    {
      "metadata": {
        "id": "AOR3AsLt6Sb8",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**4a)** Modify the create_fractal function to launch as a kernel on the GPU and call your new device function, julia. Use a 2D grid of thread blocks to launch the kernel and determine which threads are responsible for each pixel in the image."
      ]
    },
    {
      "metadata": {
        "id": "6V-cRn4Iz2CA",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "VN5xFf7U1t2K",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "w7QGkLSY1iZf",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# add code here"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "7RbFpxB17Czg",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "**4b)** How does the GPU implementation compare to the CPU version? How about to yesterday's Numba CPU version?"
      ]
    },
    {
      "metadata": {
        "id": "USY6UCyZ7q8f",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "Write your answer here."
      ]
    },
    {
      "metadata": {
        "id": "6H-ozx9JDhMu",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}