sequential_strassen_matrix_multiplication_algorithm.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "sequential-strassen-matrix-multiplication-algorithm.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyN1w3DUM++fEAELUGDCDHyL",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU",
    "gpuClass": "standard"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/AlirezaAzadbakht/sequential-strassen-matrix-multiplication-algorithm/blob/main/sequential_strassen_matrix_multiplication_algorithm.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ypvzx6lTtawx"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import argparse\n",
        "import os\n",
        "\n",
        "base_code = \"\"\"\n",
        "import numpy as np\n",
        "# Numba supports CUDA GPU programming by directly compiling a restricted subset of Python code into CUDA kernels and device functions following the CUDA execution model.\n",
        "# so we use numba to access cuda \n",
        "from numba import cuda, types, float32\n",
        "\n",
        "# just in time compiler is a low-level entry point to the CUDA features in Numba\n",
        "@cuda.jit\n",
        "def cuda_matrix_multiplication(A, B, C):\n",
        "    # each thread in each block is the number of columns in matrices\n",
        "    thread_per_block = 4\n",
        "    # define shared arrays with type of float32 for operations\n",
        "    a = cuda.shared.array(shape=(thread_per_block, thread_per_block), dtype=float32)\n",
        "    b = cuda.shared.array(shape=(thread_per_block, thread_per_block), dtype=float32)\n",
        "    # getting block corrdiantes\n",
        "    x, y = cuda.grid(2)\n",
        "    a = A\n",
        "    b = B\n",
        "    cuda.syncthreads()\n",
        "    tx = cuda.threadIdx.x\n",
        "    ty = cuda.threadIdx.y\n",
        "    block_per_grid = cuda.gridDim.x    # blocks per grid\n",
        "\n",
        "    # check if return matrix shapes are correct otherwise return\n",
        "    if x >= C.shape[0] and y >= C.shape[1]:\n",
        "        return\n",
        "\"\"\"\n",
        "\n",
        "def init_matrix(M):\n",
        "    a = np.array(range(M**2)).astype('object')\n",
        "    for i in range(M**2):\n",
        "        a[i] = f'a[{i}]'\n",
        "    a = a.reshape(M,M)\n",
        "\n",
        "    b = np.array(range(M**2)).astype('object')\n",
        "    for i in range(M**2):\n",
        "        b[i] = f'b[{i}]'\n",
        "    b = b.reshape(M,M)\n",
        "    return a, b\n",
        "\n",
        "def code_generator(M, LOOK_UP):\n",
        "    try:\n",
        "        os.remove('generated_code.py')\n",
        "    except:\n",
        "        pass\n",
        "    first = True\n",
        "    with open(\"generated_code.py\", \"a\") as myfile:\n",
        "        myfile.write(f\"{base_code}\\n\")\n",
        "\n",
        "    for i in range(M):\n",
        "        for j in range(M):\n",
        "            if first:\n",
        "                out = f'    if x == {i} and y == {j}:C[x, y] = {LOOK_UP[i,j]}'\n",
        "                first = False\n",
        "            else:\n",
        "                out = f'    elif x == {i} and y == {j}:C[x, y] = {LOOK_UP[i,j]}'\n",
        "            with open(\"generated_code.py\", \"a\") as myfile:\n",
        "                    myfile.write(f\"{out}\\n\")\n",
        "\n",
        "def matrix_multiplication(A, B):\n",
        "    C = np.zeros((A.shape[0], B.shape[1])).astype('object')\n",
        "    for i in range(A.shape[0]):\n",
        "        for j in range(B.shape[1]):\n",
        "            C[i][j] = ''\n",
        "            for k in range(A.shape[0]):\n",
        "                if C[i][j] == '':\n",
        "                    C[i][j] =  f'({A[i][k]} * {B[k][j]})'\n",
        "                else:\n",
        "                    C[i][j] =  f'({C[i][j]} + ({A[i][k]} * {B[k][j]}))'\n",
        "    return C\n",
        "\n",
        "def add(A, B):\n",
        "    C = np.zeros((A.shape[0], A.shape[1])).astype('object')\n",
        "    for i in range(A.shape[0]):\n",
        "        for j in range(A.shape[1]):\n",
        "            C[i][j] =  f'({A[i][j]} + {B[i][j]})'\n",
        "    return C\n",
        "\n",
        "def sub(A, B):\n",
        "    C = np.zeros((A.shape[0], A.shape[1])).astype('object')\n",
        "    for i in range(A.shape[0]):\n",
        "        for j in range(A.shape[1]):\n",
        "            C[i][j] =  f'({A[i][j]} - {B[i][j]})'\n",
        "    return C\n",
        "\n",
        "def strassen(x, y):\n",
        "    if len(x) == 1:\n",
        "        return matrix_multiplication(x, y)\n",
        " \n",
        "    a, b, c, d = x[:x.shape[0]//2, :x.shape[1]//2], x[:x.shape[0]//2, x.shape[1]//2:], x[x.shape[0]//2:, :x.shape[1]//2], x[x.shape[0]//2:, x.shape[1]//2:]\n",
        "    e, f, g, h = y[:y.shape[0]//2, :y.shape[1]//2], y[:y.shape[0]//2, y.shape[1]//2:], y[y.shape[0]//2:, :y.shape[1]//2], y[y.shape[0]//2:, y.shape[1]//2:]\n",
        " \n",
        "    p1 = strassen(a, sub(f, h)) \n",
        "    p2 = strassen(add(a, b), h)       \n",
        "    p3 = strassen(add(c, d), e)       \n",
        "    p4 = strassen(d, sub(g, e))       \n",
        "    p5 = strassen(add(a, d), add(e, h))       \n",
        "    p6 = strassen(sub(b, d), add(g, h)) \n",
        "    p7 = strassen(sub(a, c), add(e, f)) \n",
        "\n",
        "    c11 = add(sub(add(p5 , p4), p2) , p6 )\n",
        "    c12 = add(p1 , p2)          \n",
        "    c21 = add(p3 , p4 )          \n",
        "    c22 = sub(sub(add(p1 , p5 ), p3) , p7) \n",
        "\n",
        "    c = np.vstack((np.hstack((c11, c12)), np.hstack((c21, c22))))\n",
        "    return c\n",
        "\n",
        "def main(M):\n",
        "    a, b = init_matrix(M)\n",
        "    LOOK_UP = strassen(a, b)\n",
        "    code_generator(M, LOOK_UP)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "main(4)"
      ],
      "metadata": {
        "id": "GtQMSkQ3toKW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "from numba import cuda, types, float32\n",
        "import numpy as np\n",
        "from generated_code import cuda_matrix_multiplication\n",
        "\n",
        "M = 4\n",
        "N = 4\n",
        "\n",
        "# prepare 2 random matrix of size 32*32 and one zeros matrix of size 32*32 for results\n",
        "a = np.random.rand(M,N).flatten()\n",
        "b = np.random.rand(N,M).flatten()\n",
        "c = np.zeros((M, M))\n",
        "\n",
        "# convert host arrays to device arrays\n",
        "d_a = cuda.to_device(a)\n",
        "d_b = cuda.to_device(b)\n",
        "d_c = cuda.to_device(c)\n",
        "\n",
        "# set size of blocks and grid here we have one grid and 32*32 block size \n",
        "block_size = (N,N)\n",
        "grid_size = (int(M/N),int(M/N))\n",
        "\n",
        "# calling cuda funtion with grid and block size and device arrays\n",
        "cuda_matrix_multiplication[grid_size, block_size](d_a, d_b, d_c)\n",
        "\n",
        "# copy back result to host arrays\n",
        "c = d_c.copy_to_host()\n",
        "c"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_vmOhFOgt4CQ",
        "outputId": "a9b67171-8e5d-4810-e89f-0a8b575f5819"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[1.75945613, 1.31676909, 1.31827759, 0.90438576],\n",
              "       [1.11057839, 0.9906213 , 0.47307758, 0.43834963],\n",
              "       [1.4850823 , 1.20046921, 1.2415892 , 0.84346731],\n",
              "       [1.06511736, 0.72204581, 0.82190742, 0.65158696]])"
            ]
          },
          "metadata": {},
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        ""
      ],
      "metadata": {
        "id": "4CLfAJg-zRMq"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}