Experimental/intro_to_gpu_python.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:2569c5a47def2a2a0d5ab06c0d4fdb8e9c5dfdbb688e232f469798a428f64ee5"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# Introduction to Python GPU Programming with Numba and NumbaPro"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Misc. import\n",
      "from IPython.display import Image "
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "skip"
      }
     },
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# Numba\n",
      "* Opensource BSD license\n",
      "* Basic CUDA GPU JIT compilation\n",
      "* OpenCL support coming"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numba\n",
      "print \"numba\", numba.__version__"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "numba 0.17.0\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# NumbaPro\n",
      "* Commercial\n",
      "* Contains library extensions and extra multicore and GPU features:\n",
      "    * Parallel vectorize\n",
      "    * GPU vectorize, guvectorize\n",
      "    * CUDA library bindings: cuRAND, cuBLAS, cuFFT, cuSparse\n",
      "* NumbaPro namespace re-exports all public symbols in Numba\n",
      "* We will import from numbapro whenever the feature is NumbaPro only"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numbapro\n",
      "print \"numbapro\", numbapro.__version__"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "numbapro 0.17.1\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# The CUDA GPU\n",
      "\n",
      "- A massively parallel processor (many cores)\n",
      "    - 100 threads, 1,000 threads, and more\n",
      "- optimized for data throughput\n",
      "    - Simple (shallow) cache hierarchy\n",
      "    - Best with manual caching!\n",
      "    - Cache memory is called shared memory and it is addressable\n",
      "- CPU is latency optimized\n",
      "    - Deep cache hierarchy\n",
      "    - L1, L2, L3 caches\n",
      "- GPU execution model is different\n",
      "- GPU forces you to think and program *in parallel*"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get all the imports we need\n",
      "import numba.cuda\n",
      "import numpy as np\n",
      "import math\n",
      "\n",
      "my_gpu = numba.cuda.get_current_device()\n",
      "print \"Running on GPU:\", my_gpu.name\n",
      "cores_per_capability = {\n",
      "    1: 8,\n",
      "    2: 32,\n",
      "    3: 192,\n",
      "}\n",
      "cc = my_gpu.compute_capability\n",
      "print \"Compute capability: \", \"%d.%d\" % cc, \"(Numba requires >= 2.0)\"\n",
      "majorcc = cc[0]\n",
      "print \"Number of streaming multiprocessor:\", my_gpu.MULTIPROCESSOR_COUNT\n",
      "cores_per_multiprocessor = cores_per_capability[majorcc]\n",
      "print \"Number of cores per mutliprocessor:\", cores_per_multiprocessor\n",
      "total_cores = cores_per_multiprocessor * my_gpu.MULTIPROCESSOR_COUNT\n",
      "print \"Number of cores on GPU:\", total_cores"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Running on GPU: GeForce GTX 760\n",
        "Compute capability:  3.0 (Numba requires >= 2.0)\n",
        "Number of streaming multiprocessor: 6\n",
        "Number of cores per mutliprocessor: 192\n",
        "Number of cores on GPU: 1152\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# High-level Array-Oriented Style\n",
      "\n",
      "- Use NumPy array as a unit of computation\n",
      "- Use NumPy universal function (ufunc) as an abstraction of computation and scheduling\n",
      "- ufuncs are elementwise functions\n",
      "- If you use NumPy, you are using ufuncs"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.sin, \"is of type\", type(np.sin)\n",
      "print np.add, \"is of type\", type(np.add)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<ufunc 'sin'> is of type <type 'numpy.ufunc'>\n",
        "<ufunc 'add'> is of type <type 'numpy.ufunc'>\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "# Vectorize\n",
      "\n",
      "- generate a ufunc from a python function\n",
      "- converts scalar function to elementwise array function\n",
      "- Numba provides CPU support\n",
      "- NumbaPro provides GPU support"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# CPU version\n",
      "@numba.vectorize(['float32(float32, float32)',\n",
      "                  'float64(float64, float64)'], target='cpu')\n",
      "def cpu_sincos(x, y):\n",
      "    return math.sin(x) * math.cos(y)\n",
      "\n",
      "# CUDA version\n",
      "@numbapro.vectorize(['float32(float32, float32)',\n",
      "                     'float64(float64, float64)'], target='gpu')\n",
      "def gpu_sincos(x, y):\n",
      "    return math.sin(x) * math.cos(y)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "```\n",
      "@numba.vectorize(<list of signatures>, target=<'cpu', 'gpu'>)\n",
      "```\n",
      "\n",
      "- A ufunc can be overloaded to work on multiple type signatures\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "### Test it out\n",
      "\n",
      "- 2 input arrays\n",
      "- 1 output array\n",
      "- 1 million doubles (8 MB) per array\n",
      "- Total 24 MB of data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Generate data\n",
      "n = 1000000\n",
      "x = np.linspace(0, np.pi, n)\n",
      "y = np.linspace(0, np.pi, n)\n",
      "\n",
      "# Check result\n",
      "np_ans = np.sin(x) * np.cos(y)\n",
      "nb_cpu_ans = cpu_sincos(x, y)\n",
      "nb_gpu_ans = gpu_sincos(x, y)\n",
      "\n",
      "print \"CPU vectorize correct: \", np.allclose(nb_cpu_ans, np_ans)\n",
      "print \"GPU vectorize correct: \", np.allclose(nb_gpu_ans, np_ans)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU vectorize correct:  "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "True\n",
        "GPU vectorize correct:  "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "True\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "## Benchmark"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"NumPy\"\n",
      "%timeit np.sin(x) * np.cos(y)\n",
      "\n",
      "print \"CPU vectorize\"\n",
      "%timeit cpu_sincos(x, y)\n",
      "\n",
      "print \"GPU vectorize\"\n",
      "%timeit gpu_sincos(x, y)\n",
      "\n",
      "# Optional cleanup \n",
      "del x, y"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "NumPy\n",
        "10 loops, best of 3: 52.9 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "CPU vectorize\n",
        "10 loops, best of 3: 43.9 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "GPU vectorize\n",
        "100 loops, best of 3: 15.7 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "- CPU vectorize time is similar to pure NumPy time because ``sin()`` and ``cos()`` calls dominate the time.\n",
      "- GPU vectorize is a lot faster\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "\n",
      "## Behind the scence\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "\n",
      "### Automatic memory transfer\n",
      "\n",
      "- NumPy arrays are automatically transferred\n",
      "    - CPU -> GPU\n",
      "    - GPU -> CPU"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "\n",
      "### Automatic work scheduling\n",
      "\n",
      "- The work is distributed the across all threads on the GPU\n",
      "- The GPU hardware handles the scheduling\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "\n",
      "### Automatic GPU memory management\n",
      "\n",
      "- GPU memory is tied to object lifetime\n",
      "- freed automatically"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "# Another Vectorize Example"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@numba.vectorize(['float32(float32, float32, float32, float32)'])\n",
      "def cpu_powers(p, q, r, s):\n",
      "    return math.sqrt(p ** 2 + q ** 3 + r ** 4 + s ** 5)\n",
      "\n",
      "@numbapro.vectorize(['float32(float32, float32, float32, float32)'], target='gpu')\n",
      "def gpu_powers(p, q, r, s):\n",
      "    return math.sqrt(p ** 2 + q ** 3 + r ** 4 + s ** 5)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Generate data\n",
      "n = 5000000\n",
      "p = np.random.random(n).astype(np.float32)\n",
      "q = np.random.random(n).astype(np.float32)\n",
      "r = np.random.random(n).astype(np.float32)\n",
      "s = np.random.random(n).astype(np.float32)\n",
      "\n",
      "# Check results\n",
      "np_ans = np.sqrt(p ** 2 + q ** 3 + r ** 4 + s ** 5)\n",
      "cpu_ans = cpu_powers(p, q, r, s)\n",
      "gpu_ans = gpu_powers(p, q, r, s)\n",
      "print \"CPU vectorize correct\", np.allclose(np_ans, cpu_ans)\n",
      "print \"GPU vectorize correct\", np.allclose(np_ans, gpu_ans)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU vectorize correct "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "True\n",
        "GPU vectorize correct "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "True\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "## Benchmark"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"NumPy\"\n",
      "%timeit np.sqrt(p ** 2 + q ** 3 + r ** 4 + s ** 5)\n",
      "print \"CPU vectorize\"\n",
      "%timeit cpu_powers(p, q, r, s)\n",
      "print \"GPU vectorize\"\n",
      "%timeit gpu_powers(p, q, r, s)\n",
      "\n",
      "# Optional cleanup \n",
      "del p, q, r, s"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "NumPy\n",
        "1 loops, best of 3: 548 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "CPU vectorize\n",
        "10 loops, best of 3: 55.6 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "GPU vectorize\n",
        "10 loops, best of 3: 51.1 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "* NumPy is slower than CPU vectorize (likely due to allocation of temporaries)\n",
      "* CPU version is slower than GPU version because of multiple cores running (even though copying to card and back is costly)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# Generalized Universal Function (guvectorize)\n",
      "\n",
      "- Vectorize is limited to scalar arguments in the core function\n",
      "- GUVectorize accepts array arguments"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@numbapro.guvectorize(['void(float32[:,:], float32[:,:], float32[:,:])'],\n",
      "                      '(m, n),(n, p)->(m, p)', target='gpu')\n",
      "def batch_matrix_mult(a, b, c):\n",
      "    for i in range(c.shape[0]):\n",
      "        for j in range(c.shape[1]):\n",
      "            tmp = 0\n",
      "            for n in range(a.shape[1]):\n",
      "                 tmp += a[i, n] * b[n, j]\n",
      "            c[i, j] = tmp"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "```\n",
      "@numbapro.guvectorize(<list of function signatures>, <a string to desc the shape signature>, target=<'cpu', 'gpu'>)\n",
      "```"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = np.arange(1.0, 10.0, dtype=np.float32).reshape(3,3)\n",
      "b = np.arange(1.0, 10.0, dtype=np.float32).reshape(3,3)\n",
      "\n",
      "# Use the builtin matrix_multiply in NumPy for CPU test\n",
      "import numpy.core.umath_tests as ut\n",
      "\n",
      "# Check result\n",
      "print 'NumPy result'\n",
      "np_ans = ut.matrix_multiply(a, b)\n",
      "print np_ans\n",
      "\n",
      "print 'Numba GPU result'\n",
      "gpu_ans = batch_matrix_mult(a, b)\n",
      "print gpu_ans\n",
      "\n",
      "assert np.allclose(np_ans, gpu_ans)\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "NumPy result\n",
        "[[  30.   36.   42.]\n",
        " [  66.   81.   96.]\n",
        " [ 102.  126.  150.]]\n",
        "Numba GPU result\n",
        "[[  30.   36.   42.]\n",
        " [  66.   81.   96.]\n",
        " [ 102.  126.  150.]]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "### Test it out\n",
      "\n",
      "- batch multiply two 4 million 2x2 matrices "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 4000000\n",
      "dim = 2\n",
      "a = np.random.random(n * dim * dim).astype(np.float32).reshape(n, dim, dim)\n",
      "b = np.random.random(n * dim * dim).astype(np.float32).reshape(n, dim, dim)\n",
      "\n",
      "print 'NumPy time'\n",
      "%timeit ut.matrix_multiply(a, b)\n",
      "\n",
      "print 'Numba GPU time'\n",
      "%timeit batch_matrix_mult(a, b)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "NumPy time\n",
        "1 loops, best of 3: 205 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Numba GPU time\n",
        "1 loops, best of 3: 180 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "- GPU time seems to be similar to CPU time\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "\n",
      "## Manually Transfer the data to the GPU\n",
      "\n",
      "- This will let us see the actual compute time without the CPU<->GPU transfer"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dc = numba.cuda.device_array_like(a)\n",
      "da = numba.cuda.to_device(a)\n",
      "db = numba.cuda.to_device(b)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "* ```numba.cuda.device_array_like``` allocate without initialization with the type and shape of another array.\n",
      "    * similar to ```numpy.empty_like(a)```\n",
      "* ```numba.cuda.to_device``` create a GPU copy of the CPU array\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "\n",
      "## Pure compute time on the GPU"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def check_pure_compute_time(da, db, dc):\n",
      "    batch_matrix_mult(da, db, out=dc)\n",
      "    numba.cuda.synchronize()   # ensure the call has completed\n",
      "    \n",
      "%timeit check_pure_compute_time(da, db, dc)\n",
      "del da, db, dc"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 11 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "* Actual compute time is **a lot faster**\n",
      "* PCI-express transfer overhead"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "#### Tips\n",
      "If you have a sequence of ufuncs to apply, pin the data on the GPU by manual transfer"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# NumbaPro CUDA Libraries Bindings\n",
      "\n",
      "- Access to CUDA libraries \n",
      "- Work seamless with NumPy\n",
      "    - auto memory transfer\n",
      "    - managed memory\n",
      "\n",
      "- cuBLAS: CUDA version of BLAS\n",
      "- cuSparse: CUDA sparse matrix support\n",
      "- cuFFT: FFT on CUDA\n",
      "- cuRAND: random number generation on CUDA"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "# An Example with CUDA Libs\n",
      "## Convolution on the GPU"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Import cuFFT"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numbapro.cudalib import cufft"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "#### Misc. imports"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.signal import fftconvolve\n",
      "from scipy import misc, ndimage\n",
      "from matplotlib import pyplot as plt\n",
      "from timeit import default_timer as timer"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Build elementwise complex array multiplication CUDA function"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@numbapro.vectorize(['complex64(complex64, complex64)'], target='gpu')\n",
      "def vmult(a, b):\n",
      "    return a * b"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Prepare image and filter"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "image = misc.face(gray=True).astype(np.float32)\n",
      "\n",
      "laplacian_pts = '''\n",
      "-4 -1 0 -1 -4\n",
      "-1 2 3 2 -1\n",
      "0 3 4 3 0\n",
      "-1 2 3 2 -1\n",
      "-4 -1 0 -1 -4\n",
      "'''.split()\n",
      "\n",
      "laplacian = np.array(laplacian_pts, dtype=np.float32).reshape(5, 5)\n",
      "\n",
      "response = np.zeros_like(image)\n",
      "response[:5, :5] = laplacian\n",
      "\n",
      "print \"Image size: %s\" % (image.shape,)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Image size: (768, 1024)\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Convolution on the CPU"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ts = timer()  # Start Timer\n",
      "cvimage_cpu = fftconvolve(image, laplacian, mode='same')\n",
      "te = timer()  # Stop Timer\n",
      "print 'CPU: %.2fs' % (te - ts)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU: 0.36s\n"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Convolution on the GPU"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "image_complex = image.astype(np.complex64)\n",
      "response_complex = response.astype(np.complex64)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Trigger initialization the cuFFT system.\n",
      "# This takes significant time for small dataset.\n",
      "# We should not be including the time wasted here\n",
      "cufft.FFTPlan(shape=image.shape, itype=np.complex64, otype=np.complex64)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "skip"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 23,
       "text": [
        "<numbapro.cudalib.cufft.api.FFTPlan at 0x11604ead0>"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ts = timer()     # Start timer\n",
      "\n",
      "d_image_complex = numba.cuda.to_device(image_complex)\n",
      "d_response_complex = numba.cuda.to_device(response_complex)\n",
      "\n",
      "# Forward FFT\n",
      "cufft.fft_inplace(d_image_complex)\n",
      "cufft.fft_inplace(d_response_complex)\n",
      "\n",
      "# Multiply the image with teh filter\n",
      "vmult(d_image_complex, d_response_complex, out=d_image_complex)\n",
      "\n",
      "# Inverse FFT\n",
      "cufft.ifft_inplace(d_image_complex)\n",
      "\n",
      "cvimage_gpu = d_image_complex.copy_to_host().real / np.prod(image.shape)\n",
      "\n",
      "te = timer()   # Stop timer\n",
      "\n",
      "print 'GPU: %.2fs' % (te - ts)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "GPU: 0.03s\n"
       ]
      }
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.title('Original')\n",
      "plt.imshow(image, cmap=plt.cm.gray)\n",
      "plt.axis('off')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.title('CPU')\n",
      "plt.imshow(cvimage_cpu, cmap=plt.cm.gray)\n",
      "plt.axis('off')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.title('GPU')\n",
      "plt.imshow(cvimage_gpu, cmap=plt.cm.gray)\n",
      "plt.axis('off')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "-----\n",
      "\n",
      "# Low-Level Approach: @numba.cuda.jit\n",
      "\n",
      "- Numba can generate CUDA functions with the `@jit` decorator\n",
      "- Decorated function follows CUDA execution model "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "## CUDA Execution Model\n",
      "\n",
      "- Kernel functions\n",
      "    - visible to the host CPU\n",
      "    - cannot return any value\n",
      "        - use output argument\n",
      "    - associates to a _grid_\n",
      "- Grid\n",
      "    - a group of blocks\n",
      "    - 1D, 2D, 3D\n",
      "- Blocks\n",
      "    - a group of threads\n",
      "    - 1D, 2D, 3D  \n",
      "- Every thread executes the same kernel\n",
      "    - thread can use the grid, block, thread coordinate system to determine its ID"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Image(url='http://docs.nvidia.com/cuda/cuda-c-programming-guide/graphics/grid-of-thread-blocks.png')"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "## Compiling a CUDA Kernel"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numba import cuda\n",
      "\n",
      "@numba.cuda.jit(\"void(float32[:], float32[:], float32[:])\")\n",
      "def vadd(arr_a, arr_b, arr_out):\n",
      "    tx = cuda.threadIdx.x\n",
      "    bx = cuda.blockIdx.x\n",
      "    bw = cuda.blockDim.x    # number of threads per block\n",
      "    i = tx + bx * bw\n",
      "    if i >= arr_out.size:\n",
      "        return\n",
      "    arr_out[i] = arr_a[i] + arr_b[i]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "### Code Explained\n",
      "\n",
      "#### Define a CUDA kernel with three 1D float32 arrays as args\n",
      "\n",
      "```\n",
      "@numba.cuda.jit(\"void(float32[:], float32[:], float32[:])\")\n",
      "def vadd(arr_a, arr_b, arr_out):\n",
      "```\n",
      "\n",
      "#### Map thread, block coordinate to global position\n",
      "```\n",
      "    tx = cuda.threadIdx.x   # thread label (along x dimension)\n",
      "    bx = cuda.blockIdx.x    # block label (along x dimension)\n",
      "    bw = cuda.blockDim.x    # number of threads in each block (along x dimension)\n",
      "    i = tx + bx * bw        # flattened linear address for each thread\n",
      "```\n",
      "or simplified to:\n",
      "```\n",
      "    i = cuda.grid(1)\n",
      "``` \n",
      "\n",
      "#### Ensure global position is within array size\n",
      "\n",
      "```\n",
      "    if i >= arr_out.size:\n",
      "        return\n",
      "```\n",
      "\n",
      "#### The actual work\n",
      "\n",
      "```\n",
      "    arr_out[i] = arr_a[i] + arr_b[i]\n",
      "```"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "## Launch kernel"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "#### Prepare data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 100\n",
      "a = np.arange(n, dtype=np.float32)\n",
      "b = np.arange(n, dtype=np.float32)\n",
      "c = np.empty_like(a)                 # Must prepare the output array to hold the result\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Calculate thread, block count\n",
      "\n",
      "- thread count is set to **warp size** of the GPU\n",
      "    - Warp size is similar to SIMD vector width on the CPU\n",
      "    - **performance tips**: set thread count to multiple of warp size\n",
      "- block count is ceil(n/thread_ct)\n",
      "\n",
      "**Note:**\n",
      "This will launch more threads than there are elements in the array"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "thread_ct = my_gpu.WARP_SIZE\n",
      "block_ct = int(math.ceil(float(n) / thread_ct))\n",
      "\n",
      "print \"Threads per block:\", thread_ct\n",
      "print \"Block per grid:\", block_ct"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Launch kernel\n",
      "\n",
      "Kernel function object uses the ``__getitem__`` (indexing notation) to configure the grid and block dimensions.\n",
      "\n",
      "```\n",
      "    kernel_function[griddim, blockdim](*args)\n",
      "```\n",
      "\n",
      "- **griddim**\n",
      "    - Number of blocks per grid (grid dimension)\n",
      "    - type: int for 1d or 1,2,3-tuple of ints for 1d, 2d, or 3d respectively\n",
      "\n",
      "- **blockdim**: \n",
      "    - Number of threads per block (blockdim dimension)\n",
      "    - type: int for 1d or 1,2,3-tuple of ints for 1d, 2d, or 3d respectively\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vadd[block_ct, thread_ct](a, b, c)    # Last argument is the output array in this case\n",
      "print c"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "# Example: Matrix Matrix Multiplication\n",
      "\n",
      "- Show manual caching with shared memory\n",
      "- Not the best matrix matrix multiplication implementation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Prepare constants"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numba import float32\n",
      "\n",
      "bpg = 150\n",
      "tpb = 32\n",
      "n = bpg * tpb\n",
      "shared_mem_size = (tpb, tpb)\n",
      "griddim = bpg, bpg\n",
      "blockdim = tpb, tpb"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Naive version"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Image(url=\"http://docs.nvidia.com/cuda/cuda-c-programming-guide/graphics/matrix-multiplication-without-shared-memory.png\")"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@numba.cuda.jit(\"void(float32[:,:], float32[:,:], float32[:,:])\")\n",
      "def naive_matrix_mult(A, B, C):\n",
      "    x, y = cuda.grid(2)\n",
      "    if x >= n or y >= n:\n",
      "        return\n",
      "\n",
      "    C[y, x] = 0\n",
      "    for i in range(n):\n",
      "        C[y, x] += A[y, i] * B[i, x]\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Optimized version (shared memory + blocking)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Image(url=\"http://docs.nvidia.com/cuda/cuda-c-programming-guide/graphics/matrix-multiplication-with-shared-memory.png\")"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@numba.cuda.jit(\"void(float32[:,:], float32[:,:], float32[:,:])\")\n",
      "def optimized_matrix_mult(A, B, C):\n",
      "    # Declare shared memory\n",
      "    sA = cuda.shared.array(shape=shared_mem_size, dtype=float32)\n",
      "    sB = cuda.shared.array(shape=shared_mem_size, dtype=float32)\n",
      "    \n",
      "    tx = cuda.threadIdx.x\n",
      "    ty = cuda.threadIdx.y\n",
      "    x, y = cuda.grid(2)\n",
      "\n",
      "    acc = 0\n",
      "    for i in range(bpg):\n",
      "        if x < n and y < n:\n",
      "            # Prefill cache\n",
      "            sA[ty, tx] = A[y, tx + i * tpb]\n",
      "            sB[ty, tx] = B[ty + i * tpb, x]\n",
      "\n",
      "        # Synchronize all threads in the block\n",
      "        cuda.syncthreads()\n",
      "\n",
      "        if x < n and y < n:\n",
      "            # Compute product\n",
      "            for j in range(tpb):\n",
      "                acc += sA[ty, j] * sB[j, tx]\n",
      "\n",
      "        # Wait until all threads finish the computation\n",
      "        cuda.syncthreads()\n",
      "\n",
      "    if x < n and y < n:\n",
      "        C[y, x] = acc"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Prepare data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Prepare data on the CPU\n",
      "A = np.array(np.random.random((n, n)), dtype=np.float32)\n",
      "B = np.array(np.random.random((n, n)), dtype=np.float32)\n",
      "\n",
      "print \"N = %d x %d\" % (n, n)\n",
      "\n",
      "# Prepare data on the GPU\n",
      "dA = cuda.to_device(A)\n",
      "dB = cuda.to_device(B)\n",
      "dC = cuda.device_array_like(A)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Benchmark"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Time the unoptimized version\n",
      "s = timer()\n",
      "naive_matrix_mult[griddim, blockdim](dA, dB, dC)\n",
      "numba.cuda.synchronize()\n",
      "e = timer()\n",
      "unopt_ans = dC.copy_to_host()\n",
      "tcuda_unopt = e - s\n",
      "\n",
      "# Time the optimized version\n",
      "s = timer()\n",
      "optimized_matrix_mult[griddim, blockdim](dA, dB, dC)\n",
      "numba.cuda.synchronize()\n",
      "e = timer()\n",
      "opt_ans = dC.copy_to_host()\n",
      "tcuda_opt = e - s"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "#### Result"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "assert np.allclose(unopt_ans, opt_ans)\n",
      "print \"Without shared memory:\", \"%.2f\" % tcuda_unopt, \"s\"\n",
      "print \"With shared memory:\", \"%.2f\" % tcuda_opt, \"s\""
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# Summary\n",
      "\n",
      "- Numba\n",
      "    - opensource low-level GPU support\n",
      "    - CUDA kernel ``@numba.cuda.jit``\n",
      "- NumbaPro\n",
      "    - commerical with high-level GPU support\n",
      "    - vectorize ``@numbapro.vectorize``\n",
      "    - guvectorize ``@numbapro.guvectorize``\n",
      "    - CUDA libraries ``@numbapro.cudalib.*``\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "# Future of Numba / NumbaPro\n",
      "\n",
      "- OpenCL support in Numba (WiP)\n",
      "- Deferred array object that hides most CPU+GPU details\n",
      "    - use it like NumPy array"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}