Parallel programming, introduction #1

Contents

Why does HPC exist? 🚀
Two HPC architectures 🏗️
The four phases of designing something parallel 🪜
Shared memory in practice, OpenMP 🧵
What a directive looks like 🔍
The OpenMP memory model 💾
The parallel directive, in detail 🔀
Compiling and running on Ubuntu 🐧
Hello world, with threads 👋
The song of the post

Hellooooo, it’s me again, but with a new topic “Parallel programming”, I’m studying it at the time of writing this post for my master’s degree. So there’s no better way to learn than writing yourself about what you learned

Why does HPC exist? 🚀

HPC (High Performance Computing) is the fancy way of saying “supercomputers and clusters in research centers, companies, and universities”. But the underlying goal isn’t that different from programming on your laptop, you want the same result, just using many processors or cores at once, to get there faster or to be able to solve something a single processor couldn’t even attempt (because of time, or simply because of lack of memory).

When do you really need it? Basically in three cases, when what you need to compute is more than a normal desktop can handle, when the problem demands more precise mathematical models (and therefore more expensive in compute, memory, and disk), or when the result is only useful if it arrives on time (a 24 hour weather forecast that takes 3 days to calculate is useless to anyone).

There’s even a group of problems, the Grand Challenge applications, that can only be tackled on HPC systems, large scale molecular dynamics, global climate simulation, that kind of thing. The two usual obstacles to getting into that world are access to those machines and the difficulty of programming them. Parallelizing forces you to think differently, it’s not just “the same code but faster”.

Most real HPC applications are modeling and simulation, automotive, naval and aerospace industry, structural engineering, combustion simulation, fluid mechanics, computational chemistry, bioinformatics, weather forecasting, seismology, circuit design. And underneath all of that there’s a handful of algorithms that keep showing up over and over, so it’s worth recognizing them, Monte Carlo simulations, N-body problems, the FFT, graph partitioning, genetic algorithms, dense and sparse matrix algebra.

A random but interesting fact, there’s the TOP500, a ranking of the most powerful supercomputers on the planet. What matters there isn’t so much who’s first, but the trend, industry and web services carry more and more weight compared to the purely academic use of a decade ago. HPC stopped being just a thing for physics labs.

Now, relax, the vast majority of the problems you’re going to solve in your life don’t need anything from the TOP500, unless you’re a researcher, mind you, nobody’s limiting anybody here.

But in reality a small multicore cluster, with a few nodes each having several cores, is the HPC architecture you’re actually going to run into day to day.

Two HPC architectures 🏗️

Well, to simplify a lot of things, we can lump the whole HPC world into 2 architectures based on how the processors have to communicate with each other, do they all see the same memory, or does each one have its own?

flowchart LR
    subgraph SC["🟢 Shared memory"]
        direction TB
        P0[P0] --- M1[(shared memory)]
        P1[P1] --- M1
        P2[P2] --- M1
        P3[P3] --- M1
    end
    subgraph SD["🔴 Distributed memory"]
        direction TB
        Q0[P0] --- L0[(local mem.)]
        Q1[P1] --- L1[(local mem.)]
        Q2[P2] --- L2[(local mem.)]
        Q3[P3] --- L3[(local mem.)]
        Q0 <--> RED[interconnect network]
        Q1 <--> RED
        Q2 <--> RED
        Q3 <--> RED
    end

    classDef proc fill:#dcefdd,stroke:#1f5c2e,color:#1f5c2e,font-weight:bold;
    classDef mem fill:#4a8f57,stroke:#1f5c2e,color:#ffffff,font-weight:bold;
    classDef dproc fill:#fbd8d3,stroke:#8a241b,color:#8a241b,font-weight:bold;
    classDef dmem fill:#c0392b,stroke:#8a241b,color:#ffffff,font-weight:bold;
    class P0,P1,P2,P3 proc
    class M1 mem
    class Q0,Q1,Q2,Q3 dproc
    class RED,L0,L1,L2,L3 dmem
    style SC fill:#f2faf3,stroke:#4a8f57
    style SD fill:#fdf3f2,stroke:#c0392b

