Memory Hierarchy, Warps, and Warp Shuffles

In the first CUDA trapezoidal implementation, 1,000,000 threads simultaneously attempted atomic additions on a single memory address, causing massive hardware serialization.

To achieve maximum performance on GPUs, we must combine tree-structured reductions with knowledge of the CUDA memory hierarchy and warp-level primitives.

---

6.16 Serialized vs. Tree-Structured Global Sums

Consider an 8-thread sum where each thread adds its partial trapezoid value (my_trap) into *trap_p via atomicAdd:

Table 6.6: Serialized Global Sum with Eight Threads

Time Step	Thread Winning Access	Partial Value (my_trap)	Global Sum (*trap_p)
Start	—	—	$9$
$t_0$	Thread 5	$11$	$20$
$t_1$	Thread 2	$5$	$25$
$t_2$	Thread 3	$7$	$32$
$t_3$	Thread 7	$15$	$47$
$t_4$	Thread 4	$9$	$56$
$t_5$	Thread 6	$13$	$69$
$t_6$	Thread 0	$1$	$70$
$t_7$	Thread 1	$3$	$73$

This serialized approach requires $t$ consecutive sequential operations.

The Tree-Structured Reduction

By pairing active threads in parallel rounds, the computation shrinks from $t$ sequential additions to $\lceil \log_2(t) \rceil$ phases:

flowchart TD
  subgraph TreeRed["Figure 6.4: Tree-Structured Reduction (8 Threads)"]
      direction TB
      subgraph Initial["Initial Values"]
          T0_0["T0: 1"]
          T1_0["T1: 3"]
          T2_0["T2: 5"]
          T3_0["T3: 7"]
          T4_0["T4: 9"]
          T5_0["T5: 11"]
          T6_0["T6: 13"]
          T7_0["T7: 15"]
      end
      subgraph Stage1["Stage 1: 4 Parallel Adds"]
          T0_1["T0: 1 + 9 = 10"]
          T1_1["T1: 3 + 11 = 14"]
          T2_1["T2: 5 + 13 = 18"]
          T3_1["T3: 7 + 15 = 22"]
      end
      subgraph Stage2["Stage 2: 2 Parallel Adds"]
          T0_2["T0: 10 + 18 = 28"]
          T1_2["T1: 14 + 22 = 36"]
      end
      subgraph Stage3["Stage 3: Final Warp Sum"]
          T0_3["T0: 28 + 36 = 64"]
      end
      
      T0_0 --> T0_1
      T4_0 -.-> T0_1
      T1_0 --> T1_1
      T5_0 -.-> T1_1
      T2_0 --> T2_1
      T6_0 -.-> T2_1
      T3_0 --> T3_1
      T7_0 -.-> T3_1
      
      T0_1 --> T0_2
      T2_1 -.-> T0_2
      T1_1 --> T1_2
      T3_1 -.-> T1_2
      
      T0_2 --> T0_3
      T1_2 -.-> T0_3
  end

For $1,000,000$ values, a tree reduction drops execution from $1,000,000$ sequential steps to just $20$ parallel steps!

---

6.17 The CUDA Memory Hierarchy

Performance in CUDA is dictated by where operands are stored:

