OpenMP: The Parallel For Directive

Module 13

Let OpenMP divide the loop

In Module 12, we divided the trapezoids ourselves. Each thread asked for its number and the team size, then took every team_size-th trapezoid, just as our MPI processes did. That works, but the loop no longer looks like the serial program.

Most of the work in a scientific program happens in for loops. OpenMP’s parallel for directive starts a team of threads and divides a loop’s iterations among them. Often one line turns a serial loop into a parallel one. This module follows §5.5 of Pacheco and Malensek. The catch, and the main topic of today, is that the compiler does not check whether that line is safe.

Our three programs are:

  • omp_trap_for.c — The trapezoid rule again, with one directive on the serial loop.
  • omp_fibo.c — A loop that compiles and runs but should not be parallelized.
  • omp_pi.c — Estimate π with a series, fixing a dependence and a scope bug along the way.

As in Module 12, compile and run these on PicoCluster with -fopenmp.

First program: the trapezoid rule with one directive

Here is the core of omp_trap_for.c. Compare it with omp_trap_reduction.c from last time.

double total_area = 0.0;

double start = omp_get_wtime();
#   pragma omp parallel for num_threads(thread_count) reduction(+: total_area)
for (int i = 0; i < n; i++) {
    double left = a + i * h;
    double right = a + (i + 1) * h;
    total_area += h * (f(left) + f(right)) / 2.0;
}
double elapsed = omp_get_wtime() - start;

Delete the # pragma line and this is exactly the serial loop. Gone are my_rank, team_size, the my_area variable, and the i += team_size stride. The structured block after parallel for must be a for loop, and OpenMP decides which iterations each thread runs.

gcc -Wall -O2 -fopenmp -o omp_trap_for omp_trap_for.c
./omp_trap_for 1
./omp_trap_for 4

The estimate should match Module 12, and the times should be similar to omp_trap_reduction.

Which thread gets which iterations?

With 4 threads and 500,000,000 trapezoids, most systems give each thread one block of consecutive iterations: thread 0 takes roughly the first quarter, thread 1 the next quarter, and so on. Our Module 12 loop used a cyclic division. Either works here because every trapezoid costs about the same. We will see how to choose the division with a schedule clause later in Chapter 5.

What is shared and what is private?

The rules from Module 12 still apply, with one addition:

  • Variables declared before the loop (a, h, n) are shared.
  • Variables declared inside the loop body (left, right) are private.
  • The loop variable i is private, even if it were declared before the loop. If all threads shared one i, every i++ would be a race.
  • total_area is a reduction variable. Without the reduction clause it would be shared, and total_area += … would be an unprotected update. That is the bug from step 2 of last time’s experiment.

Which loops can OpenMP divide?

OpenMP has to know how many iterations there are before the loop starts, so it can hand them out. It only accepts for loops in canonical form (§5.5.1):

for (index = start; index < end; index++)

The test may use <, <=, >, or >=. The update may add or subtract a fixed increment. index must be an integer (or pointer), and start, end, and the increment must not change during the loop. The loop body must not change index either.

So while loops, for (;;) loops, and loops that can break early cannot use parallel for. Pacheco tries to parallelize a linear search:

#   pragma omp parallel for num_threads(thread_count)
for (i = 0; i < n; i++) {
    if (A[i] == key) return i;
}

gcc rejects it: invalid exit from OpenMP structured block. A return or break would leave the other threads with work assigned and nowhere to go.

These rules are checked by the compiler. The next problem is not.

Second program: a loop that should not be parallel

The Fibonacci numbers are \(1, 1, 2, 3, 5, 8, 13, \ldots\), where each is the sum of the two before it. The core of omp_fibo.c is Pacheco’s example:

fibo[0] = fibo[1] = 1;
#   pragma omp parallel for num_threads(thread_count)
for (int i = 2; i < N; i++) {
    fibo[i] = fibo[i - 1] + fibo[i - 2];
}

This loop is in canonical form, so it compiles without a warning. The program also fills a check array with an ordinary serial loop and reports the first value where the two arrays disagree.