In shared memory (multiprocessors, nowadays basically any multicore CPU), all processors see a single address space. You don’t have to distribute anything, you communicate implicitly, reading and writing shared variables. The downside is that you have to synchronize by hand (semaphores, critical sections, barriers) so concurrent accesses don’t step on each other.

In distributed memory (multicomputers, nowadays clusters), each processor has its own private address space. You have to split up the data yourself across the local memories, and know at all times where everything is (that’s called exploiting locality). Communication is explicit, through messages, which is more work to program, but synchronization comes almost for free because it’s implicit in the message itself.

The mental rule that stuck with me while studying this topic is very basic, shared is easy communication but manual synchronization, distributed is almost free synchronization but manual communication.

Okay now we can say in a very rough way that all the other HPC architectures (hybrid, PGAS, data parallel) are variations or mixes of these two ideas.

The four phases of designing something parallel 🪜

For me this is something new. Let’s see, I come from programming concurrency in Go but not large distributed computing systems with shared memory and so on, so what I do day to day is really pretty sequential programming, and I didn’t know there were design rules for this.

Well, anyway, no matter which paradigm you end up using, designing a parallel program always goes through these four decisions, in this order.

flowchart TD
    A["🔪 1. Decomposition
    split the computation into tasks"] --> B["📦 2. Assignment
    static or dynamic"]
    B --> C["🔗 3. Coordination
    communication and synchronization"]
    C --> D["🗺️ 4. Mapping
    which process runs on which core"]

    classDef fase fill:#dcefdd,stroke:#1f5c2e,color:#1f5c2e,font-weight:bold;
    class A,B,C,D fase

Decomposition. Split the sequential computation into tasks that can be handed out. The goal is load balancing, making sure no processor sits waiting while another keeps working.

Assignment. How you hand out those tasks among processes or threads. It can be static (decided before running, simpler) or dynamic (decided at runtime, better when the load is irregular).

Coordination. The communications and synchronizations needed between tasks. Golden rule, minimize them, because they’re pure overhead compared to the sequential version.

Mapping. Decide which process or thread runs on which physical processor. In practice it’s almost always 1:1, one process or thread per core.

A piece of advice that comes up a lot in the parallel programming world is that you shouldn’t try to optimize the whole program equally. Focus on the hot spots, the most expensive areas (typically loops that iterate over mountains of data), because that’s where parallelism really justifies the effort of programming it.

Shared memory in practice, OpenMP 🧵

The typical programming model for shared memory is thread based concurrency, following the fork join pattern, a master thread runs sequentially until it reaches a parallel region, there it “launches” (fork) a group of threads that work in parallel, and when they all finish they rejoin (join) into the master, which keeps going alone until the next parallel region.

flowchart LR
    M1((master)) --> F1{fork}
    F1 --> T1[thread]
    F1 --> T2[thread]
    F1 --> T3[thread]
    T1 --> J1{join}
    T2 --> J1
    T3 --> J1
    J1 --> M2((master)) --> F2{fork}
    F2 --> T4[thread]
    F2 --> T5[thread]
    T4 --> J2{join}
    T5 --> J2
    J2 --> M3((master))

    classDef master fill:#141414,stroke:#141414,color:#ffffff,font-weight:bold;
    classDef fj fill:#c0392b,stroke:#8a241b,color:#ffffff,font-weight:bold;
    classDef thread fill:#dcefdd,stroke:#1f5c2e,color:#1f5c2e,font-weight:bold;
    class M1,M2,M3 master
    class F1,F2,J1,J2 fj
    class T1,T2,T3,T4,T5 thread

You can program this with native threads (Pthreads, for example in C), portable but low level, and you have to create, synchronize, and destroy threads by hand with library calls and blah blah blah, things that really complicate our lives. The good thing about all this is that there’s a high level alternative that’s practically an industry standard, OpenMP. It’s a de facto standard for programming shared memory systems in Fortran, C, and C++, made of directives that get inserted on top of your sequential code without changing it almost at all. Beautiful, honestly 🤯

OpenMP relies on three elements.

Parallelism control. The parallel directive, and work sharing directives like for.
Data and communication control. The shared and private clauses.
Synchronization. Barriers, critical sections, atomic.

Somewhat important information, to keep in mind what OpenMP actually is. It’s maintained by the OpenMP Architecture Review Board (hardware, software, and computing center vendors, AMD, ARM, Intel, IBM, among others), the specs are free, you don’t need a license, and all the info is at openmp.org.

What a directive looks like 🔍

In C and C++, OpenMP directives use the language’s own #pragma mechanism.

#pragma omp directive [clauses]

Always in lowercase. And here comes a detail. If the compiler doesn’t recognize the pragma, it simply ignores it. This happens for three reasons, the compiler doesn’t support OpenMP, you forgot the flag to enable it (-fopenmp in gcc), or there’s a typo (for example you wrote #pragma openmp instead of #pragma omp). In all three cases your program compiles with no error, but runs completely sequentially. If your “parallel program” always behaves as if it had a single thread, this is the first suspect.

Most directives apply to a structured block. Code with no jumps that enter or leave the block from outside (no goto crossing the boundary), WATCH OUT FOR THIS!!!

The OpenMP memory model 💾

OpenMP lets you define two types of variables.

shared. A single copy in memory, visible to all threads.
private. Each thread has its own copy, invisible to the others.

Each thread also has its own stack (to store the arguments and local variables of the functions that thread calls), and that’s where the space for private variables lives. The size of that stack depends on the implementation, but it can be configured with the OMP_STACKSIZE environment variable.

Updating (writing) a shared variable can get expensive if there are many threads involved, or if you’re on a NUMA architecture and the access is to remote memory on another processor instead of local memory. Private variables reduce how often you have to touch shared variables (less overhead), at the cost of using more total memory for the program. In general it’s worth using shared variables when they’re read only, when different threads access different elements of the same variable, or precisely when you want to communicate a value between threads.

An important nuance here, updates to shared variables can temporarily stay in each thread’s private view (cache), so at any given moment two threads could see different values of the same shared variable. OpenMP requires shared objects to be made visible (actually written to memory) at synchronization points. If a thread needs a value that another thread wrote, there has to be a synchronization point in between. To force this manually there’s the flush directive, and if you’ve studied C you probably already know it.

The `parallel` directive, in detail 🔀

#pragma omp parallel [clauses]
{
    /* structured block */
}

It defines a parallel region, a block of code that several threads run at once. When a thread reaches this directive, it creates a team of threads, becomes the master of that team, and the number of threads in the team is decided by a clause, the OMP_NUM_THREADS environment variable, or a call to omp_set_num_threads(). At the end of the region there’s an implicit barrier, and only the master thread keeps going.

double A[1000];
int ID;
omp_set_num_threads(4);

#pragma omp parallel private(ID) shared(A)
{
    ID = omp_get_thread_num();
    compute(ID, A);
}

printf("\ndone");

With this, 4 threads are created inside the parallel region, and each one runs compute with its own ID (0, 1, 2, or 3), all sharing the same copy of A.

The most important clauses of parallel.

if(expression). The region only runs in parallel if the expression is true.
num_threads(n). Sets the number of threads in the team.
private(list). Each thread gets its own local, uninitialized copy of those variables.
firstprivate(list). Same as private, but the copy starts with the value it had before entering the region.
shared(list). A single copy, visible to everyone, with no guaranteed mutual exclusion (that’s your responsibility).
default(shared | none). Defines the default behavior. With none it forces you to explicitly specify the type of every variable, which is generally a good idea because it forces you to think about each one.
reduction(operator:list). You’ll see it in detail in the next post.
copyin(list). For threadprivate variables, you won’t need it often.

By default, global, static, and dynamic duration data are shared, local variables are private, and the index of a loop associated with a for (or parallel for) is always private. Even so, it’s strongly recommended to explicitly specify the type of every variable instead of relying on the default value.

If a thread from a team that’s already running a parallel region finds another parallel region inside it, it creates a new team and becomes its master (this is called nested parallelism). By default, nested regions come serialized, the new team has a single thread, but this can be changed with OMP_NESTED or omp_set_nested().

Compiling and running on Ubuntu 🐧

Before touching any code, on Ubuntu (or any Debian-like) the standard compiler is GCC, and it’s come with built-in OpenMP support for many versions now. You don’t need to install any separate library, just add a flag.

First, check which compiler you have.

gcc --version
g++ --version

If you get something like gcc (Ubuntu 13.2.0-...) 13.2.0, you already have a C and C++ compiler!!! Great. If it says command not found, install the package that brings both at once.

sudo apt update
sudo apt install build-essential

build-essential installs gcc, g++, make, and the standard C/C++ headers for you.

Next, confirm that compiler really supports OpenMP.

echo | gcc -fopenmp -dM -E - | grep _OPENMP

It should print something like #define _OPENMP 201511 (that number is the date, year and month, of the standard version it supports, 201511 is OpenMP 4.5). If it prints nothing, that gcc was built without OpenMP support, which is pretty rare on Ubuntu.

Check how many cores your machine has, because in practice that number is the limit of threads it makes sense to test your programs with.

nproc

And these are the four lines you’re going to reuse for every example in this series.

gcc  -fopenmp -O2 file.c   -o program      # compile C
g++  -fopenmp -O2 file.cpp -o program      # compile C++
OMP_NUM_THREADS=4 ./program                # run forcing 4 threads
./program                                  # run using all available cores

The most common mistake when starting out. Forgetting -fopenmp. The program compiles with no error at all (the compiler ignores the #pragma omp it doesn’t recognize), but runs completely sequentially. If your “parallel” program always prints that there’s only one thread running, check that flag before anything else.

Hello world, with threads 👋

Lots of theory and little practice, let’s wrap up with the mandatory example in any language, a Hello world, but a parallel region that identifies each thread.

#include <stdio.h>
#include <omp.h>

int main(void) {
    #pragma omp parallel
    {
        int id       = omp_get_thread_num();
        int nthreads = omp_get_num_threads();
        printf("Thread %d of %d: Hello world!\n", id, nthreads);
    }
    return 0;
}
/* output (nondeterministic order):
   Thread 2 of 4: Hello world!
   Thread 0 of 4: Hello world!
   Thread 3 of 4: Hello world!
   Thread 1 of 4: Hello world!         */

This is exactly the fork join from above. Up to #pragma omp parallel there’s a single thread (the master) running main sequentially. There the additional threads are created (by default as many as the machine has cores, unless you set OMP_NUM_THREADS or num_threads(n)), and all of them, including the master, run a copy of the block. When the closing brace is reached, the additional threads are destroyed and you go back to having a single thread of control.

Two things that surprise people the first time, the output order isn’t guaranteed (the operating system decides in which order each thread runs, so each run can print in a different order), and in theory, if several threads call printf at almost the same time, the output could get interleaved character by character. In practice most libc implementations protect each printf call as a unit, but that’s a guarantee from the system’s input and output implementation, not from the OpenMP standard.

With this we now have more or less an idea, we know why HPC exists and how to program with OpenMP. There’s still a lot left, but it’s a good start.

The song of the post

WHAT A GREAT SONG TO START SOMETHING WITH, I WON’T SAY ANYTHING ELSE, IT’S REALLY GOOD.

Last Nite
· The Strokes

Listening to