09_Numba.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wath is Numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numba is a library that allows translating Python functions to optimized machine code at runtime using the industry-standard LLVM compiler library. The complete and original guide can be faund at [https://nyu-cds.github.io/python-numba/](https://nyu-cds.github.io/python-numba/)\n",
    "\n",
    "Numba provides the ability to speed up applications with high-performance functions written directly in Python, rather than using language extensions such as Cython.\n",
    "\n",
    "With a few simple annotations, array-oriented and math-heavy Python code can be just-in-time (JIT) optimized to achieve performance similar to C, C++, and Fortran, without having to switch languages or Python interpreters.\n",
    "\n",
    "Numba works at the function level. From a function, Numba can generate native code for that function as well as the wrapper code needed to call it directly from Python. This compilation is done on-the-fly and in-memory.\n",
    "\n",
    "Numba’s main features are:\n",
    "\n",
    "+ On-the-fly code generation (at import time or runtime, at the user’s preference)\n",
    "+ Native code generation for the CPU (default) and GPU hardware\n",
    "+ Integration with the Python scientific software stack (thanks to NumPy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numba's central feature is the ```numba.jit()``` decoration. Using this decorator, it is possible to mark a function for optimization by Numba’s JIT compiler. Various invocation modes trigger differing compilation options and behaviours."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numba reads the Python bytecode for a decorated function and combines this with information about the types of the input arguments to the function. It analyzes and optimizes your code, and finally uses the LLVM compiler library to generate a machine code version of your function, tailored to your CPU capabilities. This compiled version is then used every time your function is called."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see Numba in action. The following is a Python implementation of bubblesort for NumPy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import math\n",
    "os.environ[\"NUMBA_ENABLE_CUDASIM\"] = \"1\"\n",
    "from numba import jit, njit, vectorize, cuda\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "original = np.arange(0.0, 10.0, 0.01, dtype='f4')\n",
    "shuffled = original.copy()\n",
    "np.random.shuffle(shuffled)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bubblesort(X):\n",
    "    N = len(X)\n",
    "    for end in range(N, 1, -1):\n",
    "        for i in range(end - 1):\n",
    "            cur = X[i]\n",
    "            if cur > X[i + 1]:\n",
    "                tmp = X[i]\n",
    "                X[i] = X[i + 1]\n",
    "                X[i + 1] = tmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "sorted_arr = shuffled.copy()\n",
    "bubblesort(sorted_arr)\n",
    "print(np.array_equal(sorted_arr, original))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "314 ms ± 5.89 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit sorted_arr[:] = shuffled[:]; bubblesort(sorted_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "@cuda.jit\n",
    "def bubblesort(X):\n",
    "    N = len(X)\n",
    "    for end in range(N, 1, -1):\n",
    "        for i in range(end - 1):\n",
    "            cur = X[i]\n",
    "            if cur > X[i + 1]:\n",
    "                tmp = X[i]\n",
    "                X[i] = X[i + 1]\n",
    "                X[i + 1] = tmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "\nKernel launch configuration was not specified. Use the syntax:\n\nkernel_function[blockspergrid, threadsperblock](arg0, arg1, ..., argn)\n\nSee https://numba.pydata.org/numba-doc/latest/cuda/kernels.html#kernel-invocation for help.\n\n",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-9-247a9a924ee3>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mget_ipython\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'timeit'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'sorted_arr[:] = shuffled[:]; bubblesort(sorted_arr)'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\IPython\\core\\interactiveshell.py\u001b[0m in \u001b[0;36mrun_line_magic\u001b[1;34m(self, magic_name, line, _stack_depth)\u001b[0m\n\u001b[0;32m   2325\u001b[0m                 \u001b[0mkwargs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'local_ns'\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_local_scope\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstack_depth\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2326\u001b[0m             \u001b[1;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 2327\u001b[1;33m                 \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   2328\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2329\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<decorator-gen-54>\u001b[0m in \u001b[0;36mtimeit\u001b[1;34m(self, line, cell, local_ns)\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\IPython\\core\\magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[1;34m(f, *a, **k)\u001b[0m\n\u001b[0;32m    185\u001b[0m     \u001b[1;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    186\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0marg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 187\u001b[1;33m         \u001b[0mcall\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    188\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    189\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0marg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\IPython\\core\\magics\\execution.py\u001b[0m in \u001b[0;36mtimeit\u001b[1;34m(self, line, cell, local_ns)\u001b[0m\n\u001b[0;32m   1167\u001b[0m             \u001b[1;32mfor\u001b[0m \u001b[0mindex\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m10\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1168\u001b[0m                 \u001b[0mnumber\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m10\u001b[0m \u001b[1;33m**\u001b[0m \u001b[0mindex\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1169\u001b[1;33m                 \u001b[0mtime_number\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtimer\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtimeit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnumber\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1170\u001b[0m                 \u001b[1;32mif\u001b[0m \u001b[0mtime_number\u001b[0m \u001b[1;33m>=\u001b[0m \u001b[1;36m0.2\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1171\u001b[0m                     \u001b[1;32mbreak\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\IPython\\core\\magics\\execution.py\u001b[0m in \u001b[0;36mtimeit\u001b[1;34m(self, number)\u001b[0m\n\u001b[0;32m    167\u001b[0m         \u001b[0mgc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdisable\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    168\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 169\u001b[1;33m             \u001b[0mtiming\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minner\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtimer\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    170\u001b[0m         \u001b[1;32mfinally\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    171\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mgcold\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<magic-timeit>\u001b[0m in \u001b[0;36minner\u001b[1;34m(_it, _timer)\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\numba\\cuda\\simulator\\kernel.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, *args)\u001b[0m\n\u001b[0;32m     66\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     67\u001b[0m         \u001b[1;31m# Ensure we've been given a valid grid configuration\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 68\u001b[1;33m         grid_dim, block_dim = normalize_kernel_dimensions(self.grid_dim,\n\u001b[0m\u001b[0;32m     69\u001b[0m                                                           self.block_dim)\n\u001b[0;32m     70\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\numba\\cuda\\errors.py\u001b[0m in \u001b[0;36mnormalize_kernel_dimensions\u001b[1;34m(griddim, blockdim)\u001b[0m\n\u001b[0;32m     47\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     48\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[1;32mNone\u001b[0m \u001b[1;32min\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mgriddim\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mblockdim\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 49\u001b[1;33m         \u001b[1;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmissing_launch_config_msg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     50\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     51\u001b[0m     \u001b[0mgriddim\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcheck_dim\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgriddim\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'griddim'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: \nKernel launch configuration was not specified. Use the syntax:\n\nkernel_function[blockspergrid, threadsperblock](arg0, arg1, ..., argn)\n\nSee https://numba.pydata.org/numba-doc/latest/cuda/kernels.html#kernel-invocation for help.\n\n"
     ]
    }
   ],
   "source": [
    "%timeit sorted_arr[:] = shuffled[:]; bubblesort(sorted_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to specify the signature of the Numba function. A function signature describes the types of the arguments and the return type of the function. This can produce slightly faster code as the compiler does not need to infer the types. However the function is no longer able to accept other types."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the example below, ```float64(int32, int32)``` is the function's signature specifying a function that takes two 32-bit integer arguments and returns a double precision float. Numba provides a shorthand notation, so the same signature can be specified as ```f8(i4, i4)```.\n",
    "\n",
    "The specialization will be compiled by the ```@jit decorator```, and no other specialization will be allowed. This is useful if you want fine-grained control over types chosen by the compiler (for example, to use single-precision floats).\n",
    "\n",
    "If you omit the return type, e.g. by writing ```(int32, int32)``` instead of ```float64(int32, int32)```, Numba will try to infer it for you. Function signatures can also be strings, and you can pass several of them as a list; see the numba.jit() and cuda.jit() documentation for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import jit, int32, float64\n",
    "\n",
    "@jit(float64(int32, int32))\n",
    "def f(x, y):\n",
    "    # A somewhat trivial example\n",
    "    return (x + y) / 3.14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2738853503184713"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(1, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numba has two compilation modes: ```nopython``` mode and ```object``` mode \n",
    "\n",
    "In ```nopython``` mode, the Numba compiler will generate code that does not access the Python C API. This mode produces the highest performance code, but requires that the native types of all values in the function can be inferred.\n",
    "\n",
    "In ```object``` mode, the Numba compiler generates code that handles all values as Python objects and uses the Python C API to perform all operations on those objects. Code compiled in object mode will often run no faster than Python interpreted code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import jit\n",
    "@jit(\"void(f4[:])\",nopython=True) #without simulator you can use cuda.jit\n",
    "def bubblesort(X):\n",
    "    N = len(X)\n",
    "    for end in range(N, 1, -1):\n",
    "        for i in range(end - 1):\n",
    "            cur = X[i]\n",
    "            if cur > X[i + 1]:\n",
    "                tmp = X[i]\n",
    "                X[i] = X[i + 1]\n",
    "                X[i + 1] = tmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "739 µs ± 15.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit sorted_arr[:] = shuffled[:]; bubblesort(sorted_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numba's ```@vectorize``` decorator allows Python functions taking scalar input arguments to be used as NumPy ufuncs. Creating a traditional ```NumPy ufunc``` is not the most straightforward process and involves writing some C code. Numba makes this easy. Using the ```@vectorize``` decorator, Numba can compile a pure Python function into a ```ufunc``` that operates over NumPy arrays as fast as traditional ufuncs written in C.\n",
    "\n",
    "The @vectorize decorator has two modes of operation:\n",
    "\n",
    "+ *Eager* or decoration-time, compilation. If you pass one or more type signatures to the decorator, you will be building a Numpy ufunc. We’re just going to consider eager compilation here.\n",
    "+ *Lazy*, or call-time, compilation. When not given any signatures, the decorator will give you a Numba dynamic universal function (DUFunc) that dynamically compiles a new kernel when called with a previously unsupported input type.\n",
    "\n",
    "Using ```@vectorize```, you write your function as operating over input scalars, rather than arrays. Numba will generate the surrounding loop (or kernel) allowing efficient iteration over the actual inputs. \n",
    "\n",
    "The following code defines a function that takes two ```float32``` arrays and returns a ```float32``` array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.9575646  1.6118802 ]\n",
      " [1.3659555  0.49743554]\n",
      " [0.7874614  1.5384891 ]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from numba import vectorize, int64, float32\n",
    "\n",
    "@vectorize([float32(float32, float32)], target='parallel', fastmath=True) \n",
    "#to use cuda target must become cuda. Since emulator has a bug we use parallel\n",
    "def vec_add(x, y):\n",
    "    c =  x + y\n",
    "    return c\n",
    "\n",
    "a = np.random.rand(3,2).astype(np.float32) \n",
    "b = np.random.rand(3,2).astype(np.float32) \n",
    "\n",
    "\n",
    "print(vec_add(a, b))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Writing CUDA kernels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "CUDA has an execution model unlike the traditional sequential model used for programming CPUs. In CUDA, the code you write will be executed by multiple threads at once (often hundreds or thousands). Your solution will be modeled by defining a thread hierarchy of grid, blocks, and threads.\n",
    "\n",
    "Numba also exposes three kinds of GPU memory:\n",
    "\n",
    "+ global device memory\n",
    "+ shared memory\n",
    "+ local memory\n",
    "For all but the simplest algorithms, it is important that you carefully consider how to use and access memory in order to minimize bandwidth requirements and contention.\n",
    "\n",
    "NVIDIA recommends that programmers focus on following those recommendations to achieve the best performance:\n",
    "\n",
    "- Find ways to parallelize sequential code\n",
    "- Minimize data transfers between the host and the device\n",
    "- Adjust kernel launch configuration to maximize device utilization\n",
    "- Ensure global memory accesses are coalesced\n",
    "- Minimize redundant accesses to global memory whenever possible\n",
    "- Avoid different execution paths within the same warp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To execute a function on GPU, you have to either define something called a *kernel function* or a *device function*\n",
    "+ kernels explicitly declare their thread hierarchy when called, i.e. the number of blocks and number of threads per block. You can compile your kernel once, and call it multiple times with different block and grid sizes.\n",
    "+ kernels cannot return a value. So, either you will have to do changes on original array, or pass another array for storing the result. For computing scalar, you will have to pass a 1 element array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "@cuda.jit\n",
    "def increment_kernel(io_array):\n",
    "    pos = cuda.grid(1)\n",
    "    # x, y = cuda.grid(2) # For 2D array\n",
    "    if pos < io_array.size:\n",
    "        io_array[pos] += 1 # do the computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.\n",
      " 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.\n",
      " 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.\n",
      " 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.\n",
      " 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.\n",
      " 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.\n",
      " 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.\n",
      " 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.\n",
      " 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.\n",
      " 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.\n",
      " 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]\n"
     ]
    }
   ],
   "source": [
    "import numpy\n",
    "\n",
    "# Create the data array - usually initialized some other way\n",
    "data = numpy.ones(256)\n",
    "\n",
    "# Set the number of threads in a block\n",
    "threadsperblock = 32 \n",
    "\n",
    "# Calculate the number of thread blocks in the grid\n",
    "blockspergrid = (data.size + (threadsperblock - 1)) // threadsperblock\n",
    "\n",
    "# Now start the kernel\n",
    "increment_kernel[blockspergrid, threadsperblock](data)\n",
    "\n",
    "# Print the result\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two main steps:\n",
    "\n",
    "+ Instantiate the kernel proper, by specifying a number of blocks per grid and a number of threads per block. The product of the two will give the total number of threads launched. Kernel instantiation is done by taking the compiled kernel function and indexing it with a tuple of integers.\n",
    "+ Running the kernel, by passing it the input array (and any separate output arrays if necessary). By default, running a kernel is synchronous: the function returns when the kernel has finished executing and the data is synchronized back."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two-level thread hierarchy is important for the following reasons:\n",
    "\n",
    "On the software side, the block size determines how many threads share a given area of shared memory.\n",
    "On the hardware side, the block size must be large enough for full occupation of execution units; recommendations can be found in the CUDA C Programming Guide.\n",
    "The block size you choose depends on a range of factors, including:\n",
    "\n",
    "+ The size of the data array\n",
    "+ The size of the shared mempory per block (e.g. 64KB)\n",
    "+ The maximum number of threads per block supported by the hardware (e.g. 512 or 1024)\n",
    "+ The maximum number of threads per multiprocessor (MP) (e.g. 2048)\n",
    "+ The maximum number of blocks per MP (e.g. 32)\n",
    "+ The number of threads that can be executed concurrently (a “warp” i.e. 32)\n",
    "+ The execution of threads in a warp has a big effect on the computational throughput. If all threads in a warp are executing the same instruction then they can all be executed in parallel. But if one or more threads is executing a different instruction, the warp has to be split into groups of threads, and these groups execute serially.\n",
    "\n",
    "Rules of thumb for threads per block:\n",
    "\n",
    "+ Should be a round multiple of the warp size (32)\n",
    "+ A good place to start is 128-512 but benchmarking is required to determine the optimal value.\n",
    "+ Each streaming multiprocessor (SP) on the GPU must have enough active warps to achieve maximum throughput. In other words, the blocksize is usually selected to maximize the “occupancy”.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When running a kernel, the kernel function's code is executed by every thread once. It therefore has to know which thread it is in, in order to know which array element(s) it is responsible for. More complex algorithms may define more complex responsibilities, but the underlying principle is the same.\n",
    "\n",
    "To help deal with multi-dimensional arrays, CUDA allows you to specify multi-dimensional blocks and grids. In the example above, you could make blockspergrid and threadsperblock tuples of one, two or three integers. Compared to 1-dimensional declarations of equivalent sizes, this doesn’t change anything to the efficiency or behaviour of generated code, but can help you write your algorithms in a more natural way.\n",
    "\n",
    "One way is for the thread to determines its position in the grid and block and manually compute the corresponding array position:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.\n",
      " 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.\n",
      " 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.\n",
      " 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.\n",
      " 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.\n",
      " 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.\n",
      " 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.\n",
      " 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.\n",
      " 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.\n",
      " 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.\n",
      " 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.]\n"
     ]
    }
   ],
   "source": [
    "@cuda.jit\n",
    "def double_kernel(io_array):\n",
    "    # Thread id in a 1D block\n",
    "    tx = cuda.threadIdx.x\n",
    "    # Block id in a 1D grid\n",
    "    ty = cuda.blockIdx.x\n",
    "    # Block width, i.e. number of threads per block\n",
    "    bw = cuda.blockDim.x\n",
    "    # Compute flattened index inside the array\n",
    "    pos = tx + ty * bw\n",
    "    if pos < io_array.size:  # Check array boundaries\n",
    "        io_array[pos] *= 2 # do the computation\n",
    "        \n",
    "        \n",
    "        \n",
    "# Now start the kernel\n",
    "double_kernel[blockspergrid, threadsperblock](data)\n",
    "\n",
    "# Print the result\n",
    "print(data)        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following special objects are provided by the CUDA backend for the sole purpose of knowing the geometry of the thread hierarchy and the position of the current thread within that geometry:\n",
    "\n",
    "+ numba.cuda.threadIdx - The thread indices in the current thread block. For 1-dimensional blocks, the index (given by the x attribute) is an integer spanning the range from 0 to numba.cuda.blockDim - 1. A similar rule exists for each dimension when more than one dimension is used.\n",
    "+ numba.cuda.blockDim - The shape of the block of threads, as declared when instantiating the kernel. This value is the same for all threads in a given kernel, even if they belong to different blocks (i.e. each block is “full”).\n",
    "+ numba.cuda.blockIdx - The block indices in the grid of threads launched a kernel. For a 1-dimensional grid, the index (given by the x attribute) is an integer spanning the range from 0 to numba.cuda.gridDim - 1. A similar rule exists for each dimension when more than one dimension is used.\n",
    "+ numba.cuda.gridDim - The shape of the grid of blocks, i.e. the total number of blocks launched by this kernel invocation, as declared when instantiating the kernel.\n",
    "\n",
    "These objects can be 1-, 2- or 3-dimensional, depending on how the kernel was invoked. To access the value at each dimension, use the x, y and z attributes of these objects, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-21-25cb615a60af>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     19\u001b[0m \u001b[0mblockspergrid\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mblockspergrid_x\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mblockspergrid_y\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     20\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 21\u001b[1;33m \u001b[0mdouble_kernel_2D\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mblockspergrid\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mthreadsperblock\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     22\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\numba\\cuda\\simulator\\kernel.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, *args)\u001b[0m\n\u001b[0;32m    104\u001b[0m                 \u001b[1;32mfor\u001b[0m \u001b[0mgrid_point\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndindex\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mgrid_dim\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    105\u001b[0m                     \u001b[0mbm\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mBlockManager\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfn\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgrid_dim\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mblock_dim\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 106\u001b[1;33m                     \u001b[0mbm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgrid_point\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0mfake_args\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    107\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    108\u001b[0m             \u001b[1;32mfor\u001b[0m \u001b[0mwb\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mretr\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\numba\\cuda\\simulator\\kernel.py\u001b[0m in \u001b[0;36mrun\u001b[1;34m(self, grid_point, *args)\u001b[0m\n\u001b[0;32m    249\u001b[0m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_f\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    250\u001b[0m             \u001b[0mt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mBlockThread\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtarget\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgrid_point\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mblock_point\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 251\u001b[1;33m             \u001b[0mt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstart\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    252\u001b[0m             \u001b[0mthreads\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0madd\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mt\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    253\u001b[0m             \u001b[0mlivethreads\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0madd\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mt\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\threading.py\u001b[0m in \u001b[0;36mstart\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    855\u001b[0m                 \u001b[1;32mdel\u001b[0m \u001b[0m_limbo\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    856\u001b[0m             \u001b[1;32mraise\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 857\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_started\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    858\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    859\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mrun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\threading.py\u001b[0m in \u001b[0;36mwait\u001b[1;34m(self, timeout)\u001b[0m\n\u001b[0;32m    556\u001b[0m             \u001b[0msignaled\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_flag\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    557\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0msignaled\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 558\u001b[1;33m                 \u001b[0msignaled\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_cond\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    559\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0msignaled\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    560\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\threading.py\u001b[0m in \u001b[0;36mwait\u001b[1;34m(self, timeout)\u001b[0m\n\u001b[0;32m    300\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m    \u001b[1;31m# restore state no matter what (e.g., KeyboardInterrupt)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    301\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mtimeout\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 302\u001b[1;33m                 \u001b[0mwaiter\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    303\u001b[0m                 \u001b[0mgotit\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    304\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "from numba import cuda\n",
    "import numpy\n",
    "import math\n",
    "\n",
    "# CUDA kernel\n",
    "@cuda.jit\n",
    "def double_kernel_2D(io_array):\n",
    "    row, col = cuda.grid(2)\n",
    "    if row < io_array.shape[0] and col < io_array.shape[1]:\n",
    "        io_array[row][col] *= 2 # do the computation\n",
    "\n",
    "# Host code   \n",
    "data = numpy.ones((1024,1024))\n",
    "\n",
    "# Configure the blocks\n",
    "threadsperblock = (64, 16) #ThreadsPerBlock must be a divisor of the data dimension: 1024/ 64 / 16 = 1\n",
    "blockspergrid_x = int(math.ceil(data.shape[0] / threadsperblock[0]))\n",
    "blockspergrid_y = int(math.ceil(data.shape[1] / threadsperblock[1]))\n",
    "blockspergrid = (blockspergrid_x, blockspergrid_y)\n",
    "\n",
    "double_kernel_2D[blockspergrid, threadsperblock](data)\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#http://numba.pydata.org/numba-doc/0.13.4/CUDAJit.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *device function* are function that can be called only inside a kernel function. In other words they can be seen as *auxiliary functions*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.ones(100, dtype=np.float64)\n",
    "b = np.ones(100, dtype=np.float64)\n",
    "\n",
    "@cuda.jit(device=True)\n",
    "def cudaDevice_func(a, b):\n",
    "  for i in range(2):\n",
    "    a = math.exp(a*b)\n",
    "  return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@cuda.jit\n",
    "def cudaKernal_func(a, b, result): # cuda.jit does not return result yet\n",
    "  pos = cuda.grid(1)\n",
    "  if (pos < a.shape[0]) and (pos < b.shape[0]):\n",
    "    result[pos] = cudaDevice_func(a[pos], b[pos])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = np.zeros(100, dtype=np.float64)\n",
    "threadsperblock = 32\n",
    "blockspergrid = (a.shape[0] + (threadsperblock - 1)) // threadsperblock\n",
    "\n",
    "\n",
    "cudaKernal_func[threadsperblock, blockspergrid](a, b, result)\n",
    "\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At least you have to pay attention that CUDA kernel execution is designed to be asynchronous with respect to the host program. \n",
    "This means that the kernel launch (```double_kernel_2D[blockspergrid, threadsperblock](data)```) returns immediately, allowing the CPU to continue executing while the GPU works in the background. \n",
    "Only host<->device memory copies or an explicit synchronization call will force the CPU to wait until previously queued CUDA kernels are complete.\n",
    "\n",
    "When you pass host NumPy arrays to a CUDA kernel, Numba has to synchronize on your behalf, but if you pass device arrays, processing will continue. If you launch multiple kernels in sequence without any synchronization in between, they will be queued up to run sequentially by the driver, which is usually what you want. If you want to run multiple kernels on the GPU in parallel (sometimes a good idea, but beware of race conditions!), take a look at CUDA streams."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import cuda\n",
    "import numpy\n",
    "import math\n",
    "\n",
    "# CUDA kernel\n",
    "@cuda.jit\n",
    "def double_kernel_2D(io_array):\n",
    "    row, col = cuda.grid(2)\n",
    "    if row < io_array.shape[0] and col < io_array.shape[1]:\n",
    "        io_array[row][col] *= 2 # do the computation\n",
    "\n",
    "# Host code   \n",
    "data = numpy.ones((1024,1024))\n",
    "cuda.synchronize()\n",
    "# Configure the blocks\n",
    "threadsperblock = (64, 16) #ThreadsPerBlock must be a divisor of the data dimension: 1024/ 64 / 16 = 1\n",
    "blockspergrid_x = int(math.ceil(data.shape[0] / threadsperblock[0]))\n",
    "blockspergrid_y = int(math.ceil(data.shape[1] / threadsperblock[1]))\n",
    "blockspergrid = (blockspergrid_x, blockspergrid_y)\n",
    "\n",
    "\n",
    "device_data = cuda.to_device(data)\n",
    "\n",
    "double_kernel_2D[blockspergrid, threadsperblock](device_data)\n",
    "cuda.synchronize()\n",
    "print(device_data.copy_to_host())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}