Chunking

The problem of false sharing can be avoided by chunking the larger calculation into much smaller ones. An additional benefit of this method is that we can use the original implementation to maximise code reuse. If we make any performance improvements to the serial algorithm, this will also improve the parallel implementation.

A very helpful method for chunking the operations is from the base Iterators library called Iterators.partition, which breaks a collection up into smaller blocks. Be careful with its second argument: Iterators.partition(collection, chunk_size) takes the size of each block, not the number of blocks. It yields blocks of chunk_size elements each, so the number of blocks you get back is length(collection)/chunk_size, rounded up. To end up with one block per thread we therefore have to work out the chunk size ourselves first, which is what the div(n, n_threads, RoundUp) line below does.

function est_pi_mc_threaded_chunked(n)
    n_threads = Threads.nthreads()
    num_inside = zeros(Float64, n_threads)
    # Calculate maximum chunk size
    chunk_size = div(n, n_threads, RoundUp)

    # Create an iterator and collect to turn into an array
    iter = collect(enumerate(Iterators.partition(1:n, chunk_size)))
    Threads.@threads for info in iter
        i, idx_range = info # Unpack the tuple from enumerate
        n_block = length(idx_range)
        pi_est = est_pi_mc_serial(n_block)
        num_inside[i] = pi_est*n_block/4
    end

    n_c = sum(num_inside)
    return 4 * n_c / n
end

Now, we can benchmark this final algorithm and compare it to the others:

function est_pi_mc_serial(n)
    n_c = zero(typeof(n))
    for _ in 1:n
        # Choose random numbers between -1 and +1 for x and y
        x = rand() * 2 - 1
        y = rand() * 2 - 1
        # Work out the distance from origin using Pythagoras
        r2 = x*x+y*y
        # Count point if it is inside the circle (r^2=1)
        if r2 <= 1
            n_c += 1
        end
    end
    return 4 * n_c / n
end

function est_pi_mc_threaded(n)
    n_cs = zeros(typeof(n), Threads.maxthreadid())
    Threads.@threads for _ in 1:n
        # Choose random numbers between -1 and +1 for x and y
        x = rand() * 2 - 1
        y = rand() * 2 - 1
        # Work out the distance from origin using Pythagoras
        r2 = x*x+y*y
        # Count point if it is inside the circle (r^2=1)
        if r2 <= 1
            n_cs[Threads.threadid()] += 1
        end
    end
    n_c = sum(n_cs)
    return 4 * n_c / n
end

function est_pi_mc_threaded_chunked(n)
    n_threads = Threads.nthreads()
    num_inside = zeros(Float64, n_threads)
    # Calculate maximum chunk size
    chunk_size = div(n, n_threads, RoundUp)

    # Create an iterator and collect to turn into an array
    iter = collect(enumerate(Iterators.partition(1:n, chunk_size)))
    Threads.@threads for info in iter
        i, idx_range = info # Unpack the tuple from enumerate
        n_block = length(idx_range)
        pi_est = est_pi_mc_serial(n_block)
        num_inside[i] = pi_est*n_block/4
    end

    n_c = sum(num_inside)
    return 4 * n_c / n
end

n = 100_000_000
mc_pi_serial_time = @belapsed est_pi_mc_serial($n)
mc_pi_threaded_time = @belapsed est_pi_mc_threaded($n)
mc_pi_threaded_chunked_time = @belapsed est_pi_mc_threaded_chunked($n)
@show mc_pi_serial_time / mc_pi_threaded_time
@show mc_pi_serial_time / mc_pi_threaded_chunked_time
@show mc_pi_threaded_time / mc_pi_threaded_chunked_time

mc_pi_serial_time / mc_pi_threaded_time = 0.8486258897387036
mc_pi_serial_time / mc_pi_threaded_chunked_time = 6.209387009511506
mc_pi_threaded_time / mc_pi_threaded_chunked_time = 7.316989835678251

We see that the chunked approach was much faster, and almost reached the theoretical maximum performance. In general, it is better split parallel tasks into large chunks that can be sequentially processed by that chunk. This avoids a lot of scheduling and orchestration overhead when managing the threads, as most variables can live inside the stack with little need to coordinate execution. Additionally, we massively reduced the number of writes to memory with the chunked approach as the count could live in registers close to the CPU and only be saved to memory once the bulk of the calculation was completed.

However, it should be noted that the overhead of parallel execution can be significant. The only way to see the effect is to measure the performance relative to the input size $n$. We can already assess that this algorithm has a time complexity of $\mathcal{O}(n)$. Instead of plotting these lines together, we will plot the speed-up $S$, compared to the theoretical maximum given by Amdahl’s law, which for an embarrassingly parallel problem split into $w$ equal chunks is $S = N / \lceil N/w \rceil$, as derived in the Theoretical Expectations chapter.

Inspecting Figure 1, we can see that the chunked implementation approaches the maximum speed-up for this algorithm, but suffers at lower values of $n$.

One can infer that the cost of using multithreading is quite high, especially when the contents of the for loop are computationally inexpensive. However, this cost is more or less constant and increasing the throughput will minimise the relative size of this cost. It takes on the order of $10^4$ samples before the threaded version even breaks even with the serial one, and $10^6$ or more before the speed-up gets close to the theoretical maximum. Below that, the overhead of the threaded approach is greater than the time taken to perform the calculations.