Convolution exploiting GPUs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convolution exploiting GPUs\n",
    "\n",
    "GPUs offer a highly parallel substrate for accelerating data-parallel computations. Most modern GPU architectures are organized as clusters, which are analogous to **Compute Unit** in OpenCL. Each cluster contains several ALUs, called **Processing Elements** in OpenCL-speak, and a **Local Memory** that is available for shared communication within the cluster. Multiple threads of work, or **work items** can be mapped to a cluster depending on the workgroup scheduling. A GPU matches the OpenCL hardware model to the closest extent among existing OpenCL-compatible hardware platforms.\n",
    "\n",
    "1. Local/Constant Memory Optimizations\n",
    "Apart from SIMD, Threading, and Unrolling optimizations, the GPU substrate provides a unique optimization opportunity in the form of memory selection. Both **constant** and **local** memories can be used effectively to help accelerate your OpenCL task. It requires writing OpenCL kernels in a specific way to exploit these memories. This optimization is typically useless on a CPU as neither of these memories are explicitly available, and default to the processor cache. \n",
    "- To use **constant** memory, we simply need to tag the relevant data structures with the **__constant** identifier. The memory allocation on the host also needs to be tagged with **CL_MEM_READ_ONLY** qualifier to help copy the data from the host to the correct RAM on the OpenCL device.\n",
    "- To use **local** memory, we can only do so by declaring a fixed-size array within the OpenCL kernel body and declaring it with the qualifier **__local**. Also, we have to explicitly copy the data from the **__global** memory to the **__local** memory structures ourselves. Depending on the size of the workgroup and the **__local** structures, we can divvy up the memory loading task across multiple work-items. \n",
    "\n",
    "For 2D convolution, we can store the kernels in **__constant** memory, and prefetch portions of the input image into **__local** memory. Alternatively, when considering 3D convolution tasks (multiple 2D convolutions), we can store the output image in **__local** memory instead. We can visualize the memory hierarchy in OpenCL below.\n",
    "\n",
    "![](memory-hierarchy.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we initialize the system with PyOpenCL headers as shown below. Make sure to choose the GPU platform choice when prompted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Choose platform:\n",
      "[0] <pyopencl.Platform 'Intel(R) OpenCL' at 0x2691b10>\n",
      "[1] <pyopencl.Platform 'Clover' at 0x7fde7cdc8ae0>\n",
      "Choice [0]:0\n",
      "Set the environment variable PYOPENCL_CTX='0' to avoid being asked again.\n"
     ]
    }
   ],
   "source": [
    "# Plots\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = (5, 5)        # medium images\n",
    "plt.rcParams['image.interpolation'] = 'nearest'  # don't interpolate: show square pixels\n",
    "plt.rcParams['image.cmap'] = 'gray'  # use grayscale output rather than a (potentially misleading) color heatmap\n",
    "\n",
    "#load the ipython extensions for pyopencl\n",
    "%load_ext pyopencl.ipython_ext\n",
    "\n",
    "from __future__ import absolute_import, print_function\n",
    "import pyopencl as cl\n",
    "\n",
    "ctx = cl.create_some_context(interactive=True)\n",
    "queue = cl.CommandQueue(ctx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now show the OpenCL kernel for 2D convolution optimized for using **constant** and **local** memories below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python2.7/dist-packages/pyopencl/__init__.py:63: CompilerWarning: Non-empty compiler output encountered. Set the environment variable PYOPENCL_COMPILER_OUTPUT=1 to see more.\n",
      "  \"to see more.\", CompilerWarning)\n"
     ]
    }
   ],
   "source": [
    "%%cl_kernel\n",
    "\n",
    "__kernel void convolve2D_gpu(\n",
    "        const __global float *in,               // W*H input images\n",
    "        __constant float *filt,                 // K*K filter kernel\n",
    "        __global float *out)                    // W*H output images\n",
    "{\n",
    "    // get image resolution\n",
    "    const int W = get_global_size(0);\n",
    "    const int H = get_global_size(1);\n",
    "    const int K = 3;\n",
    "\n",
    "    // get pixel position\n",
    "    const int x = get_global_id(0); \n",
    "    const int y = get_global_id(1);\n",
    "\n",
    "    // allocate local RAM for storing input pixels\n",
    "    // local RAM sizes must be statically determined.\n",
    "    __local float in_local[512];\n",
    "\n",
    "    // load data into the local RAM\n",
    "    in_local[x*W+y] = in[x*W+y];\n",
    "    barrier(CLK_LOCAL_MEM_FENCE);\n",
    "\n",
    "    float sum = 0;\n",
    "    int r = 0, c = 0;\n",
    "\n",
    "    // loop over rows\n",
    "    for (r = 0; r < K; r++)\n",
    "    {\n",
    "        // loop over columns\n",
    "        for(c = 0; c < K; c++)\n",
    "        {\n",
    "            sum += filt[r*K+c]*in_local[((y+r)*W+x)+c];\n",
    "        }\n",
    "    }\n",
    "    out[y*W+x] = sum;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this block of code, we tag the **filt** structure as a **constant** and let OpenCL load the data into constant memory with the host code shown below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "        d_filter = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, K*K, h_filter, &err);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The other portion of code that has been added is the allocation, loading, and use of **in_local** variable. By preloading the image into the GPU local memory, accesses to the pixels are now fast. This is important as each pixel in convolution is read $K$x$K$ times. While caching can help reduce overheads of repetitive accesses, some accesses may miss the cache resulting in wasted cycles. We extract the newly added portion of code below:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "        // allocate local RAM for storing input pixels\n",
    "        __local in_local[W*H];\n",
    "\n",
    "        // load data into the local RAM\n",
    "        in_local[x*W+y] = in[x*W+y];\n",
    "        barrier(CLK_LOCAL_MEM_FENCE);```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **barrier** call forces all workitems in the workgroup to synchornize before proceeding to the computation. This is important as all pixels must complete their memory loads before we can use them to do filtering. By enforcing (1) barrier, and (2) memory fence, both synchronization and completion of memory loads is guaranteed. Remember, this only operates on workitems inside a single workgroup on a single compute unit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We the run the OpenCL kernel on the GPU as shown below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import ndimage, misc, signal;\n",
    "import numpy as np\n",
    "\n",
    "f = misc.ascent();\n",
    "in_np = f.astype(np.float32);\n",
    "\n",
    "filt_np = np.array([[1,2,1],[2,4,2],[1,2,1]], dtype=np.float32);\n",
    "filt_np /= sum(filt_np)\n",
    "\n",
    "mf = cl.mem_flags\n",
    "in_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=in_np)\n",
    "filt_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=filt_np)\n",
    "out_g = cl.Buffer(ctx, mf.WRITE_ONLY, in_np.nbytes)\n",
    "\n",
    "# Run OpenCL convolve2D function\n",
    "print(in_np.shape);\n",
    "k_event = convolve2D_gpu(queue, in_np.shape, None, in_g, filt_g, out_g)\n",
    "k_event.wait()\n",
    "\n",
    "out_np = np.zeros_like(in_np, dtype=np.float32)\n",
    "cl.enqueue_copy(queue, out_np, out_g)\n",
    "\n",
    "fig, (a, b) = plt.subplots(1, 2)\n",
    "a.imshow(in_np, cmap='gray')\n",
    "b.imshow(out_np, cmap='gray')\n",
    "plt.show() # apparently, plt.show() works than fig.show() in juypter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}