Multi-Warp Reductions and Parallel Bitonic Sort

Restricting thread blocks to 32 threads underutilizes the GPU’s hardware occupancy. Modern GPUs support blocks with up to 1024 threads (32 warps).

In this chapter, we scale reductions across multiple warps using __syncthreads(), eliminate shared memory bank conflicts, and implement a high-performance parallel bitonic sort across multi-block grids.

---

6.22 Multi-Warp Block Reductions and __syncthreads()

When a thread block contains multiple warps (e.g., 1024 threads = 32 warps), warps execute independently. If Warp 0 attempts to aggregate the partial results of the other warps before they finish, a race condition occurs.

CUDA provides an intra-block barrier:

__device__ void __syncthreads(void);

__syncthreads() forces all threads in the thread block to wait until every thread in that block reaches the barrier.

Barriers and Conditional Code

All threads in a block must execute __syncthreads(). Placing __syncthreads() inside a conditional branch where only some threads enter (e.g., if (threadIdx.x < 16) __syncthreads();) results in hardware deadlock. Furthermore, __syncthreads() synchronizes threads only within the same block, never across the entire grid.

---

6.23 Shared Memory Banks and Bank Conflicts

Shared memory is partitioned into 32 independent memory banks organized so that consecutive 32-bit (4-byte) words map to successive banks:

$\text{Bank Number} = \left( \frac{\text{Byte Address}}{4} \right) \pmod{32}$

Table 6.9: Shared Memory Bank Mapping (thread_calcs float indices)

Memory Bank	Bank 0	Bank 1	Bank 2	…	Bank 30	Bank 31
Row 0	0	1	2	…	30	31
Row 1	32	33	34	…	62	63
Row 2	64	65	66	…	94	95
Row 3	96	97	98	…	126	127
Row 31	992	993	994	…	1022	1023

Bank Access Rules:

• Conflict-Free: If all 32 threads in a warp access different banks (e.g., elements 0..31), requests are served concurrently in a single clock cycle.
• Broadcast: If multiple threads read the exact same memory address, the value is broadcast in a single cycle.
• Bank Conflict: If multiple threads access different addresses within the same bank (e.g., threads accessing elements 0, 32, 64), requests are serialized, degrading throughput by up to $32\times$.

---

6.24 Complete Multi-Warp Kernel: Program 6.15

In this optimized implementation:

1. Each warp performs an intra-warp reduction using Shared_mem_sum.
2. Lane 0 of each warp stores its warp sum into warp_sum_arr[my_warp] (contiguous indices 0..31, preventing bank conflicts!).
3. __syncthreads() ensures all 32 warps have written their results.
4. Warp 0 performs a final reduction across warp_sum_arr.
5. Only Thread 0 of the entire block executes a single atomicAdd into global memory.

/* Program 6.15: Trapezoidal rule with multi-warp shared memory */
#define MAX_BLKSZ 1024
#define WARPSZ 32

__global__ void Dev_trap(
    const float a,
    const float b,
    const float h,
    const int   n,
    float*      trap_p
) {
    __shared__ float thread_calcs[MAX_BLKSZ];
    __shared__ float warp_sum_arr[WARPSZ];

    int my_i    = blockDim.x * blockIdx.x + threadIdx.x;
    int my_warp = threadIdx.x / warpSize;
    int my_lane = threadIdx.x % warpSize;

    float* shared_vals = thread_calcs + my_warp * warpSize;
    float  blk_result  = 0.0f;

    /* 1. Calculate trapezoid area */
    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);
    }

    /* 2. Intra-warp reduction */
    float my_result = Shared_mem_sum(shared_vals);
    if (my_lane == 0) {
        warp_sum_arr[my_warp] = my_result;
    }

    /* 3. Block synchronization barrier */
    __syncthreads();

    /* 4. Warp 0 aggregates all warp sums in the block */
    if (my_warp == 0) {
        if (threadIdx.x >= blockDim.x / warpSize) {
            warp_sum_arr[threadIdx.x] = 0.0f;
        }
        blk_result = Shared_mem_sum(warp_sum_arr);
    }

    /* 5. One atomicAdd per block */
    if (threadIdx.x == 0) {
        atomicAdd(trap_p, blk_result);
    }
}

Table 6.10: Mean Run-times Across Optimization Stages ($n = 1,048,576$)

Version	Description	Nvidia GK20A	Nvidia GTX Titan X
Original	Naive global atomicAdd	$20.7\,\text{ms}$	$3.08\,\text{ms}$
32-Thread Blocks	1 Warp per Block	$14.4\,\text{ms}$	$0.210\,\text{ms}$
1024-Thread Blocks	Multi-Warp Reduction	$12.8\,\text{ms}$	$0.141\,\text{ms}$ ($21.8\times$ Faster!)

---

6.25 Parallel Bitonic Sort on the GPU

Sorting arbitrary arrays on a GPU requires an algorithm with fixed, predictable comparison patterns. Bitonic Sort is ideal because it decouples comparisons into independent, parallel butterfly operations.

A bitonic sequence is a sequence of keys that monotonically increases and then monotonically decreases (or can be cyclically shifted to do so).

