80_gpu_aos_vs_soa_mem_layout.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/ggruszczynski/gpu_colab/blob/main/80_gpu_aos_vs_soa_mem_layout.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5Tkq6AdIwtOR"
      },
      "source": [
        "# Memory layout - porting issues\n",
        "\n",
        "Consider a set of points (x,y,z) describing some geometry...\n",
        "\n",
        "![Mesh-of-a-F1-car.jpg](https://github.com/ggruszczynski/gpu_colab/blob/main/lectures/figures/Mesh-of-a-F1-car.jpg?raw=1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "80fYwvaihxY0",
        "outputId": "a7d29557-7eb0-46bd-871f-50d8974cbe6b"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Thu Jun 29 12:43:36 2023       \n",
            "+-----------------------------------------------------------------------------+\n",
            "| NVIDIA-SMI 525.85.12    Driver Version: 525.85.12    CUDA Version: 12.0     |\n",
            "|-------------------------------+----------------------+----------------------+\n",
            "| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |\n",
            "| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |\n",
            "|                               |                      |               MIG M. |\n",
            "|===============================+======================+======================|\n",
            "|   0  Tesla T4            Off  | 00000000:00:04.0 Off |                    0 |\n",
            "| N/A   60C    P8    10W /  70W |      0MiB / 15360MiB |      0%      Default |\n",
            "|                               |                      |                  N/A |\n",
            "+-------------------------------+----------------------+----------------------+\n",
            "                                                                               \n",
            "+-----------------------------------------------------------------------------+\n",
            "| Processes:                                                                  |\n",
            "|  GPU   GI   CI        PID   Type   Process name                  GPU Memory |\n",
            "|        ID   ID                                                   Usage      |\n",
            "|=============================================================================|\n",
            "|  No running processes found                                                 |\n",
            "+-----------------------------------------------------------------------------+\n"
          ]
        }
      ],
      "source": [
        "!nvidia-smi"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "CCNRnP5xh26B",
        "outputId": "6a17eb11-8fc0-40e3-e7aa-1fb69bd21f1b"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Total number of elements: 201326592 \n",
            "Dimensions of the matrix MxM: 8192x8192 \n",
            "Memory size of array in [MB]: 1610.612736\n"
          ]
        }
      ],
      "source": [
        "from numba import cuda\n",
        "from numba import jit\n",
        "import numpy as np\n",
        "import math\n",
        "from numba import vectorize, int8, int32, int64, float32, float64\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "np.set_printoptions(precision=1, suppress=True)\n",
        "\n",
        "%matplotlib inline\n",
        "\n",
        "M = 2**13     # ~1000MB, 2**13=8192 on GPU SOA is better than AOS\n",
        "# M = 2**11       # ~100MB, 2**11=2048\n",
        "# M = 2**10     # ~10MB, 2**10=1024 on both CPU & GPU: SOA ~ AOS\n",
        "# M = 2**6      # <1MB, 2**6=64 no difference between AOS and SOA, CPU is faster\n",
        "\n",
        "\n",
        "\n",
        "N = 3*M**2\n",
        "\n",
        "a = np.arange(N, dtype=np.float64) # [0...N] on the host\n",
        "print(f\"Total number of elements: {N} \\nDimensions of the matrix MxM: {M}x{M} \\nMemory size of array in [MB]: {a.nbytes/1E6}\")\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "rvacfqFtkxvn",
        "outputId": "666c4c41-68cf-4fb7-d06c-c6153c3e355e"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Blockspergrid: (256, 256), Total number of threads in a block: 1024\n"
          ]
        }
      ],
      "source": [
        "# threadsperblock = (8,8)  # for M = 2**6 = 64\n",
        "threadsperblock = (32,32)\n",
        "blockspergrid_x = math.ceil(M / threadsperblock[0])\n",
        "blockspergrid_y = math.ceil(M / threadsperblock[1])\n",
        "# blockspergrid_z = math.ceil(d_a.shape[2] / threadsperblock[2])\n",
        "blockspergrid = (blockspergrid_x, blockspergrid_y)\n",
        "total_threads_in_block = threadsperblock[0]*threadsperblock[1]\n",
        "assert total_threads_in_block <= 1024 # hardware limit\n",
        "assert total_threads_in_block <= M # avoid useless threads\n",
        "\n",
        "print(f\"Blockspergrid: {blockspergrid}, Total number of threads in a block: {total_threads_in_block}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [],
      "metadata": {
        "id": "2mLP_01j6lkk"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Uln1RQRirPjA"
      },
      "source": [
        "# AOS (Array of Structures)\n",
        "\n",
        "The MxM matrix stores Points(x,y,z) in its most inner dimension\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "3VQD-8LmiBCs",
        "outputId": "9183495c-4f55-4120-9649-0bac756afdcf"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            " Shape: (8192, 8192, 3) \n",
            " Strides (196608, 24, 8) \n",
            " Content [[0. 1. 2.]\n",
            " [3. 4. 5.]\n",
            " [6. 7. 8.]]\n"
          ]
        }
      ],
      "source": [
        "arr_aos = a.reshape(M,M,3)        # The MxM matrix stores Points(x,y,z) in its most inner dimension\n",
        "aos_out = np.zeros_like(arr_aos)\n",
        "print(f\" Shape: {arr_aos.shape} \\n Strides {arr_aos.strides} \\n Content {arr_aos[0,0:3,:]}\")\n",
        "\n",
        "d_arr_aos = cuda.to_device(arr_aos)           # Copy array from host to the device\n",
        "d_aos_out = cuda.device_array_like(d_arr_aos) # preallocate an arracy filled with 0\n",
        "\n",
        "# print(f\"GPU result:\\n {d_aos_out.copy_to_host()[0,0:3,:]}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "hbO3ERQ8izdI"
      },
      "outputs": [],
      "source": [
        "@jit(nopython=True)\n",
        "def cpu_aos(M, A_in,  A_out):\n",
        "    for tidy in range(M):\n",
        "      for tidx in range(M):\n",
        "        # A_out[tidy][tidx][0] = A_in[tidy][tidx][0]\n",
        "        # A_out[tidy][tidx][1] = A_in[tidy][tidx][1]\n",
        "        A_out[tidy][tidx][2] = A_in[tidy][tidx][2] + 1000\n",
        "\n",
        "@cuda.jit\n",
        "def kernel_gpu_aos(M, A_in, A_out):\n",
        "    tidx, tidy = cuda.grid(2)\n",
        "    # The above is equivalent to the following 2 lines of code:\n",
        "    # x = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.x\n",
        "    # y = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y\n",
        "\n",
        "    # loop over all points in domain (except boundary)\n",
        "    if (tidx < M and tidy < M):\n",
        "        # A_out[tidx][tidy][0] = A_in[tidx][tidy][0]\n",
        "        # A_out[tidx][tidy][1] = A_in[tidx][tidy][1]\n",
        "        A_out[tidx][tidy][2] = A_in[tidx][tidy][2] + 1000\n",
        "\n",
        "@cuda.jit\n",
        "def kernel_gpu_aos2(M, A_in, A_out):\n",
        "    tidx, tidy = cuda.grid(2)\n",
        "    # The above is equivalent to the following 2 lines of code:\n",
        "    # x = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.x\n",
        "    # y = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y\n",
        "\n",
        "    # loop over all points in domain (except boundary)\n",
        "    if (tidx < M and tidy < M):\n",
        "        # A_out[tidy][tidx][0] = A_in[tidy][tidx][0]\n",
        "        # A_out[tidy][tidx][1] = A_in[tidy][tidx][1]\n",
        "        A_out[tidy][tidx][2] = A_in[tidy][tidx][2] + 1001"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "9_hI52-pkCZM",
        "outputId": "1cc128ca-06e0-4bd1-cf12-fc4499a8d168"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "25 ms ± 1.82 ms per loop (mean ± std. dev. of 7 runs, 5000 loops each)\n"
          ]
        }
      ],
      "source": [
        "%timeit -r 7 -n 5000 kernel_gpu_aos[blockspergrid, threadsperblock](M, d_arr_aos, d_aos_out)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%timeit -r 7 -n 5000 kernel_gpu_aos2[blockspergrid, threadsperblock](M, d_arr_aos, d_aos_out)"
      ],
      "metadata": {
        "id": "VuBSLTNj9QsD",
        "outputId": "d8efbfba-597b-4acc-9e7f-949f3e38e0f0",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "20.3 ms ± 398 µs per loop (mean ± std. dev. of 7 runs, 5000 loops each)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "x8htVlgVmFis"
      },
      "outputs": [],
      "source": [
        "# %timeit -r 7 -n 1000 cpu_aos(M, arr_aos, aos_out) #to slow"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "PDd3pDZUp1yC",
        "outputId": "5707bbe7-be31-433f-97ce-86728c2fed69"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "GPU result:\n",
            " [[   0.    0. 1003.]\n",
            " [   0.    0. 1006.]\n",
            " [   0.    0. 1009.]]\n"
          ]
        }
      ],
      "source": [
        "# print(f\"CPU result:\\n {aos_out[0,0:3,:]}\")\n",
        "print(f\"GPU result:\\n {d_aos_out.copy_to_host()[0,0:3,:]}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pDTTeOBArUEQ"
      },
      "source": [
        "## Cleanup"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gcPBw2whrV-e"
      },
      "outputs": [],
      "source": [
        "del arr_aos\n",
        "del aos_out\n",
        "\n",
        "del d_arr_aos\n",
        "del d_aos_out"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wEC6wPRLrKC0"
      },
      "source": [
        "# SOA (Structure of Arrays)\n",
        "\n",
        "Each of the MxM matrices stores only the x, y or z coordinates."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "gat73fHNo3U4",
        "outputId": "e4cca31a-cf0e-4747-cb12-d7857a275466"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            " Shape: (3, 8192, 8192) \n",
            " Strides (536870912, 65536, 8) \n",
            " Content [1.3e+08 1.3e+08 1.3e+08 1.3e+08 1.3e+08 1.3e+08 1.3e+08 1.3e+08 1.3e+08\n",
            " 1.3e+08]\n"
          ]
        }
      ],
      "source": [
        "arr_soa = a.reshape(3,M,M)                     # Each of the most outer dimensions stores a MxM matrix stores only x, y or z coordinates.\n",
        "out_soa = np.zeros_like(arr_soa)\n",
        "\n",
        "print(f\" Shape: {arr_soa.shape} \\n Strides {arr_soa.strides} \\n Content {arr_soa[2,0,:10]}\")\n",
        "\n",
        "d_arr_soa = cuda.to_device(arr_soa)            # Copy array from host to the device\n",
        "d_soa_out = cuda.device_array_like(d_arr_soa)  # preallocate an arracy filled with 0"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "L8pMe3CHmPMR"
      },
      "outputs": [],
      "source": [
        "@jit(nopython=True)\n",
        "def cpu_soa(M, A_in,  A_out):\n",
        "    for tidy in range(M):\n",
        "      for tidx in range(M):\n",
        "        # A_out[0][tidy][tidx] = A_in[0][tidy][tidx]\n",
        "        # A_out[1][tidy][tidx] = A_in[1][tidy][tidx]\n",
        "        A_out[2][tidy][tidx] = A_in[2][tidy][tidx] + 1000\n",
        "\n",
        "@cuda.jit\n",
        "def kernel_gpu_soa(M, A_in,  A_out):\n",
        "    tidx, tidy = cuda.grid(2)\n",
        "    # The above is equivalent to the following 2 lines of code:\n",
        "    # x = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.x\n",
        "    # y = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y\n",
        "\n",
        "    # loop over all points in domain (except boundary)\n",
        "    if (tidx < M and tidy < M):\n",
        "      #  A_out[0][tidx][tidy] = A_in[0][tidx][tidy]\n",
        "      #  A_out[1][tidx][tidy] = A_in[1][tidx][tidy]\n",
        "       A_out[2][tidx][tidy] = A_in[2][tidx][tidy] + 1000\n",
        "\n",
        "@cuda.jit\n",
        "def kernel_gpu_soa2(M, A_in,  A_out):\n",
        "    tidx, tidy = cuda.grid(2)\n",
        "    # The above is equivalent to the following 2 lines of code:\n",
        "    # x = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.x\n",
        "    # y = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y\n",
        "\n",
        "    # loop over all points in domain (except boundary)\n",
        "    if (tidx < M and tidy < M):\n",
        "      #  A_out[0][tidy][tidx] = A_in[0][tidy][tidx]\n",
        "      #  A_out[1][tidy][tidx] = A_in[1][tidy][tidx]\n",
        "       A_out[2][tidy][tidx] = A_in[2][tidy][tidx] + 1001"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jbCeO19HpRdk",
        "outputId": "00d093d5-1b4c-466b-af95-54383716d42d"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "9.59 ms ± 711 µs per loop (mean ± std. dev. of 7 runs, 5000 loops each)\n"
          ]
        }
      ],
      "source": [
        "%timeit -r 7 -n 5000 kernel_gpu_soa[blockspergrid, threadsperblock](M, d_arr_soa, d_soa_out)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%timeit -r 7 -n 5000 kernel_gpu_soa2[blockspergrid, threadsperblock](M, d_arr_soa, d_soa_out)"
      ],
      "metadata": {
        "id": "o1RhqJCD6ti8",
        "outputId": "3d045405-1042-48a9-bc09-8ab737e6a041",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "4.6 ms ± 387 µs per loop (mean ± std. dev. of 7 runs, 5000 loops each)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QfIKS9onl-Cj"
      },
      "outputs": [],
      "source": [
        "# %timeit -r 7 -n 1000 loops_no cpu_soa(M, arr_soa, out_soa)  #to slow"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_iJhdOa5l-fF",
        "outputId": "0a43debf-3c04-45a4-98f8-dd037e6a6be0"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "GPU result:\n",
            " [1.3e+08 1.3e+08 1.3e+08 1.3e+08 1.3e+08 1.3e+08 1.3e+08 1.3e+08 1.3e+08\n",
            " 1.3e+08]\n"
          ]
        }
      ],
      "source": [
        "# print(f\"CPU result:\\n {out_soa[2,0,:10]}\")\n",
        "print(f\"GPU result:\\n {d_soa_out.copy_to_host()[2,0,:10]}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m40dapYIy4nN"
      },
      "source": [
        "## Cleanup"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Q-efEasAywsE"
      },
      "outputs": [],
      "source": [
        "del arr_soa\n",
        "del out_soa\n",
        "\n",
        "del d_arr_soa\n",
        "del d_soa_out"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## timeit: General advise regarding the values of repeat and number\n",
        "\n",
        "A computer has different \"clocks\" to measure times. These clocks have different \"ticks\" (depending on the OS). For example it could measure seconds, milliseconds or nanoseconds - these ticks are called the granularity of the clock.\n",
        "\n",
        "If the duration of the execution is smaller or roughly equal to the granularity of the clock one cannot get representative timings.\n",
        "\n",
        "The `number` (`-n`) and `repeat` (`-r`) are separate arguments because they serve different purposes. The number controls how many executions (aka loops) are done for each timing and it's used to get representative timings. The repeat argument controls how many timings are done and its use is to get accurate statistics. IPython uses the mean or average to calculate the run-time of the statement of all repetitions and then divides that number with number. So it measures the average of the averages.\n",
        "\n",
        "If you want to modify either `number` or `repeat` then you should set number to the minimum value possible without running into the granularity of the timer. In my experience number should be set so that number executions of the function take at least 10 microseconds (0.00001 seconds) otherwise you might only \"time\" the minimum resolution of the \"timer\".\n",
        "\n",
        "The repeat should be set as high as possible. Having more repeats will make it more likely that you really find the real best or average. However more repeats will take longer so there's a trade-off as well.\n",
        "\n",
        "Source: <https://stackoverflow.com/questions/48258008/n-and-r-arguments-to-ipythons-timeit-magic>\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "PG57MdH3Wj0N"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Exercise\n",
        "\n",
        "Try to `timeit` the kernels without specifying the `-r` and `-n` flags explicitly."
      ],
      "metadata": {
        "id": "_tiMrYwlqFZz"
      }
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "provenance": [],
      "gpuType": "T4",
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "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.10"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}