The Trapezoidal Rule and the Reduction Clause

Having examined basic thread creation, we now apply OpenMP to numerical integration using the trapezoidal rule. This highlights fundamental shared-memory concepts: race conditions, critical sections, variable scope, and the reduction clause.

---

5.3 Numerical Integration via the Trapezoidal Rule

Recall that the trapezoidal rule approximates the definite integral $\int_a^b f(x) \, dx$ by dividing the interval $[a, b]$ into $n$ equal subintervals of width $h$:

$h = \frac{b - a}{n}, \quad x_i = a + i \cdot h \quad (0 \le i \le n)$

$\text{Area} \approx h \left[ \frac{f(x_0)}{2} + f(x_1) + f(x_2) + \dots + f(x_{n-1}) + \frac{f(x_n)}{2} \right]$

Serial Algorithm

/* Serial Trapezoidal Rule */
h = (b - a) / n;
approx = (f(a) + f(b)) / 2.0;
for (i = 1; i <= n - 1; i++) {
    x_i = a + i * h;
    approx += f(x_i);
}
approx = h * approx;

---

5.4 A First OpenMP Version: Race Conditions and Critical Sections

In a shared-memory environment, we partition the $n$ trapezoids across thread_count threads. Assuming $n$ is evenly divisible by thread_count, each thread is assigned a contiguous block of:

$\text{local\_n} = \frac{n}{\text{thread\_count}}$

