Lecture 2a -- A note on performance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# GEMM in Julia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m  Updating\u001b[22m\u001b[39m registry at `~/.julia/registries/General`\n",
      "\u001b[32m\u001b[1m  Updating\u001b[22m\u001b[39m git-repo `https://github.com/JuliaRegistries/General.git`\n",
      "\u001b[?25l\u001b[2K\u001b[?25h"
     ]
    }
   ],
   "source": [
    "using Pkg\n",
    "Pkg.activate(\"envs/lecture2a\")\n",
    "Pkg.instantiate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "N = 1024\n",
    "A = rand(N, N);\n",
    "B = rand(N, N);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mygemm! (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mygemm!(C, A, B)\n",
    "   for j in 1:size(C, 2)\n",
    "    for k in 1:size(A, 2)\n",
    "      for i in 1:size(C, 1)\n",
    "        @inbounds C[i, j] += A[i, k] * B[k, j]\n",
    "                end \n",
    "    end\n",
    "  end\n",
    "  return C\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mygemm (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mygemm(A, B)\n",
    "  @assert size(A, 2) == size(B, 1)\n",
    "  C = zeros(size(A, 1), size(B,2))\n",
    "  mygemm!(C, A, B)\n",
    "  return C\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mygemm(A, B) ≈ A*B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.419721 seconds (6 allocations: 8.000 MiB)\n"
     ]
    }
   ],
   "source": [
    "# Note matrix size 1024x1024\n",
    "@time mygemm(A, B);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.068962 seconds (6 allocations: 8.000 MiB)\n"
     ]
    }
   ],
   "source": [
    "@time A*B; # Using OpenBLAS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# So Julia is slow?\n",
    "\n",
    "- No. Our implementation is slow!\n",
    "- Today TB showed you all the tricks one needs to do to make it fast (in C)\n",
    "- Those same tricks apply to Julia\n",
    "- Performance of a pure Julia implementation (https://github.com/Sacha0/GemmDemo.jl)\n",
    "- reference is OpenBLAS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Matrix Multiply](gemm.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Is Julia fast?\n",
    "*Enough talk -- let's code*\n",
    "\n",
    "$$\n",
    "sum(a) = \\sum_i^n a_i\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- This material began life as a wonderful [lecture by Steven Johnson at MIT](https://github.com/stevengj/18S096-iap17/blob/master/lecture1/Boxes-and-registers.ipynb).\n",
    "- With apologies to the numerical computing folks, this is not the algorithm you should use!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# pick a large N to not measure call-overhead\n",
    "data = rand(10_000_000);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "C_code =  \"\"\"\n",
    "#include <stddef.h>\n",
    "double c_sum(size_t n, double *X) {\n",
    "    double s = 0.0;\n",
    "    for (size_t i = 0; i < n; ++i) {\n",
    "        s += X[i];\n",
    "    }\n",
    "    return s;\n",
    "}\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "using Libdl\n",
    "const Clib = tempname()   # make a temporary file\n",
    "\n",
    "# compile to a shared library by piping C_code to gcc\n",
    "# (works only if you have gcc installed):\n",
    "\n",
    "open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * \".\" * Libdl.dlext) -`, \"w\") do f\n",
    "    print(f, C_code) \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "c_sum (generic function with 1 method)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function c_sum(X::Array{Float64})\n",
    "    ccall((\"c_sum\", Clib), \n",
    "          Float64,\n",
    "          (Csize_t, Ref{Float64}),\n",
    "          length(X), X)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Precompiling PyCall [438e738f-606a-5dbb-bf0a-cddfbfd45ab0]\n",
      "└ @ Base loading.jl:1186\n"
     ]
    }
   ],
   "source": [
    "using PyCall\n",
    "\n",
    "# Get two python objects that represent data\n",
    "# First a list and then a numpy array\n",
    "# We do this to cut down conversion overhead\n",
    "apy_list = PyCall.array2py(data)\n",
    "apy_numpy = PyObject(data)\n",
    "\n",
    "# get the Python built-in \"sum\" function:\n",
    "pysum = pybuiltin(\"sum\");\n",
    "# get the Numpy \"sum\" function:\n",
    "numpy_sum = pyimport(\"numpy\").sum;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "py\"\"\"\n",
    "def py_sum(a):\n",
    "    s = 0.0\n",
    "    for x in a:\n",
    "        s = s + x\n",
    "    return s\n",
    "\"\"\"\n",
    "\n",
    "sum_py = py\"py_sum\";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mysum (generic function with 1 method)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mysum(data)\n",
    "  acc = zero(eltype(data))\n",
    "  for x in data\n",
    "      acc += x\n",
    "  end\n",
    "  return acc\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "sum(a::<b>AbstractArray</b>) in Base at <a href=\"https://github.com/JuliaLang/julia/tree/55e36cc308b66d3472990a06b2797f9f9154ea0a/base/reducedim.jl#L648\" target=\"_blank\">reducedim.jl:648</a>"
      ],
      "text/plain": [
       "sum(a::AbstractArray) in Base at reducedim.jl:648"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@which sum(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element BenchmarkTools.BenchmarkGroup:\n",
       "  tags: []\n",
       "  \"Julia builtin\" => Benchmark(evals=1, seconds=5.0, samples=10000)\n",
       "  \"Julia handwritten\" => Benchmark(evals=1, seconds=5.0, samples=10000)\n",
       "  ⋮"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using BenchmarkTools\n",
    "\n",
    "suite = BenchmarkGroup()\n",
    "suite[\"Julia handwritten\"]       = @benchmarkable mysum($data)\n",
    "suite[\"Julia builtin\"]           = @benchmarkable sum($data)\n",
    "suite[\"Simple C function\"]       = @benchmarkable c_sum($data)\n",
    "suite[\"Python builtin (list)\"]   = @benchmarkable $pysum($apy_list)\n",
    "suite[\"Python builtin (numpy)\"]  = @benchmarkable $numpy_sum($apy_numpy)\n",
    "suite[\"Python handwritten\"]      = @benchmarkable $sum_py($apy_list)\n",
    "\n",
    "# If a cache of tuned parameters already exists, use it, otherwise, tune and cache\n",
    "# the benchmark parameters. Reusing cached parameters is faster and more reliable\n",
    "# than re-tuning `suite` every time the file is included.\n",
    "paramspath = joinpath(@__DIR__, \"sum_bench.json\")\n",
    "\n",
    "if isfile(paramspath)\n",
    "    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);\n",
    "else\n",
    "    tune!(suite)\n",
    "    BenchmarkTools.save(paramspath, params(suite));\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element BenchmarkTools.BenchmarkGroup:\n",
       "  tags: []\n",
       "  \"Julia builtin\" => Trial(4.329 ms)\n",
       "  \"Julia handwritten\" => Trial(11.348 ms)\n",
       "  ⋮"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = run(suite)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julia builtin.........................4.33 ms\n",
      "Python builtin (numpy)................4.36 ms\n",
      "Simple C function.....................11.3 ms\n",
      "Julia handwritten....................11.35 ms\n",
      "Python builtin (list)................39.39 ms\n",
      "Python handwritten..................223.75 ms\n"
     ]
    }
   ],
   "source": [
    "for (name, trial) in sort(collect(results), by=x->time(x[2]))\n",
    "    t = time(trial) / 1e6\n",
    "    println(rpad(name, 25, \".\"), lpad(string(round(t, digits=2), \" ms\"), 20, \".\"))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "# In conclusion:\n",
    "The point is not that Julia has the fasted sum, the point is that a simple Julia implementation (that looks and feels like Python) is as fast as simple C, and that well optimised Julia code can be even faster.\n",
    "\n",
    "We can improve performance by giving the compiler more information `@inbounds` and `@simd`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Aside: Which method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "sum(a::<b>AbstractArray</b>) in Base at <a href=\"https://github.com/JuliaLang/julia/tree/55e36cc308b66d3472990a06b2797f9f9154ea0a/base/reducedim.jl#L648\" target=\"_blank\">reducedim.jl:648</a>"
      ],
      "text/plain": [
       "sum(a::AbstractArray) in Base at reducedim.jl:648"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@which sum(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "mapreduce(f, op, A::<b>AbstractArray</b>) in Base at <a href=\"https://github.com/JuliaLang/julia/tree/55e36cc308b66d3472990a06b2797f9f9154ea0a/base/reducedim.jl#L304\" target=\"_blank\">reducedim.jl:304</a>"
      ],
      "text/plain": [
       "mapreduce(f, op, A::AbstractArray) in Base at reducedim.jl:304"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@which mapreduce(identity, +, data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "_mapreduce<i>{T}</i>(f, op, ::<b>IndexLinear</b>, A::<b>AbstractArray{T,N} where N</b>) in Base at <a href=\"https://github.com/JuliaLang/julia/tree/55e36cc308b66d3472990a06b2797f9f9154ea0a/base/reduce.jl#L302\" target=\"_blank\">reduce.jl:302</a>"
      ],
      "text/plain": [
       "_mapreduce(f, op, ::IndexLinear, A::AbstractArray{T,N} where N) where T in Base at reduce.jl:302"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@which Base._mapreduce(identity, Base.add_sum, IndexStyle(data), data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "mapreduce_impl(f, op, A::<b>AbstractArray</b>, ifirst::<b>Integer</b>, ilast::<b>Integer</b>) in Base at <a href=\"https://github.com/JuliaLang/julia/tree/55e36cc308b66d3472990a06b2797f9f9154ea0a/base/reduce.jl#L178\" target=\"_blank\">reduce.jl:178</a>"
      ],
      "text/plain": [
       "mapreduce_impl(f, op, A::AbstractArray, ifirst::Integer, ilast::Integer) in Base at reduce.jl:178"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inds = LinearIndices(data)\n",
    "@which Base.mapreduce_impl(identity, Base.add_sum, data, first(inds), last(inds))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blksize = Base.pairwise_blocksize(identity, Base.add_sum) = 1024\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "mapreduce_impl(f, op, A::<b>AbstractArray</b>, ifirst::<b>Integer</b>, ilast::<b>Integer</b>, blksize::<b>Int64</b>) in Base at <a href=\"https://github.com/JuliaLang/julia/tree/55e36cc308b66d3472990a06b2797f9f9154ea0a/base/reduce.jl#L156\" target=\"_blank\">reduce.jl:156</a>"
      ],
      "text/plain": [
       "mapreduce_impl(f, op, A::AbstractArray, ifirst::Integer, ilast::Integer, blksize::Int64) in Base at reduce.jl:156"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@show blksize = Base.pairwise_blocksize(identity, Base.add_sum)\n",
    "@which Base.mapreduce_impl(identity, Base.add_sum, data, first(inds), last(inds),blksize)"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": "b5c7c8eae1dc44599b32bf8c7eebe85f",
   "lastKernelId": "f8bfe3c0-5c79-444c-9d12-5abdcea9bded"
  },
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Julia 1.1.1",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}