notebook_posts/Thread Parallelism in Julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia has 3 kinds of parallelism.\n",
    "The well known, safe, slowish and easyish, *distributed parallelism*, via `pmap`, `@spawn` and `@remotecall`.\n",
    "The wellish known, very safe, very easy, not-actually-parallelism, *asynchronous parallelism* via `@async`.\n",
    "And the more obscure, less documented, experimental, really unsafe, *shared memory parallelism* via `@threads`.\n",
    "It is the last we are going to talk about today.\n",
    "\n",
    "I'm not sure if I can actually teach someone how to write threaded code.\n",
    "Let alone efficient threaded code.\n",
    "But this is me giving it a shot.\n",
    "The example here is going to be fairly complex.\n",
    "For a much simpler example of use,\n",
    "on a problem that is more easily parallelizable,\n",
    "see my recent [stackoverflow post on parallelizing sorting](https://stackoverflow.com/a/47235391/179081).\n",
    "\n",
    "(Spoilers: in the end I don't manage to extract any serious performance gains from parallelizing this prime search. Unlike parallelizing that sorting. Paralising sorting worked out great)\n",
    "\n",
    "<!--more-->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a [previous post](http://white.ucc.asn.au/2017/11/18/Lazy-Sequences-in-Julia.html#primes),\n",
    "I used prime generation as an example to motivate the use of coroutines as generators.\n",
    "Now coroutines are neither parallelism, nor fast.\n",
    "Lets see how fast we can go if we want to crank it up using  `Base.Threading`.\n",
    "(answer: not as much as you might hope).\n",
    "\n",
    "I feel that julia threading is a bit nerfed.\n",
    "In that all threading must take place in a for-loop, where work is distributed equally to all threads.\n",
    "And the loop end blocks til all threads are done.\n",
    "You can't just fire off one thread to do a thing and then let it go.\n",
    "I spent some time a while ago trying to workout how to do that,\n",
    "and in short I found that end of thread block is hard to get around.\n",
    "`@async` on its own can't break out of it.\n",
    "Though one could rewrite ones whole program to never actually exit that loop.\n",
    "But then one ends up building own own threading system.\n",
    "And I have a thesis to finish."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Primes\n",
    "This is the same paragraph from that earlier post.\n",
    "**I'll let you know now, this is not an optimal prime number finding algorithm, by any means.**\n",
    "We're just using it for demonstration. It has a good kind of complexity for talking about shared memory parallelism.\n",
    "\n",
    "\n",
    "If a number is prime, then no prime (except the number itself), will divide it.\n",
    "Since if it has a divisor that is non-prime, then that divisor itself, will have a prime divisor that will divide the whole.\n",
    "So we only need to check primes as candidate divisors.\n",
    "Further: one does not need to check divisibility by all prior primes in order to check if a number $s$ is prime.\n",
    "One only needs to check divisibility by the primes less than or equal to $\\sqrt{x}$, since if $x=a \\times b$, for some $a>\\sqrt{x}$ that would imply that $b<\\sqrt{x}$, and so its composite nature would have been found when $b$ was checked as a divisor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Here is the channel code for before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Base.Test\n",
    "using Primes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "primes_ch (generic function with 2 methods)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function primes_ch(T=BigInt)\n",
    "    known_primes = T[2]\n",
    "    Channel(csize=256, ctype=T) do c\n",
    "        x = big\"3\"\n",
    "        put!(c, 2) # Output the first prime, as we already put int in the list of known primes\n",
    "        while true\n",
    "            for p in known_primes\n",
    "                if p > sqrt(x)\n",
    "                    # Must be prime as we have not found any divisor\n",
    "                    push!(known_primes, x)\n",
    "                    put!(c, x)\n",
    "                    break\n",
    "                end\n",
    "                if x % p == 0 # p divides\n",
    "                    # not prime\n",
    "                    break\n",
    "                end\n",
    "            end\n",
    "            x+=1            \n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.076319 seconds (6.91 M allocations: 258.664 MiB, 20.54% gc time)\n",
      "  1.139354 seconds (6.10 M allocations: 240.307 MiB, 17.69% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time collect(Iterators.take(primes_ch(UInt), 10^4));\n",
    "@time collect(Iterators.take(primes_ch(BigInt), 10^4));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Primes Array\n",
    "\n",
    "So the first and obvious thing to do is to switch to doing this eagerly with an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "primes_array (generic function with 2 methods)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function primes_array(n,T=UInt)\n",
    "    known_primes = T[2]\n",
    "    sizehint!(known_primes, n)\n",
    "    \n",
    "    x=T(3)\n",
    "    while true\n",
    "        for p in known_primes\n",
    "            if p > sqrt(x)\n",
    "                # Must be prime as we have not found any divisor\n",
    "                push!(known_primes, x)\n",
    "                break\n",
    "            end\n",
    "            if x % p == 0 # p divides\n",
    "                # not prime\n",
    "                break\n",
    "            end\n",
    "        end\n",
    "        x+=1\n",
    "        length(known_primes) == n && break\n",
    "    end\n",
    "    return known_primes\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.188296 seconds (1.70 M allocations: 26.084 MiB, 11.46% gc time)\n",
      "  1.356425 seconds (5.89 M allocations: 231.466 MiB, 15.92% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time primes_array(10^4, UInt);\n",
    "@time primes_array(10^4, BigInt);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  3.627753 seconds (26.43 M allocations: 404.349 MiB, 0.73% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time primes_array(10^5, Int);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives an improvement, but not as much as we might really hope for.\n",
    "(as you will see below getting more out of it is harder)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gutting `@threads` for\n",
    "the `@threads` macro eats for-loops and breaks up their ranges equally, one block per thread.\n",
    "That isn't very practical is your plan does not just a processing of some data that doesn't depend strongly on the order of processing.\n",
    "\n",
    "\n",
    "We don't plan on sequentially processing the data, since breaking all numbers into equal blocks, would result the final thread not being able to do anything until almost all the other threads were done.\n",
    "For this algorithm we need to know all the prime numbers less than $\\sqrt{x}$ before we can check if $x$ is prime.\n",
    "So we have a sequential component.\n",
    "\n",
    "So we gut the `@threads` macro, taking the core functionality,\n",
    "and we will manage giving work to the threads ourselves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nthreads() = 10\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "everythread"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Base.Threads\n",
    "@show nthreads()\n",
    "\n",
    "\"\"\"\n",
    "    everythread(fun)\n",
    "Run `fun` on everythread.\n",
    "Returns when every instance of `fun` completes\n",
    "\"\"\"\n",
    "everythread(fun) = ccall(:jl_threading_run, Ref{Void}, (Any,), fun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to check it is working:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bool[true, true, true, true, true, true, true, true, true, true]"
     ]
    }
   ],
   "source": [
    "called = fill(false, nthreads())\n",
    "everythread() do \n",
    "    called[threadid()] = true\n",
    "end\n",
    "called |> showcompact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### push_sorted!\n",
    "Before we can get into actually working on the parallelism, we need another part.\n",
    "\n",
    "Pushing to the end of our list of `known_primes` is no longer going to guarantee order.\n",
    "One thing we will need is the ability to `push!` that does maintain order.\n",
    "Because otherwise we could endup thinking we have checked enough factors but actually we skipped over one.\n",
    "(I made that mistake in an earlier version of this code).\n",
    "\n",
    "We could use a priority queue for this, but since known primes will always be almost sorted,\n",
    "I think it is going to be faster just to insert the elements into a normal vector.\n",
    "Less pointer dereferencing than using a heap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[37mTest Summary: | \u001b[39m\u001b[22m\u001b[1m\u001b[32mPass  \u001b[39m\u001b[22m\u001b[1m\u001b[36mTotal\u001b[39m\u001b[22m\n",
      "push_sorted!  | \u001b[32m   3  \u001b[39m\u001b[36m    3\u001b[39m\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "inserts `xs` at the last position such that all elements prior to it are smaller than it.\n",
    "And all after (if any) are larger.\n",
    "Assumes `xs` is sorted (increasing),\n",
    "\"\"\"\n",
    "function push_sorted!(xs, y)\n",
    "\tfor ii in endof(xs):-1:1\n",
    "\t\t@inbounds prev = xs[ii]\n",
    "\t\tif prev < y\n",
    "\t\t\tinsert!(xs, ii+1, y) # ii+1 is the index in the resulting array for y\n",
    "\t\t\treturn xs\n",
    "\t\tend\n",
    "\tend\n",
    "\t# If got to here then y must be the smallest element, so put at start\n",
    "\tunshift!(xs, y)\n",
    "end\n",
    "\n",
    "using Base.Test\n",
    "@testset \"push_sorted!\" begin\n",
    "    xs = [1,2,4]\n",
    "    @test push_sorted!(xs, 5) == [1,2,4,5]\n",
    "    @test push_sorted!(xs, 3) == [1,2,3,4,5]\n",
    "    @test push_sorted!(xs, 0) == [0,1,2,3,4,5]\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# primes threaded\n",
    "\n",
    "Ok so here we go with the main content of this piece.\n",
    "\n",
    "Here is our plan.\n",
    "\n",
    " - We will have a shared list of all `known_primes` recording what we know.\n",
    " - Each thread will grab a number and check if it is prime, if it is prime it will add it to that list.\n",
    " - To check if it is prime, it needs to check for if there are no primes that divide it.\n",
    " - This means that if there is a prime divisor that is not yet ready it must wait until it is.\n",
    "\n",
    "\n",
    "So what can go wrong?\n",
    "The most important part of getting share memory parallism correct is making sure **at no point is the same piece of memory being both written and read** (or written and written).\n",
    "There is no promise that any operation is actually atomic, except **atomic** operations, and the setting of **locks**.\n",
    "Which brings me to our two tools for dealing with ensuring that memory is not dual operated on.\n",
    "\n",
    "Atomic operators are a small set of operations available on primitive types.\n",
    "They run on atomic types.\n",
    "They might not perform quiet the operation you expect (so [read the docs](https://docs.julialang.org/en/release-0.6/stdlib/parallel/#Base.Threads.Atomic)).\n",
    "For example `atomic_add!(a::Atomic{T}, b::T)`  updates the value of `a`, but returns its old value, as type `T`.\n",
    "Julia's atomics come out of LLVM, more or less directly.\n",
    "\n",
    "Then there are locks.\n",
    "These are what you use if you want to make a block of code (which might modify non-primitively typed memory) not run at the same time as some other block of code.\n",
    "Julia has two kinds of locks `TatasLock`/`SpinLock`, and `Mutex`.\n",
    "We're going to use the first kind, they are based on **atomics**.\n",
    "The second kind (the `Mutex`) is based on lib_uv's OS agnostic wrapper of they OS's locking system.\n",
    "\n",
    "So what do we need to restrict:\n",
    "\n",
    " - `next_check`: the integer that keeps track what is the next. If we let multiple threads read it at the same time then they will initially keep checking the same numbers as each other. Once they get out of sync bad things will happen. Since it is a primitive type (unless a BigInt or similar is passed as the type) we can use an atomic.\n",
    " - `known_primes`: the list of primes we know about. Here are the operations we need to prevent against:\n",
    "    - Reading an element while it is being written (obvious reasons)\n",
    "    - Reading the length of the vector while something is being inserted (may return corrupt overly high value, leading to incorrect state flow in the program, and/or a segfault)\n",
    "    - Reading **any** element while an element is being inserted. This one caught me out, badly, a lot. Even if the element you are reading isn't being touched, it can still fail. The reason for this is that the `Vector` basically reserves (and uses) the right to move itself in memory whenever an element is added, even if you `sizehint!` it. If this occurs in the middle of a `getindex` then the value you think you are reading might not be there any more.\n",
    "    \n",
    "    \n",
    "The other thing we have going on is that we want to sleep our current thread if we are blocked by waiting for a missing prime.\n",
    "This is done using `Condition`, `wait` and `notify`  ([docs](https://docs.julialang.org/en/stable/manual/control-flow/#Tasks-and-events-1)).\n",
    "The advantage of sleeping the thread while it is waiting is that if oversubscribed (or you are doing other things on the computer), any threads currently waiting for a turn on the CPU can get it. I'm not oversubscribing here so it doesn't really matter. If anything it is slowing it down.\n",
    "Still it is good practice, and makes you a polite multi-threading citizen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "primes_threaded (generic function with 2 methods)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function primes_threaded(n,T=UInt128)\n",
    "    known_primes_lock = SpinLock()\n",
    "    prime_added = Condition()\n",
    "    known_primes = Vector{T}()\n",
    "    sizehint!(known_primes, n + nthreads()) #Allocate extra memory inc\n",
    "    push!(known_primes, 2)\n",
    "    \n",
    "    function safe_length_known_primes()\n",
    "        lock(known_primes_lock) do\n",
    "            length(known_primes)\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    function ith_prime(ii) # try and read the ith prime, if it is available. If not theen wait til it is\n",
    "        while(safe_length_known_primes()<ii)\n",
    "            wait(prime_added)\n",
    "        end\n",
    "        local ret\n",
    "        lock(known_primes_lock) do\n",
    "            @inbounds ret = known_primes[ii]\n",
    "            # we need this lock incase it is being reshuffled\n",
    "        end\n",
    "        ret\n",
    "    end\n",
    "    \n",
    "    function add_prime!(p) # Add a prime to our list and let anyone why was waiting for it know \n",
    "        lock(known_primes_lock) do\n",
    "            push_sorted!(known_primes, p)\n",
    "        end\n",
    "        notify(prime_added, all=true)\n",
    "    end\n",
    "\n",
    "    \n",
    "    next_check = Atomic{T}(3) # This is the (potentially prime) number the next thread that asked for something to check will et\n",
    "    everythread() do\n",
    "        while(true)\n",
    "            x=atomic_add!(next_check, T(1)) #atomic_add! returns the *old* value befoe the addition\n",
    "            for ii in 1:x #Not going to get up to this but it will be fine (except at x=2, got to watch that, goot thing we already have 2 covered)\n",
    "                p = ith_prime(ii) \n",
    "                if p > sqrt(x)\n",
    "                    # Must be prime as we have not found any divisor\n",
    "                    add_prime!(x)\n",
    "                    break\n",
    "                end\n",
    "                if x % p == 0 # p divides\n",
    "                    # not prime\n",
    "                    break\n",
    "                end\n",
    "            end\n",
    "\n",
    "            if safe_length_known_primes() >= n\n",
    "                return\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return known_primes\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0-element Array{Int64,1}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ps = primes_threaded(10^5, Int)\n",
    "fails = ps[.!isprime.(ps)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.781142 seconds (3.67 M allocations: 68.560 MiB)\n",
      "  0.475077 seconds (902.99 k allocations: 13.904 MiB, 49.03% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time primes_threaded(10^4, Int);\n",
    "@time primes_array(10^4, Int);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "#### Wait it is slower?\n",
    "That is right,\n",
    "this multi-threaded code is much slower that the array code.\n",
    "Getting performance out of multi-threading is hard.\n",
    "\n",
    "I can't teach it.\n",
    "But I can show you want I am going to work it out next.\n",
    "My theory is that there is too much lock contention."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Reducing Lock Contention\n",
    "\n",
    "Working in blocks reduces contention, it also results in more cache-friendly code.\n",
    "Instead of each thread asking for one number to check then checking it,\n",
    "then asking for another,\n",
    "each thread asks for a bunch to check at a time.\n",
    "The obvious contention reduction is with the atomic `next_check`.\n",
    "The less obvious is in the lock for `known_primes` which is checked ever time one wants to know how long it is to test if it is time to exit the loop.\n",
    "\n",
    "In the code that follows, while each thread asks for a block of numbers to check at a time, it reports found primes individually. I looked at having each thread collect them up localizing in a block and then inserting them into the main-list all at once. But I found that actually slowed things down. It mean allocating a lot more memory, and (I guess) the longer consecutive time in which `known_primes` was locked for the big inserts was problematic.\n",
    "Delaying checks for longer.\n",
    "\n",
    "The really big cause of contention, I feel is the time to read `known_primes`.\n",
    "Especially, the earlier elements.\n",
    "The smaller the prime the more likely it is to be a factor.\n",
    "So we would like to at least be able to check these early primes without worrying about getting locks.\n",
    "To do that we need to maintain a separate list of them.\n",
    "\n",
    "I initially, just wanted to have an atomic value keeping track of up to how far in `known_primes`, was safe to read, without having to worry about the elements changing.\n",
    "Such that everything was in one array; and we knew which were safe to read.\n",
    "But we can't do that, because inserting elements can cause the array to reallocate, so requires it to be locked.\n",
    "So we just use a second array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "primes_avoid (generic function with 3 methods)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function primes_avoid(n,T=UInt128, blocksize=256)\n",
    "    known_primes_lock = SpinLock()\n",
    "    prime_added = Condition()\n",
    "    \n",
    "    pre_known_primes = primes_array(blocksize, T) # most common factors, stored so we don't need to lock to check them\n",
    "    #Unfortunately we need to store a separate list, and can't just have a nonmutating part of the mainlist, as even with sizehinting it sometimes deallocates during an insert\n",
    "    \n",
    "    known_primes = Vector{T}()\n",
    "    # Need to initialize it with enough primes that no block is ever waiting for a prime it itself produces   \n",
    "    sizehint!(known_primes, n + nthreads()*ceil(Int, log(blocksize))) #Allocate extra memory inc\n",
    "    \n",
    "    function safe_length_known_primes()\n",
    "        lock(known_primes_lock) do\n",
    "            length(known_primes)\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    function ith_prime(ii) # try and read the ith prime, if it is available. If not theen wait til it is\n",
    "        local ret\n",
    "        if ii < length(pre_known_primes)\n",
    "            @inbounds ret = pre_known_primes[ii] \n",
    "        else\n",
    "            ii-=length(pre_known_primes) # reindex for main list\n",
    "        \n",
    "            while(safe_length_known_primes()<ii)\n",
    "                wait(prime_added)\n",
    "            end\n",
    "\n",
    "            lock(known_primes_lock) do\n",
    "                @inbounds ret = known_primes[ii] \n",
    "            end\n",
    "        end\n",
    "        ret\n",
    "    end\n",
    "       \n",
    "    function add_prime!(p) # Add a prime to our list and let anyone why was waiting for it know \n",
    "        lock(known_primes_lock) do\n",
    "            push_sorted!(known_primes, p)\n",
    "        end\n",
    "        notify(prime_added, all=true)\n",
    "    end\n",
    "    \n",
    "    next_check = Atomic{T}(blocksize + 1) #already checked the first block during initialisation\n",
    "    everythread() do\n",
    "    #(f->f())() do\n",
    "        while(true)\n",
    "            x_start=atomic_add!(next_check, T(blocksize)) #atomic_add! returns the *old* value befoe the addition\n",
    "            x_end = x_start + blocksize\n",
    "            for x in x_start:x_end\n",
    "            \n",
    "                for ii in 1:x #Not going to get up to this but it will be fine (except at x=2, got to watch that, goot thing we already have 2 covered)\n",
    "                    p = ith_prime(ii) \n",
    "                    if p > sqrt(x)\n",
    "                        # Must be prime as we have not found any divisor\n",
    "                        add_prime!(x)\n",
    "                        break\n",
    "                    end\n",
    "                    if x % p == 0 # p divides\n",
    "                        # not prime\n",
    "                        break\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "            safe_length_known_primes() + length(pre_known_primes) > n && return\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    return append!(pre_known_primes, known_primes)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.235902 seconds (889.55 k allocations: 18.294 MiB)\n",
      "  0.162847 seconds (902.99 k allocations: 13.904 MiB)\n"
     ]
    }
   ],
   "source": [
    "@time primes_avoid(10^4, Int, 256);\n",
    "@time primes_array(10^4, Int);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  6.886030 seconds (20.65 M allocations: 436.509 MiB)\n",
      "  4.135483 seconds (26.43 M allocations: 404.349 MiB)\n"
     ]
    }
   ],
   "source": [
    "gc()\n",
    "@time primes_avoid(10^5, Int, 256);\n",
    "@time primes_array(10^5, Int);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So with that we've manage to scrape in a bit closers, but we are still losing to the single threaded array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flaw in this algorithm when running in parallel\n",
    "There is actually a flaw in this algorithm, I think.\n",
    "Potentially, if your threads are far enough out of sync,\n",
    "one could be waiting for a prime potential factor,\n",
    "and the prime factor that arrives next, is not actually the next prime;\n",
    "and further more that prime arriving early is larger than $\\sqrt{x}$, so terminates the search;\n",
    "incorrectly reporting $x$ as prime.\n",
    "Which if the next prime to arrive was smaller than $\\sqrt{x}$ and was a prime factor of $x$, that would make $x$ not a prime.\n",
    "\n",
    "One solution would be to keep trace of which indices are actually stable.\n",
    "We know an index is stable if it every thread is now working on checking a number that is greater than the prime at that index.\n",
    "\n",
    "Pretty sure it is super unlikely and never happens,\n",
    "but a fix for it gives me an idea for how to go faster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working in blocks seriously\n",
    "\n",
    "Before I said we were working in blocks but we were still pushing everything into a single array at the end.\n",
    "We could actually work in Vector of Vectors,\n",
    "it makes indexing harder but lets us be fine grained with our locks.\n",
    "\n",
    "So what we are going to do is at the start of each block,\n",
    "we are going to reserve a point in out Vector of Vectors known primes,\n",
    "as well as what numbers we are going to check.\n",
    "\n",
    "We need to allocate all the block locations at the start,\n",
    "because increasing the size of an array is not threadsafe.\n",
    "A big complication is we don't know how many blocks we are going to need.\n",
    "It took me a long time to workout the solution to this.\n",
    "What we do is when we run out of allocated memory we let the block of code that is running on all threads terminate,\n",
    "then we allocate more memory and restart it.\n",
    "\n",
    "This code is pretty complex.\n",
    "As you can see from all the assert statements it took me a fair bit of debugging to get it right.\n",
    "Its still not much (if at all) better than serial.\n",
    "But I think it well illustrates how you have to turn problems around to eak out speed when trying to parallelize them.\n",
    "\n",
    "Note in particular how `reserved_blocks` is a vector of atomics indexed by `threadid()` to keep track of what memory is being held by what thread."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "primes_blockmore (generic function with 3 methods)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function primes_blockmore(n,T=UInt128, blocksize=256)\n",
    "    reserved_blocks = [Atomic{Int}(typemax(T)) for _ in 1:nthreads()]\n",
    "    reserved_conds = [Condition() for _ in 1:nthreads()]\n",
    "    \n",
    "    function inc_prime_pointer(block, ind) # try and read the ith prime, if it is available. If not theen wait til it is       \n",
    "        #@assert block<minimum(getindex.(reserved_blocks))\n",
    "        #@assert(block<=safe_length_known_primes[], \"1 rbs= $(getindex.(reserved_blocks)), kp=$(safe_length_known_primes[]), block=$(block)\")           \n",
    "        #@assert(isassigned(known_primes, block), \"block not assigned $(threadid()) $(block)\")\n",
    "        if length(known_primes[block])>ind\n",
    "            (block, ind+1)\n",
    "        else\n",
    "            #time to move to the next block\n",
    "            block += 1\n",
    "            for (owner, rb) in enumerate(reserved_blocks)\n",
    "                while true\n",
    "                    # Check to make sure the block we want to read isn't still reserved.\n",
    "                    if block == rb[]\n",
    "                        # Technically I think I actually need to synconize here,\n",
    "                        # against the lock being released in between me looking at it \n",
    "                        # and me wanting to wait for it's condition.\n",
    "                        @inbounds wait(reserved_conds[owner]) #wait til that block is ready\n",
    "                        break\n",
    "                    end\n",
    "                    length((known_primes[block]))>0 & break # skip empty blocks\n",
    "                end\n",
    "\n",
    "            end\n",
    "            #@assert length(known_primes[block])>0\n",
    "            #@assert block<minimum(getindex.(reserved_blocks))\n",
    "            #@assert(block<=safe_length_known_primes[], \"2 rbs= $(getindex.(reserved_blocks)), kp=$(safe_length_known_primes[]), block=$(block)\")\n",
    "            (block, 1)\n",
    "        end\n",
    "    end\n",
    "      \n",
    "    reserving_block_lock = SpinLock()\n",
    "    next_check = blocksize # Not an atomic as we are already protecting it with a lock\n",
    "    function get_next_block()\n",
    "        reservation = Vector{T}()\n",
    "        lock(reserving_block_lock) do\n",
    "            \n",
    "            cur_len = true_length_known_primes[]\n",
    "            out_of_allocation = cur_len == max_true_length_known_primes\n",
    "            @assert cur_len < max_true_length_known_primes\n",
    "            if out_of_allocation\n",
    "                (true, (reservation, -1,-1)) # could be using a nullable here, but I don't want the pointer\n",
    "            else\n",
    "                atomic_add!(true_length_known_primes, 1)\n",
    "                cur_len += 1\n",
    "                reserved_blocks[threadid()][] = cur_len\n",
    "                @inbounds known_primes[cur_len] = reservation\n",
    "                cur_check = next_check            \n",
    "                next_check+=blocksize\n",
    "                \n",
    "                (false, (reservation, cur_check, cur_check + blocksize))\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    ## Setup initial known_primes\n",
    "    known_primes = Vector{Vector{T}}(1)\n",
    "    max_true_length_known_primes = 1;\n",
    "    @inbounds known_primes[1] = primes_array(blocksize, T)\n",
    "    \n",
    "    @inbounds total_known = Atomic{Int}(length(known_primes[1]))\n",
    "    true_length_known_primes = Atomic{Int}(1) #The number of blocks that are started, exactly\n",
    "    safe_length_known_primes = Atomic{Int}(1) #The number of blocks that are done, a lower bound\n",
    "    \n",
    "    \n",
    "    blocksize÷=0.5\n",
    "    while(total_known[] < n) # This outerloop is to add more memory when we run out of blocks, so everything must sync up\n",
    "        blocksize*=2 # double the size each round as primes are getting rarer.\n",
    "        @show total_known[], next_check[]\n",
    "        flush(STDOUT)\n",
    "        max_true_length_known_primes += n;\n",
    "        # The upper bound on how many blocks we will allow\n",
    "        # unfortunately this is more than 1 block per prime\n",
    "        # without solving for the inverse Prime number theorem it is hard to bound\n",
    "        # have to preallocate it, AFAICT `push!` is never threadsafe.\n",
    "        append!(known_primes, Vector{Vector{T}}(max_true_length_known_primes))\n",
    "        # We are now in a position to reallocated it\n",
    "        \n",
    "        everythread() do\n",
    "        #(f->f())() do #this line is useful instead of everythread for debugging\n",
    "            while(true)\n",
    "                (done, (local_primes, x_start, x_end)) = get_next_block()\n",
    "                done && return # quit now, we are out of allocated memory\n",
    "                for x in x_start:x_end\n",
    "                    pp_block, pp_ind = (1, 1)\n",
    "                    while(true)\n",
    "\n",
    "                        # if we have an index for it we know it is safe to read\n",
    "                        @inbounds p = known_primes[pp_block][pp_ind]\n",
    "                        if p > sqrt(x)\n",
    "                            # Must be prime as we have not found any divisor\n",
    "                            push!(local_primes, x)\n",
    "                            break\n",
    "                        end\n",
    "                        if x % p == 0 # p divides\n",
    "                            # not prime\n",
    "                            break\n",
    "                        end\n",
    "                        pp_block, pp_ind = inc_prime_pointer(pp_block, pp_ind)\n",
    "                    end\n",
    "                end\n",
    "\n",
    "                #End of block stuff.\n",
    "                @inbounds notify(reserved_conds[threadid()], all=true)\n",
    "                atomic_add!(safe_length_known_primes, 1)\n",
    "                @assert length(local_primes) > 0\n",
    "                atomic_add!(total_known, length(local_primes))\n",
    "                total_known[] > n && return\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    known_primes=known_primes[1:true_length_known_primes[]]\n",
    "    all_primes = T[]\n",
    "    sizehint!(all_primes, n)\n",
    "    reduce(append!, T[], known_primes)    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(total_known[], next_check[]) = (256, 256)\n",
      "  0.169644 seconds (847.93 k allocations: 17.723 MiB)\n"
     ]
    }
   ],
   "source": [
    "gc()\n",
    "@time primes_blockmore(10^4, Int, 256);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(total_known[], next_check[]) = (256, 256)\n",
      "  4.106018 seconds (16.94 M allocations: 355.002 MiB)\n"
     ]
    }
   ],
   "source": [
    "gc()\n",
    "pr = @time primes_blockmore(10^5, Int, 256);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(total_known[], next_check[]) = (256, 256)\n",
      " 64.109311 seconds (226.03 M allocations: 4.581 GiB, 7.13% gc time)\n"
     ]
    }
   ],
   "source": [
    "gc()\n",
    "@time primes_blockmore(5*10^5, Int, 256);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 39.164720 seconds (348.88 M allocations: 5.203 GiB, 1.81% gc time)\n"
     ]
    }
   ],
   "source": [
    "gc()\n",
    "@time primes_array(5*10^5, Int);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "So with all that, \n",
    "We are still losing to the single threaded code.\n",
    "Maybe if we were using more threads.\n",
    "Or if the code was smarter,\n",
    "we could pull ahead and go faster.\n",
    "But today, I am willing to admit defeat.\n",
    "It is really hard to make this kinda code actually speed-up.\n",
    "\n",
    "If you can do better, I'ld be keen to know.\n",
    "You can get the notebook that is behind this post from\n",
    "[github](https://github.com/oxinabox/oxinabox.github.io/blob/master/notebook_posts/Thread%20Parallelism%20in%20Julia.ipynb),\n",
    "you could even fork it and make a PR and I'll regenerate this blog post (:-D).\n",
    "\n",
    "One way to make it much much faster is to use a different algorithm.\n",
    "I'm sure there actually exist well documented parallel prime finders."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}