flowchart TD
  subgraph MemHier["Figure 6.1B: CUDA Memory Hierarchy"]
      direction TB
      REG["Registers: ~1 cycle latency
Private per thread (fastest)"]
      SMEM["Shared Memory: ~a few cycles latency
On-chip, shared per block"]
      GMEM["Global Memory: 200 - 400 cycles latency
Off-chip DRAM, accessible across entire grid"]
      REG <==> SMEM
      SMEM <==> GMEM
  end

Table 6.7: Memory Statistics for Selected Nvidia GPUs

GPU Model	Compute Capability	Registers per Thread	Shared Memory per Block	Global VRAM per GPU
Quadro 600	2.1	504 Bytes	48 KB	1 GB
GK20A (Jetson TK1)	3.2	504 Bytes	48 KB	2 GB
GeForce GTX Titan X	5.2	504 Bytes	48 KB	12 GB

• Registers: Reside inside each Streaming Processor (SP). Accessing a register takes $\approx 1$ clock cycle. If a kernel uses too many local variables, they “spill” into slow local memory (backed by off-chip DRAM).
• Shared Memory: Fast, user-managed on-chip SRAM shared exclusively among threads in the same thread block. Access takes only a few cycles.
• Global Memory: Large off-chip GDDR/HBM DRAM. Accessible by all threads in the grid, but incurs high latency ($200 - 400$ cycles).

---

6.18 Warps and Warp Shuffles

In Nvidia hardware, threads within a block are scheduled and executed in groups of 32 consecutive threads called a warp:

$\text{warpSize} = 32$

A thread’s position within its warp is called its lane:

$\text{lane} = \text{threadIdx.x} \pmod{\text{warpSize}}$

Threads within a warp execute instructions in SIMT lockstep. When all threads follow the same path, execution is maximally efficient.

Warp Shuffle: __shfl_down_sync

Introduced in Kepler (Compute Capability $\ge 3.0$) and modernized in CUDA 9, warp shuffles permit threads within a warp to read registers directly from other threads in the same warp without using shared memory or memory buses:

__device__ float __shfl_down_sync(
    unsigned mask,             /* in: active thread bitmask (0xffffffff) */
    float    var,              /* in: register variable to share        */
    unsigned diff,             /* in: lane offset to read from           */
    int      width = warpSize  /* in: sub-warp size (default 32)         */
);

• mask = 0xffffffff: A 32-bit hex mask with all 32 bits set to 1, asserting that all lanes in the warp participate.
• A calling thread with lane $l$ receives the value of var held by the thread with lane: $\text{target\_lane} = l + \text{diff}$
• If $l + \text{diff} \ge \text{warpSize}$, the function returns the calling thread’s own var.

---

Implementing Intra-Warp Sum: Warp_sum

Using __shfl_down_sync, a warp of 32 threads performs a tree reduction in just 5 steps ($\text{diff} = 16, 8, 4, 2, 1$):

__device__ float Warp_sum(float var) {
    unsigned mask = 0xffffffff;

    for (int diff = warpSize / 2; diff > 0; diff = diff / 2) {
        var += __shfl_down_sync(mask, var, diff);
    }

    return var; /* Result is mathematically correct on Lane 0 */
}

flowchart TD
  subgraph ShuffleFlow["Figure 6.5: Warp Shuffle Down (diff = 4, 2, 1 on 8 lanes)"]
      direction TB
      subgraph Step1["diff = 4: Lanes 0-3 receive from Lanes 4-7"]
          S1["Lanes 0..3 add values from lanes 4..7
Lanes 4..7 retain their values"]
      end
      subgraph Step2["diff = 2: Lanes 0-1 receive from Lanes 2-3"]
          S2["Lanes 0..1 add values from lanes 2..3"]
      end
      subgraph Step3["diff = 1: Lane 0 receives from Lane 1"]
          S3["Lane 0 adds value from lane 1
Lane 0 now holds total warp sum!"]
      end
      Step1 --> Step2 --> Step3
  end

---

6.19 Alternative: Shared Memory Dissemination Sum

For architectures without warp shuffles (or as a general reduction primitive), threads can exchange values through on-chip __shared__ memory:

__device__ float Shared_mem_sum(float shared_vals[]) {
    int my_lane = threadIdx.x % warpSize;

    for (int diff = warpSize / 2; diff > 0; diff = diff / 2) {
        int source = (my_lane + diff) % warpSize;
        shared_vals[my_lane] += shared_vals[source];
    }

    return shared_vals[my_lane];
}

Because threads in a single warp operate in lockstep, all threads read shared_vals[source] before any thread overwrites shared_vals[my_lane], preventing race conditions without requiring barriers.

---

6.20 Trapezoidal Rule with 32-Thread Blocks

We restructure the trapezoidal kernel so that each block contains exactly warpSize = 32 threads:

1. Each thread computes its local trapezoid area.
2. The warp performs a fast intra-warp reduction via Warp_sum() or Shared_mem_sum().
3. Only Lane 0 of each block executes atomicAdd into the global total!

Program 6.13: Kernel with Warp Shuffle

/* Program 6.13: Trapezoidal rule with Warp_sum */
__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;
    float my_trap = 0.0f;

    if (0 < my_i && my_i < n) {
        float my_x = a + my_i * h;
        my_trap = f(my_x);
    }

    /* 32 threads reduce their values in registers */
    float result = Warp_sum(my_trap);

    /* Only Thread 0 in the block calls atomicAdd */
    if (threadIdx.x == 0) {
        atomicAdd(trap_p, result);
    }
}

Program 6.14: Kernel with Shared Memory

/* Program 6.14: Trapezoidal rule with Shared_mem_sum */
#define WARPSZ 32

__global__ void Dev_trap(
    const float a,
    const float b,
    const float h,
    const int   n,
    float*      trap_p
) {
    __shared__ float shared_vals[WARPSZ];
    int my_i = blockDim.x * blockIdx.x + threadIdx.x;
    int my_lane = threadIdx.x % warpSize;

    shared_vals[my_lane] = 0.0f;
    if (0 < my_i && my_i < n) {
        float my_x = a + my_i * h;
        shared_vals[my_lane] = f(my_x);
    }

    float result = Shared_mem_sum(shared_vals);

    if (threadIdx.x == 0) {
        atomicAdd(trap_p, result);
    }
}

---

6.21 Benchmark Comparison: Naive vs. Warp Reduction

Integrating $f(x) = x^2 + 1$ with $n = 1,048,576$ trapezoids across 32,768 blocks of 32 threads:

Table 6.8: Mean Run-times for Trapezoidal Rule Using Block Size of 32 Threads (ms)

Implementation	Nvidia GK20A (1 SM, 192 SPs)	Nvidia GeForce GTX Titan X (24 SMs, 3072 SPs)
Original Naive atomicAdd	$20.7\,\text{ms}$	$3.08\,\text{ms}$
Warp Shuffle (Warp_sum)	$14.4\,\text{ms}$	$0.210\,\text{ms}$ ($14.6\times$ Faster!)
Shared Memory (Shared_mem_sum)	$15.0\,\text{ms}$	$0.206\,\text{ms}$ ($15.0\times$ Faster!)

By reducing the number of global atomic collisions by a factor of 32 (from $1,048,576$ calls down to $32,768$), runtime on the Titan X plummeted from $3.08\,\text{ms}$ to just $0.21\,\text{ms}$!