Predict: With 2 threads and N = 50, thread 0 gets roughly i = 2 through 25 and thread 1 gets roughly 26 through 49. What does thread 1 need when it computes fibo[26]? Who computes it, and when?

gcc -Wall -O2 -fopenmp -o omp_fibo omp_fibo.c
./omp_fibo 1
./omp_fibo 2
./omp_fibo 4

Run the 2- and 4-thread versions many times. Some runs may be correct, because thread 0 can finish its short block before thread 1 starts. When a run is wrong, the wrong values are often zeros: fibo is a global array, so it starts filled with zeros, and a thread read entries no one had written yet. (A one-line shell loop helps: for t in 1 2 3 4 5 6 7 8; do ./omp_fibo 4; done.)

Loop-carried dependences

Iteration i reads fibo[i - 1] and fibo[i - 2], which are written by earlier iterations. When one iteration uses a result computed in a different iteration, the loop has a loop-carried dependence (§5.5.2). parallel for assumes the iterations are independent and can run in any order on any thread. A loop with a loop-carried dependence breaks that assumption.

Two points from the book:

  1. OpenMP does not check for loop-carried dependences. Finding them is our job.
  2. A loop like this one cannot, in general, be fixed by adding a clause. The calculation itself is sequential.

Not every dependence matters. In this loop

for (int i = 0; i < n; i++) {
    x[i] = a + i * h;
    y[i] = exp(x[i]);
}

the second line uses the result of the first, but in the same iteration, which always runs on one thread. That loop is safe to parallelize.

To look for trouble, check each variable the loop body writes: is it read or written in a different iteration? (§5.5.3)

Third program: estimating π

A classic (and slow) series for π is

\[ \pi = 4\left(1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots\right) = 4\sum_{k=0}^{\infty} \frac{(-1)^k}{2k+1}. \]

The serial code flips the sign each time through:

double factor = 1.0;
double sum = 0.0;
for (long long k = 0; k < n; k++) {
    sum += factor / (2 * k + 1);
    factor = -factor;
}
pi_approx = 4.0 * sum;

Discuss: We need reduction(+: sum). Is that enough? Look for a variable that one iteration writes and another iteration reads.

A loop-carried dependence. factor = -factor in iteration \(k\) sets the sign that iteration \(k + 1\) uses. If those iterations are on different threads, the signs are wrong. The fix is to compute the sign from \(k\) alone: it is \(+1\) when \(k\) is even and \(-1\) when \(k\) is odd.

factor = (k % 2 == 0) ? 1.0 : -1.0;
sum += factor / (2 * k + 1);

A scope bug. Even with that fix, factor was declared before the loop, so it is shared. Thread 0 can set it to 1.0, then thread 1 sets it to −1.0 before thread 0 uses it. Each thread needs its own factor.

The core of the corrected omp_pi.c is:

double factor;
double sum = 0.0;

#   pragma omp parallel for num_threads(thread_count) \
        default(none) reduction(+: sum) private(factor) shared(n)
for (long long k = 0; k < n; k++) {
    factor = (k % 2 == 0) ? 1.0 : -1.0;
    sum += factor / (2 * k + 1);
}
double pi_approx = 4.0 * sum;

The backslash continues the directive onto the next line. The new clauses are:

Clause Meaning
private(factor) Each thread gets its own factor.
shared(n) All threads use the one n.
default(none) Make the compiler require a scope for every variable from outside the loop that the loop uses.

With default(none), forgetting factor is a compile error instead of a silent race. Pacheco recommends it as a good habit (§5.5.5). The loop variable k and anything declared inside the loop are private automatically, so they need no clause.

A private copy starts uninitialized. It does not start with the value from before the loop, and its value is not copied back afterward. That is fine here, because every iteration assigns factor before using it.

Another fix: declare it inside the loop body as double factor = (k % 2 == 0) ? 1.0 : -1.0;. Then it is private automatically, and the private clause is not needed.

Build and check it

