Lecture 5 -- Multithreading.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes about yesterday:\n",
    "- Simon Danish gave a related talk at Curry On\n",
    "  - https://nextjournal.com/sdanisch/curry-on-julia\n",
    "  - https://www.youtube.com/watch?v=DCs0_T9BRp0\n",
    "  - Also NextJournal is cool\n",
    "- Zygote paper about differential programming hit the arxiv yesterday\n",
    "  - https://twitter.com/ChrisRackauckas/status/1151793248626774016\n",
    "  - https://arxiv.org/abs/1907.07587\n",
    "  - https://github.com/MikeInnes/zygote-paper\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multithreading in Julia\n",
    "For best effect use Julia 1.2 or even better 1.3\n",
    "\n",
    "MT support is under heavy development right now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NCPU = 4\n",
      "nthreads() = 1\n"
     ]
    }
   ],
   "source": [
    "NCPU = Sys.CPU_THREADS\n",
    "using Base.Threads\n",
    "@show NCPU\n",
    "@show nthreads();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I = 0.0\n",
      "cpus = 1.0\n",
      "I = 1.0\n",
      "cpus = 2.0\n",
      "I = 2.0\n",
      "cpus = 4.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Installing Julia (1.0 threads) kernelspec in /home/vchuravy/.local/share/jupyter/kernels/julia-(1.0-threads)-1.1\n",
      "└ @ IJulia /home/vchuravy/.julia/packages/IJulia/gI2uA/deps/kspec.jl:72\n",
      "┌ Info: Installing Julia (2.0 threads) kernelspec in /home/vchuravy/.local/share/jupyter/kernels/julia-(2.0-threads)-1.1\n",
      "└ @ IJulia /home/vchuravy/.julia/packages/IJulia/gI2uA/deps/kspec.jl:72\n",
      "┌ Info: Installing Julia (4.0 threads) kernelspec in /home/vchuravy/.local/share/jupyter/kernels/julia-(4.0-threads)-1.1\n",
      "└ @ IJulia /home/vchuravy/.julia/packages/IJulia/gI2uA/deps/kspec.jl:72\n"
     ]
    }
   ],
   "source": [
    "using IJulia\n",
    "for I in 0:1.0:log(2, NCPU)\n",
    "    cpus = 2^round(I)\n",
    "    @show cpus\n",
    "    installkernel(\"Julia ($cpus threads)\", env=Dict(\"JULIA_NUM_THREADS\"=>\"$cpus\"))\n",
    "end\n",
    "# Now restart Jupyter and switch kernel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your hardware (on linux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Architecture:                    x86_64\n",
      "CPU op-mode(s):                  32-bit, 64-bit\n",
      "Byte Order:                      Little Endian\n",
      "Address sizes:                   39 bits physical, 48 bits virtual\n",
      "CPU(s):                          4\n",
      "On-line CPU(s) list:             0-3\n",
      "Thread(s) per core:              2\n",
      "Core(s) per socket:              2\n",
      "Socket(s):                       1\n",
      "NUMA node(s):                    1\n",
      "Vendor ID:                       GenuineIntel\n",
      "CPU family:                      6\n",
      "Model:                           142\n",
      "Model name:                      Intel(R) Core(TM) i7-7660U CPU @ 2.50GHz\n",
      "Stepping:                        9\n",
      "CPU MHz:                         3438.326\n",
      "CPU max MHz:                     4000.0000\n",
      "CPU min MHz:                     400.0000\n",
      "BogoMIPS:                        4993.00\n",
      "Virtualization:                  VT-x\n",
      "L1d cache:                       64 KiB\n",
      "L1i cache:                       64 KiB\n",
      "L2 cache:                        512 KiB\n",
      "L3 cache:                        4 MiB\n",
      "NUMA node0 CPU(s):               0-3\n",
      "Vulnerability L1tf:              Mitigation; PTE Inversion; VMX conditional cache flushes, SMT vulnerable\n",
      "Vulnerability Mds:               Mitigation; Clear CPU buffers; SMT vulnerable\n",
      "Vulnerability Meltdown:          Mitigation; PTI\n",
      "Vulnerability Spec store bypass: Mitigation; Speculative Store Bypass disabled via prctl and seccomp\n",
      "Vulnerability Spectre v1:        Mitigation; __user pointer sanitization\n",
      "Vulnerability Spectre v2:        Mitigation; Full generic retpoline, IBPB conditional, IBRS_FW, STIBP conditional, RSB filling\n",
      "Flags:                           fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf tsc_known_freq pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch cpuid_fault epb invpcid_single pti ssbd ibrs ibpb stibp tpr_shadow vnmi flexpriority ept vpid ept_ad fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt intel_pt xsaveopt xsavec xgetbv1 xsaves dtherm ida arat pln pts hwp hwp_notify hwp_act_window hwp_epp md_clear flush_l1d\n"
     ]
    }
   ],
   "source": [
    ";lscpu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fork-join parallelism\n",
    "Julia threading model is based on a fork-join approach and is still considered experimental.\n",
    "\n",
    "Fork-join describes the control flow that a group of threads undergoes. Execution is forked and a anonymous function is then run across all threads. \n",
    "\n",
    "All threads have to join together and serial execution continues.\n",
    "\n",
    "Special care needs to be taken if the loop body access has side-effects or accesses global state. (This includes IO and random numbers)\n",
    "\n",
    "https://github.com/JuliaLang/julia/pull/32477\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10-element Array{Int64,1}:\n",
       " 1\n",
       " 1\n",
       " 1\n",
       " 1\n",
       " 1\n",
       " 1\n",
       " 1\n",
       " 1\n",
       " 1\n",
       " 1"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = zeros(Int, nthreads()*10)\n",
    "@threads for i in 1:length(a)\n",
    "    a[i] = threadid()\n",
    "end\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "A = zeros(10_000)\n",
    "@threads for i in 1:length(A)\n",
    "    A[i] = rand()\n",
    "end\n",
    "\n",
    "# This might change soon https://github.com/JuliaLang/julia/pull/32407\n",
    "const rngs = [MersenneTwister(rand(UInt64)) for i in 1:nthreads()]\n",
    "\n",
    "@threads for i in 1:length(A)\n",
    "    A[i] = rand(rngs[threadid()])\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Atomics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ff (generic function with 1 method)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ff()\n",
    "    acc = 0\n",
    "    @threads for i in 1:10_000\n",
    "        acc += 1\n",
    "    end\n",
    "    return acc\n",
    "end\n",
    "# What is the result?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "f_ref (generic function with 1 method)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function f_ref()\n",
    "    acc = Ref(0)\n",
    "    @threads for i in 1:10_000\n",
    "        acc[] += 1\n",
    "    end\n",
    "    return acc\n",
    "end\n",
    "# What is the result?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "f_acc (generic function with 1 method)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function f_acc()\n",
    "    acc = Atomic{Int64}(0)\n",
    "    @threads for i in 1:10_000\n",
    "        atomic_add!(acc, 1)\n",
    "    end\n",
    "    acc\n",
    "end\n",
    "#What is acc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Locks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "struct Accumulator{T, L}\n",
    "    x::Base.RefValue{T}\n",
    "    lock::L\n",
    "end\n",
    "\n",
    "Base.lock(a::Accumulator) = lock(a.lock)\n",
    "Base.unlock(a::Accumulator) = unlock(a.lock)\n",
    "Base.getindex(a::Accumulator) = a.x[]\n",
    "Base.setindex!(a::Accumulator, val) = a.x[] = val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "f (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function f(accum)\n",
    "    @threads for i in 1:10_000\n",
    "        lock(accum)\n",
    "        accum[] += 1\n",
    "        unlock(accum)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "acc_m = Accumulator(Ref(0), Mutex())\n",
    "f(acc_m)\n",
    "acc_m[]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "acc_sl = Accumulator(Ref(0), SpinLock())\n",
    "f(acc_sl)\n",
    "acc_sl[]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "acc_rsl = Accumulator(Ref(0), RecursiveSpinLock())\n",
    "f(acc_rsl)\n",
    "acc_rsl[]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Recompiling stale cache file /home/vchuravy/.julia/compiled/v1.1/BenchmarkTools/ZXPQo.ji for BenchmarkTools [6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf]\n",
      "└ @ Base loading.jl:1184\n"
     ]
    }
   ],
   "source": [
    "using BenchmarkTools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  203.075 μs (1 allocation: 32 bytes)\n"
     ]
    }
   ],
   "source": [
    "@btime f($acc_m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  100.056 μs (1 allocation: 32 bytes)\n"
     ]
    }
   ],
   "source": [
    "@btime f($acc_sl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  124.969 μs (1 allocation: 32 bytes)\n"
     ]
    }
   ],
   "source": [
    "@btime f($acc_rsl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  136.899 μs (9491 allocations: 148.31 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "10000"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime ff()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  45.171 μs (2 allocations: 48 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Atomic{Int64}(10000)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime f_acc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  2.641 μs (2 allocations: 48 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Base.RefValue{Int64}(10000)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime f_ref()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Useful trick \n",
    "Benchmarks for this were run on:\n",
    "```\n",
    "> lscpu\n",
    "...\n",
    "Thread(s) per core:  2\n",
    "Core(s) per socket:  2\n",
    "...\n",
    "Model name:          Intel(R) Core(TM) i7-7660U CPU @ 2.50GHz\n",
    "...\n",
    "L1d cache:           32K\n",
    "L1i cache:           32K\n",
    "L2 cache:            256K\n",
    "L3 cache:            4096K\n",
    "...\n",
    "Flags: ... avx2 ...\n",
    "```\n",
    "\n",
    "and\n",
    "```\n",
    "> lscpu\n",
    "...\n",
    "Thread(s) per core:  2\n",
    "Core(s) per socket:  12\n",
    "...\n",
    "Model name:          AMD Ryzen Threadripper 1920X 12-Core Processor\n",
    "...\n",
    "L1d cache:           32K\n",
    "L1i cache:           64K\n",
    "L2 cache:            512K\n",
    "L3 cache:            8192K\n",
    "...\n",
    "Flags: ... avx2 ...\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```julia\n",
    "@threads for id in 1:nthreads()\n",
    "    #each thread does something\n",
    "end\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "threaded_sum (generic function with 1 method)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function threaded_sum(arr)\n",
    "   @assert length(arr) % nthreads() == 0\n",
    "   results = zeros(eltype(arr), nthreads())\n",
    "   @threads for tid in 1:nthreads()\n",
    "       # split work\n",
    "       acc = zero(eltype(arr))\n",
    "       len = div(length(arr), nthreads())\n",
    "       domain = ((tid-1)*len +1):tid*len\n",
    "       @inbounds for i in domain\n",
    "           acc += arr[i]    \n",
    "       end\n",
    "       results[tid] = acc\n",
    "   end\n",
    "   sum(results)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"sl-block is-focused\" data-block-type=\"table\" style=\"height: auto; min-width: 120px; min-height: 30px; width: 471px; left: 489px; top: 140px;\" data-origin-id=\"41b373377b462cab80c0707d4d5def85\"><div class=\"sl-block-content fragment visible\" style=\"z-index: 16;\" data-table-cols=\"3\" data-table-rows=\"8\" data-table-has-header=\"true\" data-fragment-index=\"3\"><table><tbody>\n",
    "<tr>\n",
    "<th style=\"width: 91px;\">NT</th>\n",
    "<th>Skylake</th>\n",
    "<th>AMD TR</th>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">sum</td>\n",
    "<td>514.476 μs</td>\n",
    "<td>430.409 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">1</td>\n",
    "<td>1.578 ms</td>\n",
    "<td>1.206 ms</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">2</td>\n",
    "<td>831.411 μs</td>\n",
    "<td>575.872 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">4</td>\n",
    "<td>417.656 μs</td>\n",
    "<td>294.724 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">6</td>\n",
    "<td>X</td>\n",
    "<td>215.986 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">12</td>\n",
    "<td>X</td>\n",
    "<td>109.536 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">24</td>\n",
    "<td>X</td>\n",
    "<td>57.197 μs</td>\n",
    "</tr>\n",
    "</tbody></table><div class=\"editing-ui sl-table-column-resizer\" data-column-index=\"0\" style=\"left: 91px;\"></div><div class=\"editing-ui sl-table-column-resizer\" data-column-index=\"1\" style=\"left: 280px;\"></div></div></div>\n",
    "\n",
    "If your `@threads` performance with one thread is not as fast as a non `@threads` version something is off..., but yeah for linear scaling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```julia\n",
    "function threaded_sum(arr)\n",
    "       ...\n",
    "       @inbounds for i in domain\n",
    "           acc += arr[i]    \n",
    "       end\n",
    "       ...\n",
    "end\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "threaded_sum2 (generic function with 1 method)"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function threaded_sum2(arr)\n",
    "   @assert length(arr) % nthreads() == 0\n",
    "   results = zeros(eltype(arr), nthreads())\n",
    "   @threads for tid in 1:nthreads()\n",
    "       # split work\n",
    "       acc = zero(eltype(arr))\n",
    "       len = div(length(arr), nthreads())\n",
    "       domain = ((tid-1)*len +1):tid*len\n",
    "       @inbounds @simd for i in domain\n",
    "           acc += arr[i]    \n",
    "       end\n",
    "       results[tid] = acc\n",
    "   end\n",
    "   sum(results)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"sl-block is-focused\" data-block-type=\"table\" style=\"height: auto; min-width: 120px; min-height: 30px; width: 471px; left: 480px; top: 126px;\" data-origin-id=\"6e5b20edb25d4ecc3bac7f5355b40906\"><div class=\"sl-block-content fragment visible\" style=\"z-index: 16;\" data-table-cols=\"3\" data-table-rows=\"8\" data-table-has-header=\"true\" data-fragment-index=\"2\"><table><tbody>\n",
    "<tr>\n",
    "<th style=\"width: 91px;\">NT</th>\n",
    "<th>Skylake</th>\n",
    "<th>AMD TR</th>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">sum</td>\n",
    "<td>514.476 μs</td>\n",
    "<td>430.409 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">1</td>\n",
    "<td>493.384 μs</td>\n",
    "<td>401.755 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">2</td>\n",
    "<td>282.030 μs</td>\n",
    "<td>73.408 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">4</td>\n",
    "<td>230.988 μs</td>\n",
    "<td>37.541 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">6</td>\n",
    "<td>X</td>\n",
    "<td>29.185 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">12</td>\n",
    "<td>X</td>\n",
    "<td>16.491 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">24</td>\n",
    "<td>X</td>\n",
    "<td>17.693 μs</td>\n",
    "</tr>\n",
    "</tbody></table><div class=\"editing-ui sl-table-column-resizer\" data-column-index=\"0\" style=\"left: 91px;\"></div><div class=\"editing-ui sl-table-column-resizer\" data-column-index=\"1\" style=\"left: 280px;\"></div></div></div>\n",
    "\n",
    "We note:\n",
    "1. Hyperthreading\n",
    "2. Superlinear speedup from 1-2 threads on Threadripper\n",
    "  - Cache effect\n",
    "  - Data is 12MB, 2xL3 = 16MB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $3*2^{21} = 48MB$\n",
    "<div class=\"sl-block is-focused\" data-block-type=\"table\" style=\"height: auto; width: 452px; left: 28px; top: 235px;\" data-origin-id=\"bca56d5bb6bb71fa96a5eef8fc66f63b\"><div class=\"sl-block-content fragment visible\" style=\"z-index: 12;\" data-table-cols=\"3\" data-table-rows=\"8\" data-table-has-header=\"true\" data-fragment-index=\"0\"><table><tbody>\n",
    "<tr>\n",
    "<th style=\"width: 76px;\">NT</th>\n",
    "<th style=\"width: 183px;\">Skylake</th>\n",
    "<th>AMD TR</th>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 76px;\">sum</td>\n",
    "<td style=\"width: 183px;\">2.423 ms</td>\n",
    "<td>1.723 ms</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 76px;\">1</td>\n",
    "<td style=\"width: 183px;\">2.295 ms</td>\n",
    "<td>1.610 ms</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 76px;\">2</td>\n",
    "<td style=\"width: 183px;\">1.307 ms</td>\n",
    "<td>1.158 ms</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 76px;\">4</td>\n",
    "<td style=\"width: 183px;\">1.106 ms</td>\n",
    "<td>582.885 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 76px;\">6</td>\n",
    "<td style=\"width: 183px;\">X</td>\n",
    "<td>470.023 μs&nbsp;</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 76px;\">12</td>\n",
    "<td style=\"width: 183px;\">X</td>\n",
    "<td>627.699 μs</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 76px;\">24</td>\n",
    "<td style=\"width: 183px;\">X</td>\n",
    "<td>595.068 μs</td>\n",
    "</tr>\n",
    "</tbody></table><div class=\"editing-ui sl-table-column-resizer\" data-column-index=\"0\" style=\"left: 76px;\"></div><div class=\"editing-ui sl-table-column-resizer\" data-column-index=\"1\" style=\"left: 259px;\"></div></div></div>\n",
    "\n",
    "### $3*2^{22} = 96MB$\n",
    "\n",
    "<div class=\"sl-block is-focused\" data-block-type=\"table\" style=\"height: auto; width: 452px; left: 508px; top: 235px;\" data-origin-id=\"c9b30544172708e13701b826a53c57b4\"><div class=\"sl-block-content fragment visible\" style=\"z-index: 11;\" data-table-cols=\"3\" data-table-rows=\"8\" data-table-has-header=\"true\" data-fragment-index=\"0\"><table><tbody>\n",
    "<tr>\n",
    "<th style=\"width: 91px;\">NT</th>\n",
    "<th>Skylake</th>\n",
    "<th>AMD TR</th>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">sum</td>\n",
    "<td>4.970 ms</td>\n",
    "<td>3.477 ms</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">1</td>\n",
    "<td>4.996 ms</td>\n",
    "<td>3.249 ms</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">2</td>\n",
    "<td>3.867 ms</td>\n",
    "<td>2.241 ms</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">4</td>\n",
    "<td>2.742 ms</td>\n",
    "<td>1.195 ms</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">6</td>\n",
    "<td>X</td>\n",
    "<td>1.143 ms</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">12</td>\n",
    "<td>X</td>\n",
    "<td>1.225 ms</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td style=\"width: 91px;\">24</td>\n",
    "<td>X</td>\n",
    "<td>1.305 ms</td>\n",
    "</tr>\n",
    "</tbody></table><div class=\"editing-ui sl-table-column-resizer\" data-column-index=\"0\" style=\"left: 91px;\"></div><div class=\"editing-ui sl-table-column-resizer\" data-column-index=\"1\" style=\"left: 271px;\"></div></div></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## False sharing\n",
    "\n",
    "https://software.intel.com/en-us/articles/avoiding-and-identifying-false-sharing-among-threads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sum_cacheaware (generic function with 1 method)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function sum_simple()\n",
    "    acc = zeros(Int64, nthreads())\n",
    "    @threads for tid in 1:nthreads()\n",
    "        for i in 1:10_000\n",
    "             acc[tid] += 1\n",
    "        end\n",
    "    end\n",
    "    sum(acc)\n",
    "end\n",
    "\n",
    "function sum_cacheaware()\n",
    "    CACHE_LINE = 64 #bytes\n",
    "    elems = div(CACHE_LINE, sizeof(Int64))\n",
    "    acc = zeros(Int64, nthreads()*elems)\n",
    "    @threads for tid in 1:nthreads()\n",
    "        store = (tid-1)*elems+1\n",
    "        for i in 1:10_000\n",
    "             acc[store] += 1\n",
    "        end\n",
    "    end\n",
    "    sum(acc)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  128 bytes\n",
       "  allocs estimate:  2\n",
       "  --------------\n",
       "  minimum time:     12.630 μs (0.00% GC)\n",
       "  median time:      13.824 μs (0.00% GC)\n",
       "  mean time:        15.444 μs (0.00% GC)\n",
       "  maximum time:     86.649 μs (0.00% GC)\n",
       "  --------------\n",
       "  samples:          10000\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark sum_simple()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  192 bytes\n",
       "  allocs estimate:  2\n",
       "  --------------\n",
       "  minimum time:     14.048 μs (0.00% GC)\n",
       "  median time:      17.409 μs (0.00% GC)\n",
       "  mean time:        18.417 μs (0.00% GC)\n",
       "  maximum time:     87.734 μs (0.00% GC)\n",
       "  --------------\n",
       "  samples:          10000\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark sum_cacheaware()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A fun user submission\n",
    "Help me! this code is slow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  4.32 GiB\n",
       "  allocs estimate:  40000002\n",
       "  --------------\n",
       "  minimum time:     4.174 s (6.70% GC)\n",
       "  median time:      4.299 s (6.27% GC)\n",
       "  mean time:        4.299 s (6.27% GC)\n",
       "  maximum time:     4.425 s (5.86% GC)\n",
       "  --------------\n",
       "  samples:          2\n",
       "  evals/sample:     1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of threads: 1\n"
     ]
    }
   ],
   "source": [
    "using BenchmarkTools\n",
    "using Base.Threads\n",
    "using LinearAlgebra\n",
    "using Random\n",
    "println(\"Number of threads: \", nthreads())\n",
    "\n",
    "function myfun(rng::MersenneTwister)\n",
    "    s = 0.0\n",
    "    N = 10000\n",
    "    for i=1:N\n",
    "        s+=det(randn(rng, 3,3))\n",
    "    end\n",
    "    s/N\n",
    "end\n",
    "\n",
    "rgi   = [MersenneTwister(rand(UInt)) for _ in 1:nthreads()]\n",
    "function bench(rgi)\n",
    "    a  = zeros(1000)\n",
    "    @threads for i=1:length(a)\n",
    "        a[i] = myfun(rgi[threadid()])\n",
    "    end\n",
    "end\n",
    "\n",
    "result = @benchmark bench($rgi)\n",
    "display(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Thread 1\n",
    "\n",
    "```\n",
    "BenchmarkTools.Trial: \n",
    "  memory estimate:  4.32 GiB\n",
    "  allocs estimate:  40000002\n",
    "  --------------\n",
    "  minimum time:     4.063 s (4.03% GC)\n",
    "  median time:      4.217 s (3.57% GC)\n",
    "  mean time:        4.217 s (3.57% GC)\n",
    "  maximum time:     4.371 s (3.14% GC)\n",
    "  --------------\n",
    "  samples:          2\n",
    "  evals/sample:     1\n",
    "```\n",
    "\n",
    "#### Thread 4\n",
    "\n",
    "```\n",
    "BenchmarkTools.Trial: \n",
    "  memory estimate:  3.52 GiB\n",
    "  allocs estimate:  34212074\n",
    "  --------------\n",
    "  minimum time:     3.346 s (0.00% GC)\n",
    "  median time:      3.960 s (10.85% GC)\n",
    "  mean time:        3.960 s (10.85% GC)\n",
    "  maximum time:     4.574 s (18.78% GC)\n",
    "  --------------\n",
    "  samples:          2\n",
    "  evals/sample:     1\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How did I optimise this code?\n",
    "1. Memory allocations in hot-loop\n",
    "2. Eliminate allocs caused by rand\n",
    "3. Investigate how det is implemented\n",
    "4. Implement det!\n",
    "5. Remove overhead to library call\n",
    "6. Use profiling tools\n",
    "7. Start using StaticArrays\n",
    "\n",
    "Full story here: https://hackmd.io/@dLigA9a4SwKmdcaQtloXXw/BkyZ5Mmbb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```julia\n",
    "@edit det(zeros(3, 3)) -> det(lufact(A))\n",
    "lufact(A, pivot = true) = lufact!(copy(A), pivot)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```julia\n",
    "det!(A) = det(lufact!(A))\n",
    "det!(A) = det(LinearAlgebra.generic_lufact!(A))\n",
    "```\n",
    "\n",
    "`det!` originally was calling a `lufact!` from LAPACK,\n",
    "which is overkill for the matrix size. First attempt switch to a pure Julia implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  7.98 KiB\n",
       "  allocs estimate:  2\n",
       "  --------------\n",
       "  minimum time:     553.192 ms (0.00% GC)\n",
       "  median time:      748.534 ms (0.00% GC)\n",
       "  mean time:        741.813 ms (0.00% GC)\n",
       "  maximum time:     910.000 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          7\n",
       "  evals/sample:     1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of threads: 1\n"
     ]
    }
   ],
   "source": [
    "using BenchmarkTools\n",
    "using StaticArrays\n",
    "using Base.Threads\n",
    "println(\"Number of threads: \", nthreads())\n",
    "\n",
    "function myfun(rng::MersenneTwister)\n",
    "    s = 0.0\n",
    "    N = 10000\n",
    "    for i=1:N\n",
    "        s += det(randn(rng, SMatrix{3, 3}))\n",
    "    end\n",
    "    s/N\n",
    "end\n",
    "\n",
    "rgi   = [MersenneTwister(abs(rand(Int))) for s in 1:nthreads()]\n",
    "\n",
    "function bench(rgi)\n",
    "    a  = zeros(1000)\n",
    "    @threads for i=1:length(a)\n",
    "        @inbounds a[i] = myfun(rgi[threadid()])\n",
    "    end\n",
    "end\n",
    "\n",
    "result = @benchmark bench($rgi)\n",
    "display(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The future is near! Partr is coming\n",
    "- Julia 1.2 and 1.3\n",
    "- https://github.com/JuliaLang/julia/pull/32600\n",
    "- https://github.com/JuliaLang/julia/pull/32477\n",
    "- https://github.com/NHDaly/CspExamples.jl/blob/master/src/CspExamples.jl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "@par (macro with 1 method)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "macro par(expr)\n",
    "    thunk = esc(:(()->($expr)))\n",
    "    quote\n",
    "        local task = Task($thunk)\n",
    "        task.sticky = false\n",
    "        schedule(task)\n",
    "        task\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "ename": "ErrorException",
     "evalue": "type Task has no field sticky",
     "output_type": "error",
     "traceback": [
      "type Task has no field sticky",
      "",
      "Stacktrace:",
      " [1] setproperty!(::Task, ::Symbol, ::Bool) at ./sysimg.jl:19",
      " [2] top-level scope at In[17]:5",
      " [3] top-level scope at In[56]:1"
     ]
    }
   ],
   "source": [
    "@par println(\"Hello!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Julia 1.3 task can be executed on multiple worker-threads allowing fine-grained control. This is concurrency ala Go/CSP.\n",
    "\n",
    "Our handy trick from above can then simply be written as:\n",
    "\n",
    "```\n",
    "tasks = Task[]\n",
    "for tid in 1:nthreads()\n",
    "    task = @par begin\n",
    "    ### some work\n",
    "    end\n",
    "    push!(tasks, task)\n",
    "end\n",
    "for task in tasks\n",
    "    wait(task)\n",
    "end\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": "7da284ef7d1b463aa8d5703488914791",
   "lastKernelId": "6462d7ea-92cb-417f-8bfa-26bbcf6bf3c5"
  },
  "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": 2
}