flowchart LR
  subgraph Butterflies["Figures 6.7 & 6.8: Butterfly Compare-Swap Building Blocks"]
      direction TB
      subgraph TwoElt["2-Element Butterfly"]
          A["(a0, a1)"] -->|Compare & Swap| B["Increasing: (min, max)
Decreasing: (max, min)"]
      end
      subgraph FourElt["4-Element Butterfly"]
          C["(a0, a2) and (a1, a3)"] -->|Stage A: Halve| D["(a0, a1) and (a2, a3)"]
          D -->|Stage B: Sort Halves| E["Globally Sorted 4-Element List"]
      end
  end

Table 6.11: Bitonic Sort of a 4-Element List ({40, 20, 10, 30})

Step	Subscript 0	Subscript 1	Subscript 2	Subscript 3
Initial List	$40$	$20$	$10$	$30$
2-Element Butterfly	Compare-Swap (Incr) $\to 20$	$40$	Compare-Swap (Decr) $\to 30$	$10$
Bitonic Sequence	$20$	$40$	$30$	$10$
4-Element Stage A	Compare $(20, 30) \to 20$	Compare $(40, 10) \to 10$	$30$	$40$
4-Element Stage B	Compare $(20, 10) \to \mathbf{10}$	$\mathbf{20}$	Compare $(30, 40) \to \mathbf{30}$	$\mathbf{40}$

---

Bit Manipulation in Bitonic Sort

By analyzing binary indices, a thread can determine its paired comparison element using bitwise XOR (^):

$\text{my\_elt2} = \text{my\_elt1} \oplus \text{stage}$

Whether elements are sorted in increasing or decreasing order is determined by checking:

$\text{is\_increasing} = (\text{my\_elt1} \ \& \ \text{bf\_sz}) == 0$

/* Program 6.17: Serial Bitonic Sort Kernel Logic */
for (bf_sz = 2; bf_sz <= n; bf_sz *= 2) {
    for (stage = bf_sz / 2; stage > 0; stage /= 2) {
        which_bit = Which_bit(stage);
        for (th = 0; th < n / 2; th++) {
            my_elt1 = Insert_zero(th, which_bit);
            my_elt2 = my_elt1 ^ stage;
            Compare_swap(list, my_elt1, my_elt2, my_elt1 & bf_sz);
        }
    }
}

---

6.26 Parallel Bitonic Sort across a Grid

Because current GPU thread blocks are limited to 1024 threads, a single block can sort at most $2 \times 1024 = 2048$ elements.

To sort millions of elements across a multi-block grid, we overcome the block synchronization limit through host-coordinated kernel relaunching:

/* Program 6.19: Multi-block Bitonic Sort Host Driver */
void Parallel_bitonic_sort(float list[], int n, int blk_ct, int th_per_blk) {
    /* Step 1: Each block independently sorts sublists of size 2 * th_per_blk */
    Pbitonic_start<<<blk_ct, th_per_blk>>>(list, n);

    /* Step 2: Global butterflies spanning multiple blocks */
    for (int bf_sz = 4 * th_per_blk; bf_sz <= n; bf_sz *= 2) {
        /* Stages requiring inter-block communication are launched as separate kernels */
        for (int stage = bf_sz / 2; stage >= 2 * th_per_blk; stage /= 2) {
            Pbutterfly_one_stage<<<blk_ct, th_per_blk>>>(list, n, bf_sz, stage);
            /* Implicit grid barrier occurs between kernel launches */
        }
        /* Finish remaining intra-block stages with __syncthreads */
        Pbutterfly_finish<<<blk_ct, th_per_blk>>>(list, n, bf_sz);
    }
}

Kernel completions act as implicit grid-wide barriers, enabling safe cross-block data exchanges without specialized hardware locks.

---

6.27 Empirical Benchmarks: GPU vs. CPU Sorting

We compare serial CPU sorting algorithms against CUDA Parallel Bitonic Sort:

Table 6.18: Small List Sorting Performance ($n = 2048$ integers)

System	CPU: Intel Xeon 4116 (qsort)	CPU: Intel Xeon 4116 (Bitonic)	GPU: Nvidia Quadro P5000 (1 Block, 1024 Threads)
Run-time	$0.116\,\text{ms}$	$1.19\,\text{ms}$	$0.128\,\text{ms}$

For tiny datasets ($n = 2048$), kernel launch overhead neutralizes GPU parallelism.

Table 6.19: Large List Sorting Performance ($n = 2^{21} = 2,097,152$ integers)

System	CPU: Intel Xeon 4116 (qsort)	CPU: Intel Xeon 4116 (Bitonic)	GPU: Nvidia Quadro P5000 (1024 Blocks, 1024 Threads)
Run-time	$539.0\,\text{ms}$	$3710.0\,\text{ms}$	$13.1\,\text{ms}$ ($41.1\times$ Faster than CPU qsort!)

At scale, the GPU finishes in $13.1\,\text{ms}$, delivering an astounding $41\times$ speedup over optimized sequential C qsort.

---

6.28 Chapter Summary

• SIMT Execution Model: Hardware threads execute in 32-wide warps. Latency is hidden by scheduling ready warps while stalled warps wait for memory.
• Unified Memory (cudaMallocManaged): Automatically migrates memory pages between host CPU and device GPU on demand.
• Warp Shuffles (__shfl_down_sync): Ultra-fast intra-warp tree reductions executed directly within registers with zero memory overhead.
• Shared Memory Banks: 32 banks provide parallel access; conflict-free indexing patterns must be preserved to prevent memory stalls.
• Block Barriers (__syncthreads): Coordinates warps within a thread block. Global grid barriers are achieved by relaunching kernels from the host.
• Massive Parallel Speedup: Scaled algorithms (such as multi-block Bitonic Sort) surpass optimized CPU baselines by over $40\times$.