JM

High Performance Computing in Julia

Course Info

  1. 1. Introduction

    1. 1.1. Setup

  2. 2. Foundations

    1. 2.1. Overview
    2. 2.2. von Neumann Architecture
    3. 2.3. Data Types
    4. 2.4. Operating Systems
    5. 2.5. Compilers and Machine Code
    6. 2.6. Hardware
    7. 2.7. Summary

  3. 3. Julia

    1. 3.1. Overview

  4. 4. Measuring Performance

    1. 4.1. Overview
    2. 4.2. Timing
    3. 4.3. Benchmarking
    4. 4.4. Profiling
    5. 4.5. Memory Allocations
    6. 4.6. Computational Complexity

  5. 5. Optimisation

    1. 5.1. Overview
    2. 5.2. Caching and Memory Locality
    3. 5.3. Stack vs Heap Memory
    4. 5.4. Reducing Allocations
    5. 5.5. Memoisation
    6. 5.6. Vectorisation vs Loops
    7. 5.7. Views
    8. 5.8. Speed of Operations
    9. 5.9. Type Stability
    10. 5.10. Performance Annotations
    11. 5.11. Conclusion

  6. 6. Parallel Programming

    1. 6.1. Overview
    2. 6.2. Dependency Graphs
    3. 6.3. Theoretical Expectations
    4. 6.4. Maps and Reductions
    5. 6.5. Thread Safety and Race Conditions
    6. 6.6. Producer and Consumer
    7. 6.7. When to Parallelise?

  7. 7. Multithreading

    1. 7.1. Overview
    2. 7.2. Monte Carlo Pi Estimation
    3. 7.3. Naive Parallel Implementation
    4. 7.4. False Sharing
    5. 7.5. Chunking
    6. 7.6. Task-based Parallelism
    7. 7.7. Packages
    8. 7.8. Summary
    9. 7.9. Exercises

  8. 8. Multiprocessing

    1. 8.1. Overview
    2. 8.2. Multiprocessing in Julia
    3. 8.3. Multiprocessing on a Cluster
    4. 8.4. Message Passing Interface
    5. 8.5. Summary
    6. 8.6. Exercises

  9. 9. GPU Programming

    1. 9.1. Introduction to GPU Programming
    2. 9.2. Modern GPU Hardware
    3. 9.3. GPU Vendor Overview
    4. 9.4. CUDA.jl
    5. 9.5. AMDGPU.jl
    6. 9.6. Portable GPU Kernels
    7. 9.7. CUDA Libraries
    8. 9.8. Benchmarking & Profiling
    9. 9.9. Tips
    10. 9.10. Case Study - Monte-Carlo Simulations
    11. 9.11. Exercises

Table of Contents

Naive Parallel Implementation

Each dart thrown requires two randomly generated numbers, mapped to be between -1 and +1. Additionally, they must calculate the square of the distance from the dart to the origin and check to see whether it is less than 1. All of these steps can be done in parallel, however, the variable n_c is shared between the threads and so provides a race condition if multiple workers were to try and alter it at the same time. In order to solve this, we can have a variable for each thread and then separately sum these at the end.

julia 
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

Here, we have used very little memory to keep track of the counts, and made sure the implementation was thread-safe as each thread accessed a different part of memory. It should be noted here that Threads.threadid() returns an ID for the currently executing thread within the threaded loop. Additionally, the rand() function call is thread-safe on recent versions of Julia (v1.6+). Note that this implementation, even when seeded, will not be deterministic. Even if each thread had its own seeded RNG, then additional randomness is introduced when scheduling the threads for execution. This is usually not an issue, but it is something that should be kept in mind, especially when one would like to write testable code.

Now, we can benchmark these algorithms and see which is fastest. We can calculate the efficiency based on the number of threads:

julia 
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

num_threads = Threads.nthreads()

The benchmarks:

julia 
n = 100_000_000;
mc_pi_serial_time = @belapsed est_pi_mc_serial($n)
mc_pi_threaded_time = @belapsed est_pi_mc_threaded($n)
@show mc_pi_serial_time
@show mc_pi_threaded_time
speedup = mc_pi_serial_time / mc_pi_threaded_time
efficiency = speedup / num_threads
@show speedup
@show efficiency
Output 
mc_pi_serial_time = 0.232457664
mc_pi_threaded_time = 0.239279096
speedup = 0.9714917344890003
efficiency = 0.12143646681112504

One would expect the ratio of those times to approach the number of available threads, however, even for this very large value of n, the parallel efficiency is tiny. In fact, on the machine these notes were built on, the threaded version is actually slower than the serial one. This is very likely to be a case of false sharing.

Back Next