flowchart TD
  subgraph TrapPart["Figure 5.4: Partitioning Trapezoids Among 4 Threads"]
      direction LR
      subgraph T0["Thread 0"]
          S0["Interval: [a, a + local_n * h]
Compute local sum"]
      end
      subgraph T1["Thread 1"]
          S1["Interval: [local_a, local_b]
Compute local sum"]
      end
      subgraph T2["Thread 2"]
          S2["Interval: [local_a, local_b]
Compute local sum"]
      end
      subgraph T3["Thread 3"]
          S3["Interval: [local_a, b]
Compute local sum"]
      end
      
      G["Shared Variable: global_result"]
      
      T0 -->|Accumulate| G
      T1 -->|Accumulate| G
      T2 -->|Accumulate| G
      T3 -->|Accumulate| G
  end

Each thread computes its private contribution (my_result). To compute the total integral, each thread must add its result to a shared accumulator variable:

global_result += my_result;

---

The Race Condition

In assembly language, the operation global_result += my_result is not atomic; it consists of three distinct machine steps:

1. Load global_result from main memory into a CPU register.
2. Add my_result to the register.
3. Store the updated value from the register back into global_result in memory.

If two threads execute this simultaneously, their operations can interleave:

Table 5.0: Race Condition Interleaving Timeline

Time	Thread 0	Thread 1	Value of global_result in Memory
0	Reads global_result ($0$) into register	Computes my_result ($2$)	$0$
1	Adds my_result ($1$) to register $\to 1$	Reads global_result ($0$) into register	$0$
2	Writes $1$ back to global_result	Adds my_result ($2$) to register $\to 2$	1
3	Continues	Writes $2$ back to global_result	2 (Error: lost update!)

Because Thread 1 read global_result before Thread 0 finished writing, Thread 0’s computation was overwritten. The final result is $2$ instead of the correct sum $1 + 2 = 3$.

This error is a race condition: multiple concurrent threads access shared memory, at least one access is a write, and the outcome depends on the nondeterministic order of execution.

---

Mutex Protection with #pragma omp critical

A block of code that accesses a shared resource and must be executed by only one thread at a time is called a critical section.

OpenMP provides the #pragma omp critical directive to enforce mutual exclusion:

#pragma omp critical
global_result += my_result;

When a thread encounters a critical directive, it waits until no other thread is executing within the protected block.

---

Complete Implementation: Program 5.2

Below is the complete source code for Program 5.2 (omp_trap1.c):

/* Program 5.2: First OpenMP trapezoidal rule program */
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

void Trap(double a, double b, int n, double* global_result_p);

/* Function to integrate: f(x) = x * x */
double f(double x) {
    return x * x;
}

int main(int argc, char* argv[]) {
    double global_result = 0.0;
    double a, b;
    int n;
    int thread_count;

    thread_count = strtol(argv[1], NULL, 10);
    printf("Enter a, b, and n\n");
    scanf("%lf %lf %d", &a, &b, &n);

    if (n % thread_count != 0) {
        fprintf(stderr, "Error: n must be evenly divisible by thread_count\n");
        exit(1);
    }

    #pragma omp parallel num_threads(thread_count)
    Trap(a, b, n, &global_result);

    printf("With n = %d trapezoids, our estimate\n", n);
    printf("of the integral from %f to %f = %.14e\n", a, b, global_result);

    return 0;
} /* main */

void Trap(double a, double b, int n, double* global_result_p) {
    double h, x, my_result;
    double local_a, local_b;
    int i, local_n;
    int my_rank = omp_get_thread_num();
    int thread_count = omp_get_num_threads();

    h = (b - a) / n;
    local_n = n / thread_count;
    local_a = a + my_rank * local_n * h;
    local_b = local_a + local_n * h;

    my_result = (f(local_a) + f(local_b)) / 2.0;
    for (i = 1; i <= local_n - 1; i++) {
        x = local_a + i * h;
        my_result += f(x);
    }
    my_result = my_result * h;

    #pragma omp critical
    *global_result_p += my_result;
} /* Trap */

---

5.5 Scope of Variables: Shared vs. Private

In OpenMP, variable scope determines which threads can access a given memory variable inside a parallel region:

• Shared Scope: A single memory location is accessible by all threads in the team.
  • Rule: Variables declared before a parallel region default to shared scope (e.g., a, b, n, global_result, thread_count).
  • Any modification by one thread is visible to all other threads.
• Private Scope: Each thread possesses its own distinct instance of the variable.
  • Rule: Variables declared inside the parallel block or inside functions called from within the block reside on the thread’s private stack (e.g., local_a, local_b, local_n, my_result, my_rank).
  • Changes to a private variable by one thread do not affect any other thread.

Pointers to Shared Variables

In Trap(..., &global_result), the parameter global_result_p is a private pointer variable on each thread’s stack, but it points to the memory address of global_result, which is shared across all threads. Therefore, dereferencing *global_result_p += my_result modifies shared memory and requires #pragma omp critical.

---

5.6 The Reduction Clause

Suppose we want to avoid passing pointers and write a cleaner function Local_trap(a, b, n) that returns a double.

A naive attempt might look like this:

/* SEVERE PERFORMANCE BUG: Accidental Serialization */
global_result = 0.0;
#pragma omp parallel num_threads(thread_count)
{
    #pragma omp critical
    global_result += Local_trap(a, b, n);
}

Serializing Parallel Execution

Because the call to Local_trap is placed directly inside #pragma omp critical, only one thread at a time can calculate its subinterval! The entire trapezoidal computation is executed sequentially, and thread synchronization overhead makes it run slower than single-threaded code.

To fix this, the computation must remain outside the critical section:

global_result = 0.0;
#pragma omp parallel num_threads(thread_count)
{
    double my_result = 0.0; /* private */
    my_result += Local_trap(a, b, n);

    #pragma omp critical
    global_result += my_result;
}

---

OpenMP Reductions

OpenMP provides a declarative mechanism to automate this pattern: the reduction clause.

A reduction repeatedly applies an associative binary operator (e.g., addition, multiplication) to a collection of values, accumulating the result into a single variable:

$\text{Syntax:} \quad \text{reduction}(<\text{operator}> : <\text{variable list}>)$

global_result = 0.0;
#pragma omp parallel num_threads(thread_count) \
    reduction(+: global_result)
global_result += Local_trap(a, b, n);

How the Runtime Executes Reductions

1. OpenMP creates a private copy of global_result for each thread in the team.
2. Each private copy is automatically initialized to the identity value of the specified operator (e.g., 0 for addition, 1 for multiplication).
3. Threads accumulate into their private copies completely independently and in parallel without locking.
4. At the end of the parallel block, OpenMP automatically aggregates the private variables into the shared global_result inside an internal critical section.

---

Reduction Operators and Identity Values

Table 5.1: Identity Values for OpenMP Reduction Operators

Operator	Mathematical Operation	Initial Identity Value
+	Addition	0
*	Multiplication	1
-	Subtraction	0 (Partial results added across threads)
&	Bitwise AND	~0 (All bits set to 1)
|	Bitwise OR	0
^	Bitwise XOR	0
&&	Logical AND	1
||	Logical OR	0

Floating-Point Non-Associativity

While mathematical addition is associative, floating-point addition in IEEE 754 arithmetic is not strictly associative due to roundoff error: $(a + b) + c \neq a + (b + c)$. Therefore, changing the number of threads in an OpenMP reduction may cause minor differences in the least-significant bits of float or double results.