The program takes the thread count and, optionally, the number of terms (default 100,000,000):

gcc -Wall -O2 -fopenmp -o omp_pi omp_pi.c
./omp_pi 1 1000
./omp_pi 4 1000
./omp_pi 1
./omp_pi 4

With 1,000 terms, both runs should print 3.14059265383979, matching the one-thread result in the book. The series converges slowly: with \(n\) terms, the error is about \(1/n\). So 100,000,000 terms gives about eight correct decimal places. Record your times for 1, 2, and 4 threads, as we did for the trapezoid rule.

Experiment: break the π program on purpose

Start from the working program. For each change, predict the result first, then recompile and run with ./omp_pi 1 1000 and ./omp_pi 4 1000 a few times.

  1. Remove private(factor) but keep default(none). What does the compiler say?
  2. Also remove default(none). Now the program compiles. Is the 4-thread estimate right? Is it the same every run?
  3. Restore both clauses, then put back the serial sign flip: replace the factor = (k % 2 … line with factor = -factor; placed after the sum line. Why can’t this work, even though factor is private? (Hint: what is a private variable’s value at the start of each thread’s first iteration?)
ImportantWrong answers that look reasonable

In step 2, the book’s two-thread runs gave estimates such as 2.97 and 3.22. Those are close enough to look plausible if you do not compare them with a one-thread run. Always compare with a trusted serial result.

Restore the original program when you finish.

Practice: can this loop be parallelized?

For each loop, decide whether adding # pragma omp parallel for is safe as written. If not, can a clause or a small rewrite fix it? Assume the arrays do not overlap.

/* A */  for (i = 0; i < n; i++)  a[i] = b[i] + c[i];
/* B */  for (i = 1; i < n; i++)  a[i] = a[i - 1] + b[i];
/* C */  for (i = 0; i < n; i++)  total += a[i];
/* D */  for (i = 0; i < n - 1; i++)  a[i] = 2 * a[i + 1];
/* E */  for (i = 0; i < n; i++)  { t = a[i] * a[i];  b[i] = t + 1; }
/* F */  for (i = 0; i < n; i++)  if (a[i] > biggest) biggest = a[i];
  • A: Safe. Each iteration reads and writes only its own elements.
  • B: Not safe. Iteration i reads a[i - 1], written by the previous iteration. This running total is a loop-carried dependence like Fibonacci’s.
  • C: Not safe as written: every iteration updates the shared total. Add reduction(+: total).
  • D: Not safe. Iteration i reads a[i + 1], which iteration i + 1 overwrites. The result depends on which happens first. Writing the results to a second array removes the dependence.
  • E: Not safe as written, because a shared t is a scope bug. Add private(t), or declare t inside the loop.
  • F: Not safe as written: every thread may read and update biggest. Add reduction(max: biggest).

Check your understanding

  1. What does parallel for do that parallel alone does not?
  2. Why must the loop variable be private? Do you have to declare it private yourself?
  3. Why does OpenMP refuse to parallelize a loop that might break?
  4. What is a loop-carried dependence? Does the compiler warn you about one?
  5. Why is x[i] = a + i * h; y[i] = exp(x[i]); safe, even though the second statement depends on the first?
  6. What does default(none) do, and why is it a good habit?
  7. What is the value of a private variable when a thread starts its first iteration?
  1. It also divides the iterations of the following for loop among the threads. With parallel alone, every thread runs the whole block, and we divide the work ourselves.
  2. If one i were shared, every thread’s i++ would race with the others. OpenMP makes the loop variable of a parallel for private automatically.
  3. OpenMP must know the number of iterations before the loop starts so it can divide them. An early exit also leaves the structured block from the middle.
  4. One iteration uses a value written by a different iteration. No: the program compiles and may give wrong answers only sometimes.
  5. Both statements are in the same iteration, which always runs on one thread, in order.
  6. It makes the compiler require a scope for every outside variable used in the block. A forgotten private then becomes a compile error, not a race.
  7. Unspecified (uninitialized). It does not inherit the value from before the loop.
Back to top