pi_monte_carlo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "662e62cf-497e-4fc8-bffb-9b2272b75b75",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "06bbb65e-eead-424c-ac73-e7666d5ae0be",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import cuda\n",
    "import pycuda\n",
    "import pycuda.driver as cuda\n",
    "from pycuda.compiler import SourceModule\n",
    "import pycuda.autoinit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "993c47d5-91e3-4129-893a-549064bc57c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_pi_cpu(n_points):\n",
    "    #10 000 000 points gave 3.141932\n",
    "    #100 000 000 points gave 3.14141676\n",
    "    #n_total = 1000000\n",
    "\n",
    "    #First, generate random points\n",
    "    x_rand = np.random.rand(n_points)\n",
    "    y_rand = np.random.rand(n_points)\n",
    "\n",
    "    #Compute radius from origin\n",
    "    inside = np.sqrt(x_rand**2+y_rand**2) <= 1.0\n",
    "    #Count number of points inside\n",
    "    n_inside = np.sum(inside)\n",
    "    \n",
    "    #n_inside = 0\n",
    "    #for i in range(n_points):\n",
    "    #    n_inside += np.sqrt(x_rand[i]**2+y_rand[i]**2) <= 1.0\n",
    "\n",
    "    #We can estimate pi by the following formula:\n",
    "    #pi = 4 * n_inside / n_total\n",
    "    pi = 4*n_inside/n_points\n",
    "    \n",
    "    return pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5bbe45c6-eb3a-4179-a1b0-8c5a6328b35e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.1408359375\n",
      "Time taken: 0.025424 seconds\n"
     ]
    }
   ],
   "source": [
    "tic = time.time()\n",
    "print(compute_pi_cpu(512000))\n",
    "toc = time.time() \n",
    "\n",
    "#for loop: 1.84 seconds for 512 000 elements\n",
    "#vectorized: 0.018 seconds for 512 000 elements => 100 times faster!!!\n",
    "\n",
    "print(\"Time taken: {:f} seconds\".format(toc-tic))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "75fce2ec-ed1a-4650-8ef4-473c38d9933f",
   "metadata": {},
   "outputs": [],
   "source": [
    "pi_kernel_src = \"\"\"\n",
    "//Based on Stroustrup, adapted for CUDA\n",
    "//pseudorandom numbers\n",
    "__device__ float generateRandomNumber(long& last_draw) {\n",
    "    last_draw = last_draw*1103515245 + 12345;\n",
    "    long abs = last_draw & 0x7fffffff;\n",
    "    return abs / 2147483648.0; \n",
    "}\n",
    "\n",
    "\n",
    "__global__ void computePi(unsigned int* inside, unsigned int num_iterations, unsigned int seed) {\n",
    "    __shared__ unsigned int inside_shared[512];\n",
    "    \n",
    "    unsigned int tid = threadIdx.x;\n",
    "    unsigned int bid = blockIdx.x;\n",
    "\n",
    "    //1 generate random numbers\n",
    "    unsigned int num_inside = 0;\n",
    "    for (int i=0; i<num_iterations; ++i) {\n",
    "        long rand_seed = seed + blockIdx.x*blockDim.x + threadIdx.x;\n",
    "        float x = generateRandomNumber(rand_seed);\n",
    "        float y = generateRandomNumber(rand_seed);\n",
    "\n",
    "        //2 compute radius from origin\n",
    "        float r = sqrt(x*x+y*y);\n",
    "\n",
    "        //3 check if inside circle and write to memory\n",
    "        if (r <= 1) {\n",
    "            num_inside += 1;\n",
    "        }\n",
    "    }\n",
    "    inside_shared[tid] = num_inside;\n",
    "    \n",
    "    /////////////////////////\n",
    "    //Shared memory reduction\n",
    "    /////////////////////////\n",
    "    \n",
    "    // Synchronze so that all thread see the same shared memory\n",
    "    __syncthreads();\n",
    "    \n",
    "    // Find the sum in shared memory\n",
    "    //Reduce from 512 to 256 elements\n",
    "    if (threadIdx.x < 256) {\n",
    "        inside_shared[threadIdx.x] = inside_shared[threadIdx.x] + inside_shared[threadIdx.x + 256];\n",
    "    }\n",
    "    __syncthreads();\n",
    "    \n",
    "    //Reduce from 256 to 128 elements\n",
    "    if (threadIdx.x < 128) {\n",
    "        inside_shared[threadIdx.x] = inside_shared[threadIdx.x] + inside_shared[threadIdx.x + 128];\n",
    "    }\n",
    "    __syncthreads();\n",
    "    \n",
    "    //Reduce from 128 to 64 elements\n",
    "    if (threadIdx.x < 64) {\n",
    "        inside_shared[threadIdx.x] = inside_shared[threadIdx.x] + inside_shared[threadIdx.x + 64];\n",
    "    }\n",
    "    __syncthreads();\n",
    "    \n",
    "    //Reduce from 32 to 16 elements\n",
    "    //Since we here have only one active warp (threadIdx.x > 32)\n",
    "    //we do not need to call syncthreads anymore\n",
    "    volatile unsigned int* p = &inside_shared[0]; //To help the compiler not cache this variable...\n",
    "    if (threadIdx.x < 32) {\n",
    "        p[threadIdx.x] = p[threadIdx.x] + p[threadIdx.x + 32];\n",
    "        p[threadIdx.x] = p[threadIdx.x] + p[threadIdx.x + 16];\n",
    "        p[threadIdx.x] = p[threadIdx.x] + p[threadIdx.x + 8];\n",
    "        p[threadIdx.x] = p[threadIdx.x] + p[threadIdx.x + 4];\n",
    "        p[threadIdx.x] = p[threadIdx.x] + p[threadIdx.x + 2];\n",
    "        p[threadIdx.x] = p[threadIdx.x] + p[threadIdx.x + 1];\n",
    "    }\n",
    "    \n",
    "    // Finally write out to output\n",
    "    // NOTE: We have 512 threads, but only thread 0 writes to memory\n",
    "    if (threadIdx.x == 0) {\n",
    "        inside[bid] = p[0];\n",
    "    }\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "mod = SourceModule(pi_kernel_src)\n",
    "compute_pi_gpu_kernel = mod.get_function(\"computePi\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "11992b59-f7de-4f0a-98f7-272fe9cf6d25",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pycuda._driver.Function at 0x7f7893e49180>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#prepare a function call: give pycuda instructions\n",
    "#on what type of parameters to send to CUDA\n",
    "#__global__ void computePi(\n",
    "# unsigned int* inside, <= P\n",
    "# unsigned int num_iterations, <= i\n",
    "# unsigned int seed <= i\n",
    "# );\n",
    "compute_pi_gpu_kernel.prepare(\"Pii\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "936395a1-3289-4b9b-8610-6db971f196d6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.1416140625\n",
      "Time taken: 2.166616 seconds\n"
     ]
    }
   ],
   "source": [
    "def compute_pi_gpu(n_points, iterations_per_thread=1000, threads_per_block=512):\n",
    "    #10 000 000 points gave 3.141932\n",
    "    #100 000 000 points gave 3.14141676\n",
    "    #n_total = 1000000\n",
    "    \n",
    "    assert(n_points % (threads_per_block*iterations_per_thread) == 0)\n",
    "\n",
    "    #Allocate output data on the GPU\n",
    "    #Bytes per unsigned int: \n",
    "    bytes_per_uint = 4\n",
    "    inside_gpu = cuda.mem_alloc(bytes_per_uint*(n_points//threads_per_block)//iterations_per_thread)\n",
    "    \n",
    "    #Execute the pi-kernel\n",
    "    num_blocks = (n_points // threads_per_block)//iterations_per_thread\n",
    "    block=(threads_per_block,1,1)\n",
    "    grid=(num_blocks,1,1)\n",
    "    #compute_pi_gpu_kernel(inside_gpu, np.uint32(iterations_per_thread), np.uint32(time.time()), block=(threads_per_block,1,1), grid=(num_blocks,1,1))\n",
    "    compute_pi_gpu_kernel.prepared_call(grid, block, inside_gpu, np.uint32(iterations_per_thread), np.uint32(time.time()))\n",
    "    \n",
    "    #Allocate memory to download to on the CPU\n",
    "    inside_cpu = np.empty((n_points//threads_per_block)//iterations_per_thread, dtype=np.uint32)\n",
    "    \n",
    "    #Download from the GPU to the CPU\n",
    "    cuda.memcpy_dtoh(inside_cpu, inside_gpu)\n",
    "    \n",
    "    #Count number of points inside\n",
    "    # Version 6: move this reduction to the GPU, and only transfer a single number to the CPU.\n",
    "    n_inside = np.sum(inside_cpu)\n",
    "\n",
    "    #We can estimate pi by the following formula:\n",
    "    pi = 4*n_inside/n_points\n",
    "    \n",
    "    return pi\n",
    "\n",
    "#CPU version: about 18 seconds for 512 000 000 points\n",
    "#version 1: maximum 1024 threads (one block with 1024 threads)\n",
    "#version 2: Slow with 512 000 000 threads, 1.5 seconds (n blocks with 512 threads each)\n",
    "#version 3: less memory. 512 000 000 in 0.029076 seconds (shared memory reduction). About factor 50 speedup\n",
    "#version 4: more work per thread (less memory).  512 000 000 in 0.003880 seconds (for-loop in kernel). About factor 10 speedup (or more!)\n",
    "#           uses 2.16 seconds for 51 200 000 000 000 darts\n",
    "#version 5: prepared call kernel launch??? => no benefit right now. We run the kernel only once!\n",
    "#version 6: sum/reduction on the GPU???\n",
    "#version 7: ???\n",
    "#version 8: ???\n",
    "#profile driven development. \n",
    "\n",
    "tic = time.time()\n",
    "print(compute_pi_gpu(51200000000000, 10000000, 512))\n",
    "toc = time.time()\n",
    "\n",
    "print(\"Time taken: {:f} seconds\".format(toc-tic))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ed16a63-e7fd-48fe-a1fe-3d001d835c7c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}