In the previous post we left off at fork join, shared versus distributed memory, and a first “hello world” with OpenMP. Now it’s time for the part you’re actually going to use every day in this parallel programming world, splitting a loop across threads, and not shooting yourself in the foot trying, because trust me it happens 😭
The for directive 🔁
#pragma omp for [clauses]
for-loopThis tells OpenMP to split the iterations of the following loop among the threads that already exist. Watch out for that detail, for doesn’t create new threads, it has to be inside a parallel region for there to be something to split work among. There’s no implicit barrier on entry, but there is one on exit.
#pragma omp parallel
{
#pragma omp for
for (i = 0; i < 1000; i++)
A[i] = B[i];
}With two threads, that loop of 1000 iterations automatically gets split into something like i=0..500 for Thread 1 and i=500..1000 for Thread 2. The variables A and B are shared, but i is local to each thread (remember, the index of an associated for is always private), and it takes a different initial and final value depending on the thread.
for versus working it out yourself 🧮
These two pieces of code are equivalent.
/* using the for directive */
#pragma omp parallel
#pragma omp for
for (i = 0; i < n; i++)
z[i] = a*x[i] + y;/* splitting by hand, with just parallel */
#pragma omp parallel private(id, num, istart, iend)
{
id = omp_get_thread_num();
num = omp_get_num_threads();
istart = id * n / num;
iend = min(n, (id + 1) * n / num);
for (i = istart; i < iend; i++)
z[i] = a*x[i] + y;
}The difference is pure manual work without for, we’re the ones who are going to calculate which range of i belongs to each thread using its id and the total number of threads. With for, OpenMP does that calculation for you. In practice you’ll almost always prefer for, unless you need very fine grained control over the split that the schedule clause doesn’t give you (more on that further down).
Exercise, before moving on. What’s the difference between these two programs?
/* program A */
#pragma omp parallel private(i)
{
for (i = 0; i < 10; i++)
printf("Hello world %d\n", i);
}/* program B */
#pragma omp parallel private(i)
{
#pragma omp for
for (i = 0; i < 10; i++)
printf("Hello world %d\n", i);
}Before scrolling down to look at the answer, please think about it for a bit, don’t read what I’ll write below, or well, sometimes a hint is fair too, but it really would be good if you tried to solve it yourself.
The answer. In program A, each thread runs the entire loop, all 10 iterations, because there’s no work sharing directive inside the parallel region. With 4 threads, you’ll see the messages “Hello world 0” through “Hello world 9” printed 4 times, one full batch per thread. In program B, #pragma omp for splits those 10 iterations among the available threads, so each line gets printed exactly once overall, split across threads. It’s the difference between “each thread does all the work” and “the work gets divided among threads”, and it’s exactly the kind of mistake you make when you forget the for inside a parallel.
Not every loop can be parallelized ⚠️
For for to work, the loop has to meet three restrictions.
- Be a structured block. No
breakorgototo exit early. - The number of iterations has to be computable ahead of time, before entering the loop.
- There can’t be dependencies between iterations.
The first two are mechanical, but the third one is the one that really matters to understand. The rule for spotting dependencies is simple to state. If you run the loop in reverse order and the result changes, the loop isn’t parallel. If it comes out the same regardless of order, then it is.
flowchart LR
subgraph OK["✅ A[i] = A[i] + B[i] → parallel"]
direction LR
a0[A0] -.-> a0
a1[A1] -.-> a1
a2[A2] -.-> a2
end
classDef ok fill:#dcefdd,stroke:#1f5c2e,color:#1f5c2e,font-weight:bold;
class a0,a1,a2 ok
style OK fill:#f2faf3,stroke:#4a8f57flowchart LR
subgraph FWD["🚫 A[i] = A[i+1] + B[i] → not parallel, forward dependency"]
direction RL
b1[A1] --> b0[A0]
b2[A2] --> b1
b3[A3] --> b2
end
classDef bad fill:#fbd8d3,stroke:#8a241b,color:#8a241b,font-weight:bold;
class b0,b1,b2,b3 bad
style FWD fill:#fdf3f2,stroke:#c0392bflowchart LR
subgraph BACK["🚫 A[i] = A[i-1] + B[i] → not parallel, backward dependency"]
direction LR
c0[A0] --> c1[A1]
c1 --> c2[A2]
c2 --> c3[A3]
end
classDef bad fill:#fbd8d3,stroke:#8a241b,color:#8a241b,font-weight:bold;
class c0,c1,c2,c3 bad
style BACK fill:#fdf3f2,stroke:#c0392bWith A[i] = A[i] + B[i], each iteration reads and writes only its own element, it doesn’t care at all what the others do, so you can run them in any order. With A[i] = A[i+1] + B[i], iteration i needs the original value of A[i+1], which iteration i+1 is about to overwrite. If i+1 runs first, i reads the wrong value. And with A[i] = A[i-1] + B[i], each iteration needs the result that the previous iteration just wrote, so they flatly have to run in order, one after the other. Neither of these last two cases can be split with for as it stands.
Second exercise. I know I’m starting to sound like a teacher, with quite a few exercises per lesson, but we really do need to exercise the mind a bit with this topic to be able to understand it since it can seem a bit like magic sometimes, now yes, the exercise, work out whether these three snippets are parallel, looking at which elements of A get read and written on each iteration.
/* snippet 1 */
for (i = 0; i < n; i++)
A[i] = A[i+n] + B[i];/* snippet 2 */
for (i = 0; i < n; i++)
A[i] = A[i+m] + B[i];/* snippet 3 */
for (i = 0; i < n; i++) {
C[i] = A[3*i+1] + 1;
A[2*i+7] = B[i] - 3;
}Snippet 1 is indeed parallel. Since i goes from 0 to n-1, the index i+n always falls outside the range the loop itself writes (A[i+n] never matches any A[i] that another iteration is modifying), so there’s no real conflict. Snippet 2 depends on m. If m is greater than or equal to n, it’s the same case as snippet 1 and it’s parallel. If m is smaller, A[i+m] can overlap with some A[j] that another iteration writes, and at that point you can’t guarantee anything without knowing the exact value. Snippet 3 is the most interesting one. It writes to A[2*i+7] and reads from A[3*i+1], two indices with different progressions. For there to be a real dependency, there would have to exist some pair of iterations i, j where 2*i+7 == 3*j+1, in general, for arbitrary values of n, that coincidence can indeed happen, so the snippet isn’t safe to parallelize without analyzing the concrete range of i. This last one is exactly the kind of case where “it looks parallel at a glance” but you have to do the math before trusting it.
Going deeper, splitting work without stepping on each other 🕵️
We now know how to identify what can be parallelized. Now let’s look at the tools for doing it right when the split isn’t as straightforward as adding two vectors.
Let’s look at the following practical example, you might have to read it a few times and work through it on your own, I didn’t get it on the first try, if I’m honest.
reduction, adding 100,000 numbers in one line ➕
#include <stdio.h>
#include <omp.h>
#define N 100000
int main(void) {
double A[N], sum = 0.0;
for (int i = 0; i < N; i++) A[i] = 1.0;
#pragma omp parallel for reduction(+:sum)
for (int i = 0; i < N; i++)
sum += A[i];
printf("Total sum = %.0f\n", sum); /* 100000 */
return 0;
}Under the hood, reduction(+:sum) does three things. It gives each thread its own private copy of sum, initialized to the neutral element of the operation (0 for +, 1 for *). Each thread accumulates into its copy only the elements it was assigned, without stepping on anyone else’s copy (which is why you don’t need critical or atomic inside the loop). And at the end of the region it combines all the copies into the original variable, with an internal tree reduction scheme, at each step pairs of values get added together and the number of values still alive gets cut roughly in half, so with P threads it takes on the order of log₂(P) steps instead of P-1 additions one by one.
flowchart TD
S0[sum₀] --> R1
S1[sum₁] --> R1[+]
S2[sum₂] --> R2
S3[sum₃] --> R2[+]
R1 --> R3[+]
R2 --> R3
R3 --> T[total]
classDef leaf fill:#dcefdd,stroke:#1f5c2e,color:#1f5c2e,font-weight:bold;
classDef op fill:#4a8f57,stroke:#1f5c2e,color:#ffffff,font-weight:bold;
classDef final fill:#c0392b,stroke:#8a241b,color:#ffffff,font-weight:bold;
class S0,S1,S2,S3 leaf
class R1,R2,R3 op
class T finalreduction also supports other operators, * (product), max/min, &&/||, &/|/^, each with its corresponding neutral element. One nuance worth keeping in mind. Floating point addition isn’t strictly associative because of rounding, so adding the same numbers in a different order (as reduction does when splitting them across threads, compared to adding them one by one sequentially) can produce a tiny difference in the last decimal places. It’s not a bug, it’s the price of parallelizing something that’s associative in pure math but not quite in real arithmetic.
Adding 100,000 numbers, step by step 🌳
reduction is convenient, but it feels almost magical if you’ve never seen what it does under the hood. So let’s unfold the same problem by hand, element by element, until it’s crystal clear why the code above has exactly that shape.
The problem. 10 threads add up those same 100,000 numbers stored in vector A. Each thread Pn (Pn = 0…9) handles its portion of 10,000 consecutive elements and stores the partial result in sum[Pn]. Up to this point there’s nothing to coordinate, each one adds its chunk completely independently.
sum[Pn] = 0;
for (i = 10000*Pn; i < 10000*(Pn+1); i++)
sum[Pn] = sum[Pn] + A[i]; /* sum of the local portion, no coordination needed */The interesting part starts now. Those 10 partial results have to be combined into a single total. The naive way would be for thread 0 to add sum[1]+sum[2]+…+sum[9] one at a time, 9 sequential additions, without taking advantage of having 10 free threads to help. The smarter trick is a tree reduction. At each step, half of the values still alive get added to their pair, and the number of values still alive is cut roughly in half. With P values it takes on the order of log₂(P) steps instead of P-1 additions one after another.
With 10 threads the difference (3 steps versus 9) isn’t too impressive, but it really shows when you scale up. With P = 1,024 threads, the naive sum would need 1,023 steps, one after another, while the tree does it in 10 steps (log₂(1024) = 10), because at each step all the threads still alive work at the same time instead of waiting their turn. That difference, linear versus logarithmic, is the underlying reason no parallel library adds up partial results one by one.
flowchart TB
subgraph I1["Start, 10 partial sums"]
direction LR
s0[S0] & s1[S1] & s2[S2] & s3[S3] & s4[S4] & s5[S5] & s6[S6] & s7[S7] & s8[S8] & s9[S9]
end
subgraph I2["Iteration 1, half, 10 → 5"]
direction LR
t0[S0+S5] & t1[S1+S6] & t2[S2+S7] & t3[S3+S8] & t4[S4+S9]
end
subgraph I3["Iteration 2, half odd, S4 gets absorbed, then 5 → 2"]
direction LR
u0["S0+S4+S2"] & u1["S1+S3"]
end
subgraph I4["Iteration 3, half, 2 → 1"]
v0["total in sum[0]"]
end
I1 --> I2 --> I3 --> I4
classDef alive fill:#dcefdd,stroke:#1f5c2e,color:#1f5c2e,font-weight:bold;
classDef dead fill:#f2f2f2,stroke:#9a9a9a,color:#9a9a9a;
classDef total fill:#c0392b,stroke:#8a241b,color:#ffffff,font-weight:bold;
class s0,s1,s2,s3,s4,t0,t1,t2,t3,t4,u0,u1 alive
class s5,s6,s7,s8,s9 dead
class v0 totalThe complete code, with the variable half keeping track of how many values are still alive at each moment.
half = 10; /* number of threads */
repeat
synch(); /* barrier, waits for the partial sum to finish */
if (half % 2 != 0 && Pn == 0)
sum[0] = sum[0] + sum[half-1];
half = half / 2;
if (Pn < half) sum[Pn] = sum[Pn] + sum[Pn+half];
until (half == 1); /* the final result ends up in sum[0] */synch() is a barrier. Nobody moves on to the next iteration until every thread has reached that point. And it’s not just for show. Without it, nothing stops a fast thread from entering iteration 2 while a slower one still hasn’t finished writing its result from iteration 1. A concrete interleaving that breaks the result, without synch(), in the first iteration (half going from 10 to 5, sum0 += sum5).
| Moment | P0 (computes sum0 += sum5) | P5 (still finishing its local sum) |
|---|---|---|
| t1 | reads sum5 → old value, not yet updated |
still in the loop adding up its portion of A |
| t2 | writes sum0 = sum0 + (old value of sum5) |
finishes adding and only then writes the correct sum5 |
The result. sum0 ends up with a total that doesn’t include all the elements P5 was supposed to add. A silent error, with no message and no program crash, and on top of that not reliably reproducible (sometimes P5 does make it in time and the result comes out right, depending on how loaded the machine is). That’s exactly what the barrier prevents, it forces all 10 threads to reach the same point before any of them reads another one’s result.
The complete trace, iteration by iteration.
| Iteration | half before → after | Operation after the barrier |
|---|---|---|
| 1 | 10 → 5 | sum0+=sum5; sum1+=sum6; sum2+=sum7; sum3+=sum8; sum4+=sum9 |
| 2 | 5 → 2 | half is odd, first sum0+=sum4 (correction), then sum0+=sum2; sum1+=sum3 |
| 3 | 2 → 1 | sum0+=sum1 → done, total in sum[0] |
And there it is, that, exactly that, is what the reduction(+:sum) clause saves you from writing by hand.
The most typical race condition 🏁
By default, a variable declared outside the parallel region is shared. Writing to it from several threads without protection is a race condition. The result depends on the order in which the executions interleave, and it can change from one run to the next.
int maximum = a[0];
#pragma omp parallel for
for (int i = 1; i < N; i++) {
if (a[i] > maximum) /* wrong! read+write without protection */
maximum = a[i]; /* on a shared variable */
}A concrete interleaving that produces a wrong result, with maximum at 5 and two threads processing a[i]=8 and a[j]=6 at almost the same time.
| Moment | Thread A (processes a[i]=8) | Thread B (processes a[j]=6) |
|---|---|---|
| t1 | reads maximum → 5 |
reads maximum → 5 |
| t2 | compares 8 > 5 → true | compares 6 > 5 → true |
| t3 | writes maximum = 6 |
|
| t4 | writes maximum = 8 |
With this order the result comes out right by chance (maximum = 8), but if the writes from t3 and t4 get swapped, Thread A writes 8 first and Thread B overwrites it with 6 afterward. The value 8 gets lost, even though both threads “saw” at some point that 8 was greater. This is called a lost update. The problem isn’t that a thread reads garbage, it’s that one thread’s read compare write sequence gets interleaved with another’s without either one noticing the other’s progress.
The correct fix isn’t to declare maximum as private (it would lose its value on exiting the loop), but to use reduction here too, with the max operator.
maximum = a[0];
#pragma omp parallel for reduction(max:maximum)
for (int i = 1; i < N; i++)
if (a[i] > maximum) maximum = a[i]; /* correct now */private(x) is still necessary for auxiliary variables that each iteration recalculates from scratch (a temporary index, a local accumulator that isn’t read from one iteration to the next) and that precisely because of that shouldn’t be shared between threads.
critical and atomic, protecting a shared update 🔒
atomic protects a single memory operation (reading, modifying, and writing a variable) and usually translates to an atomic hardware instruction. It’s the fastest option when it applies. critical protects an arbitrary block of code, with more overhead but more flexibility.
long counter = 0;
#pragma omp parallel for
for (int i = 0; i < 1000000; i++) {
#pragma omp atomic
counter++; /* a single operation → atomic is enough */
}
printf("counter = %ld (expected: 1000000)\n", counter);#pragma omp parallel for
for (int i = 0; i < n; i++) {
double local = f(i); /* expensive computation, outside the critical section */
#pragma omp critical
{ /* if + two writes, doesn't fit in an atomic */
if (local > best_value) {
best_value = local;
best_index = i;
}
}
}atomic only supports very specific ways of updating a single scalar variable, x++, x--, x op= expr, or x = x op expr. That’s why the second example can’t use atomic. There’s an if and two different variables at once, it needs critical. And a performance warning worth remembering. Putting a critical inside a heavily iterated loop serializes exactly that part. If most of the loop’s time is spent in the critical section, parallelism doesn’t help you much, no matter how fast the rest of the body is.
With this we now have the essentials of OpenMP, how to split a loop, how to know if it can be split, and how to protect yourself from stepping on your own feet when several threads touch the same data.
The song of the post
Ghost Town by Kanye West. Honestly it’s nothing special, but I’ve been listening to it a ton lately and it’s a pretty good song, in fact I listened to it several times while writing this post haha.