The Trapezoidal Rule and CUDA Atomic Operations

Having mastered vector addition and CUDA memory management, we apply GPU parallelism to numerical integration via the trapezoidal rule.

This problem exposes a central challenge in massively parallel computing: aggregating partial results into a single global total without causing catastrophic serialization bottlenecks.

---

6.11 Numerical Integration on the GPU

The trapezoidal rule approximates $\int_a^b f(x) \, dx$ across $n$ trapezoids of uniform width $h$:

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

$\text{Area} \approx h \left[ \dfrac{f(a) + f(b)}{2} + \sum_{i=1}^{n-1} f(x_i) \right]$

Serial Reference Implementation: Program 6.11

/* Program 6.11: Serial Trapezoidal Rule on CPU */
float Serial_trap(const float a, const float b, const int n) {
    float x, h = (b - a) / n;
    float trap = 0.5f * (f(a) + f(b));

    for (int i = 1; i <= n - 1; i++) {
        x = a + i * h;
        trap += f(x);
    }
    trap = trap * h;
    return trap;
}

---

6.12 Designing the CUDA Trapezoidal Kernel

When $n$ is large, over 99.9% of the computational work occurs inside the loop computing $f(x_i)$. We assign each iteration $i \in [1, n-1]$ to an individual CUDA thread:

$\text{my\_i} = \text{blockDim.x} \times \text{blockIdx.x} + \text{threadIdx.x}$

$\text{my\_x} = a + \text{my\_i} \times h, \quad \text{my\_trap} = f(\text{my\_x})$

Five Implementation Challenges:

1. Parameter Privacy: In CUDA, function arguments passed to a kernel reside on the executing thread’s private stack. A thread cannot simply initialize a local variable and have it visible to other threads.
2. Boundary Filtering: The loop only sums interior points ($1 \le i \le n-1$). Threads with $\text{my\_i} = 0$ or $\text{my\_i} \ge n$ must not participate in the summation.
3. Accumulation Race Condition: Multiple threads writing to *trap_p += my_trap simultaneously produce undefined race conditions.
4. Return Value: Kernels return void; the final sum must be written to global/managed memory.
5. Final Scaling: The aggregated sum must be multiplied by $h$ once all thread additions complete.

---

6.13 Hardware Atomics: atomicAdd()

To prevent race conditions without software locks, CUDA provides hardware-accelerated atomic memory operations.

An operation is atomic if it executes as an indivisible unit: while one thread is reading, modifying, and writing the memory location, all other threads are locked out of that address.

__device__ float atomicAdd(
    float* address  /* in/out: memory location to update */,
    float  val      /* in: value to add                  */
);

atomicAdd adds val to *address, stores the sum back into *address, and returns the original value stored in memory prior to the addition.

---

6.14 Complete Implementation: Program 6.12

To cleanly separate setup and kernel invocation, we encapsulate execution inside a host wrapper function (Trap_wrapper):

/* Program 6.12: First CUDA Trapezoidal Rule Implementation */
#include <stdio.h>
#include <cuda.h>

/* Mathematical target function: f(x) = x^2 + 1 */
__device__ float f(float x) {
    return x * x + 1.0f;
}

/* CUDA Kernel: Executed by thousands of GPU threads */
__global__ void Dev_trap(
    const float a,
    const float b,
    const float h,
    const int   n,
    float*      trap_p
) {
    int my_i = blockDim.x * blockIdx.x + threadIdx.x;

    /* Compute f(x_1), f(x_2), ..., f(x_{n-1}) */
    if (0 < my_i && my_i < n) {
        float my_x = a + my_i * h;
        float my_trap = f(my_x);
        atomicAdd(trap_p, my_trap);
    }
} /* Dev_trap */

/* Host Wrapper Function */
void Trap_wrapper(
    const float a,
    const float b,
    const int   n,
    float*      trap_p,
    const int   blk_ct,
    const int   th_per_blk
) {
    /* Initialize trap_p with endpoint contributions */
    *trap_p = 0.5f * (f(a) + f(b));
    float h = (b - a) / n;

    /* Launch kernel */
    Dev_trap<<<blk_ct, th_per_blk>>>(a, b, h, n, trap_p);
    cudaDeviceSynchronize();

    /* Final scaling on host CPU */
    *trap_p = h * (*trap_p);
} /* Trap_wrapper */

---

6.15 Empirical Performance and the Atomic Bottleneck

To benchmark the code, we integrate $f(x) = x^2 + 1$ over $[-3, 3]$ with $n = 2^{20} = 1,048,576$ trapezoids using 1024 blocks of 1024 threads ($1,048,576$ threads):

Table 6.5: Mean Run-times for Serial CPU vs. First CUDA Trapezoidal Rule ($n = 1,048,576$)

System	CPU: ARM Cortex-A15	GPU: Nvidia GK20A	CPU: Intel Core i7	GPU: Nvidia GeForce GTX Titan X
Clock Frequency	2.3 GHz	852 MHz	3.5 GHz	1.08 GHz
Cores / SMs / SPs	4 Cores	1 SM, 192 SPs	4 Cores	24 SMs, 3072 SPs
Execution Time	$33.6\,\text{ms}$	$20.7\,\text{ms}$	$4.48\,\text{ms}$	$3.08\,\text{ms}$

---

Diagnosing the Performance Bottleneck

Why did a top-tier Nvidia GeForce GTX Titan X with 3072 CUDA Cores run only $1.45\times$ faster than a single CPU core?

flowchart TD
  subgraph Bottleneck["The Global Atomic Bottleneck"]
      direction TB
      T0["Thread 0"]
      T1["Thread 1"]
      T2["Thread 2"]
      TK["Thread 1,048,575"]
      
      MEM["Single Memory Location: *trap_p
(Global Memory DRAM)"]
      
      T0 ==>|atomicAdd| MEM
      T1 -.->|Blocked / Queued| MEM
      T2 -.->|Blocked / Queued| MEM
      TK -.->|Blocked / Queued| MEM
  end

• Massive Serialization: Over 1,000,000 threads are attempting to perform atomicAdd on the exact same memory location in high-latency global memory!
• While atomicAdd guarantees mathematical correctness, the GPU memory controller must serialize all 1,000,000 writes one by one.
• Instead of running in parallel, the accumulation phase becomes purely sequential.

In the next section, we explore tree-structured reductions and warp shuffles to eliminate this bottleneck.