CUDA Parallel Execution and Optimization

Writing CUDA code is easy, but making it run faster than a CPU is hard. Without a structured workflow, developers often optimize the wrong code or introduce subtle bugs.

This post covers APOD (Assess $\rightarrow$ Parallelize $\rightarrow$ Optimize $\rightarrow$ Deploy), NVIDIA’s standard optimization cycle, and applies it step-by-step to Parallel Reduction.

Table of contents

Quick Warm-up 🧠

Let’s review key concepts from the previous lecture:

Question 1: Block to Warps

A block has 256 threads. How many warps does the GPU hardware create?

Question 2: The Divergence Trap

How many warps diverge in a 256-thread block with this branch?

if (threadIdx.x % 2 == 0)
    a();
else
    b();

Question 3: Latency Hiding

Why does an SM need many resident warps?

Part 1: The APOD Engineering Process 🔄

Professional GPU optimization is not guesswork—it is a repeatable engineering loop.

1. Why Do We Need a Process?

Developers rewriting CPU loops in CUDA often hit common traps:

• The GPU version runs slower than single-threaded CPU code.
• Days are spent optimizing a kernel that takes only 2% of total runtime.
• Fast code produces wrong results on edge-case inputs.

These stem from four mistakes:

1. Guessing instead of measuring.
2. Optimizing the wrong bottleneck.
3. Changing too many variables at once.
4. Skipping correctness checks.

2. The 4 Stages of the APOD Cycle

Figure 1: The cyclic APOD optimization workflow.

• Assess: Profile code, identify hotspots, and estimate possible gains.
• Parallelize: Build the simplest working GPU kernel (1 thread per output). Correctness first.
• Optimize: Measure $\rightarrow$ identify limiter $\rightarrow$ change one thing $\rightarrow$ verify speedup.
• Deploy: Validate on production data and ship gains early ($3\times$ shipped beats $30\times$ unfinished).

3. Assess & Amdahl’s Law

Before writing CUDA, evaluate the theoretical ceiling. Amdahl’s Law states that program speedup is strictly limited by its serial (unparallelized) fraction:

• $P$: Parallelizable fraction of execution time.
• $1 - P$: Serial fraction (file I/O, sequential CPU work).
• $N$: Speedup of the parallelized portion.

4. Profiling Cheatsheet & Key Metrics 💻

Never guess—use standard profiling tools to answer specific questions:

Measuring Accurate Time (cudaEvent)

cudaEvent_t start, stop;
cudaEventCreate(&start); cudaEventCreate(&stop);

// Warm-up launch (driver setup overhead)
kernel<<<grid, block>>>(d_in, d_out, n);

cudaEventRecord(start);
kernel<<<grid, block>>>(d_in, d_out, n);
cudaEventRecord(stop);
cudaEventSynchronize(stop);

float ms = 0.0f;
cudaEventElapsedTime(&ms, start, stop);

Application vs Kernel Profiling CLI

# 1. Nsight Systems: Whole-program timeline (CPU, PCIe copies, GPU kernels)
nsys profile --stats=true -o app_report ./my_app
nsys stats app_report.nsys-rep

# 2. Nsight Compute: Single-kernel hardware diagnostics
ncu --section SpeedOfLight --section Occupancy --section WarpStateStats -k regex:myKernel -c 1 ./my_app

Essential Metrics & Formulas

5. The APOD Optimization Log 📝

Track every optimization round with one hypothesis and one change at a time:

Part 2: Parallel Reduction Through APOD 🧩

Now we apply APOD to optimize array summation across multiple rounds.

1. The Reduction Task & General Form

A Reduction combines an array of $n$ elements into a single scalar value using an associative binary operator:

• Sum: $\sum x_i$ (identity: $0$, operator: $+$)
• Product: $\prod x_i$ (identity: $1$, operator: $\times$)
• Min / Max: $\min x_i$ or $\max x_i$ (identity: $\pm\infty$)

2. Assess: Sequential CPU Baseline ($V_0$)

int reduceOnHost(int *in, int n) {
    int s = in[0];
    for (int i = 1; i < n; i++) {
        s += in[i];
    }
    return s;
}

• Complexity: $n - 1$ steps $\rightarrow \mathcal{O}(n)$ time, $n - 1$ additions $\rightarrow \mathcal{O}(n)$ work.
• Dataset: $n = 16,777,217$ ints ($2^{24} + 1 \approx 67.1\text{ MB}$).
• CPU time ($V_0$): $\approx 47.2\text{ ms}$.
• Hardware Speed-of-Light target (Tesla T4 peak 320 GB/s):

Our optimization goal is to push kernel execution toward this $0.21\text{ ms}$ ceiling.

3. Parallelize: The Reduction Tree Concept

Instead of adding numbers sequentially, threads pair up elements in a binary tree:

Figure 2: Parallel reduction tree pairwise summation across time steps.

The Numbers Behind the Tree

• Total Work: $\frac{N}{2} + \frac{N}{4} + \dots + 1 = \mathbf{N - 1\text{ additions}}$ (Work-efficient).
• Execution Steps: $\mathbf{\log_2(N)\text{ steps}}$ (For $N = 1,000,000$, only 20 steps vs 1,000,000 on CPU).
• Peak Resource Demand: $\frac{N}{2} = \mathbf{500,000\text{ threads}}$ at Step 1.
• Average Parallelism:

When is the Tree Allowed?

The parallel tree changes the evaluation order:

• Sequential: $(((a \oplus b) \oplus c) \oplus d)$
• Parallel: $((a \oplus b) \oplus (c \oplus d))$

The operation must satisfy:

• Associative: $(a \oplus b) \oplus c = a \oplus (b \oplus c)$

• Commutative: $a \oplus b = b \oplus a$ (required when reordering operands across threads)

• Allowed: Integer $+$, $\times$, $\min$, $\max$, AND, OR, XOR.

• Disallowed: Subtraction, division (non-associative).

• Note on Floats: Floating-point addition is not strictly associative due to rounding, so GPU and CPU sums can differ slightly.

Parallel Reduction in Real Life

• Sports & Tournaments: Single-elimination brackets to find a tournament winner ($\max$ reduction).
• Big Data & MapReduce: Google and Hadoop partition massive datasets into chunks and use reduction trees to compute global sums, word counts, and aggregate logs.
• Machine Learning & Vision: Batch normalization, loss reduction (MSE, cross-entropy), dot products, vector norms, and image histograms.

4. A Simple Single-Block Reduction Kernel (reduceBlksKernel0)

Before scaling up to millions of elements across multiple blocks, let’s examine the baseline implementation for a single block (Slide 34–35).

Figure 3: Thread-to-data mapping in a single-block reduction kernel across stride doubling steps.

Implementation

__global__ void reduceBlksKernel0(int* in, int* out, int n) {
    int i = 2 * threadIdx.x;

    for (int stride = 1; stride <= blockDim.x; stride *= 2) {
        if (threadIdx.x % stride == 0)
            in[i] += in[i + stride];
        __syncthreads();
    }

    if (threadIdx.x == 0)
        *out = in[0];
}

Step-by-Step Breakdown

1. Thread-to-Data Mapping (i = 2 * threadIdx.x):
   • Each thread initially manages two contiguous elements.
   • Thread 0 starts at index 0, Thread 1 at index 2, Thread 2 at index 4, etc.
2. Iterative Stride Doubling (stride *= 2):
   • At stride = 1: all threads are active (tid % 1 == 0), pairing indices $(0, 1), (2, 3), (4, 5), \dots$.
   • At stride = 2: only even threads participate (tid % 2 == 0), pairing $(0, 2), (4, 6), \dots$.
   • At stride = 4: only threads divisible by 4 participate, pairing $(0, 4), (8, 12), \dots$.
   • After $\log_2(\text{blockDim.x})$ steps, the block’s sum accumulates at in[0].
3. Barrier Synchronization (__syncthreads()):
   • Forces all threads in the block to finish reading and writing their pair before advancing to the next stride.
   • Prevents data races where a fast thread in step $k+1$ reads an index before a slower thread in step $k$ writes to it.
4. Final Output (threadIdx.x == 0):
   • Thread 0 writes the reduced block result in[0] to *out.

Critical Limitations

5. Hierarchical Reduction for Arbitrary Input Lengths 🌐

In real-world applications, arrays contain millions of elements ($n \gg 2048$). Because __syncthreads() cannot synchronize across different blocks, we solve this with Hierarchical Reduction:

1. Partition: Divide the array into independent chunks, assigning one chunk of $2 \times \text{blockSize}$ elements to each block.
2. Local Block Reduction: Each block computes its own partial sum using tree reduction.
3. Global Accumulation: Combine the partial sums into the final result.

Figure 4: Hierarchical reduction partitioning input into blocks to produce partial sums.

Kernel Implementation (V1_naive)

__global__ void reduceV1_naive(int *in, int *partial, int n) {
    // Starting index of this block (each block handles 2 * blockDim.x elements)
    int start = blockIdx.x * 2 * blockDim.x;
    int i = start + 2 * threadIdx.x;

    // Tree reduction inside the block
    for (int stride = 1; stride <= blockDim.x; stride *= 2) {
        if (threadIdx.x % stride == 0) {
            // Guard boundary for arbitrary array sizes
            if (i + stride < n) {
                in[i] += in[i + stride];
            }
        }
        __syncthreads();
    }

    // Thread 0 writes the block's partial sum
    if (threadIdx.x == 0) {
        partial[blockIdx.x] = in[start];
    }
}

Host-Side Launch Configuration

To cover arbitrary array lengths $n$, calculate gridSize using ceiling division:

dim3 blockSize(512);
dim3 gridSize((n + 2 * blockSize.x - 1) / (2 * blockSize.x));

// Launch kernel to compute partial sums per block
reduceV1_naive<<<gridSize, blockSize>>>(d_in, d_partial, n);

// Copy partial sums back to Host and accumulate sequentially on CPU
cudaMemcpy(h_partial, d_partial, gridSize.x * sizeof(int), cudaMemcpyDeviceToHost);
int total_sum = 0;
for (int b = 0; b < gridSize.x; b++) {
    total_sum += h_partial[b];
}

6. Optimize Round 1: Reducing Warp Divergence ($V_1 \rightarrow V_2$) ⚡

Why is $V_1$ Slow? The Warp Divergence Limiter

In $V_1$, the condition if (threadIdx.x % stride == 0) causes severe warp divergence. Consider a block of 256 threads (8 warps):

• stride = 1: All 256 threads are active $\rightarrow$ 0 diverged warps.
• stride = 2: Threads with tid % 2 == 0 run (128 threads). Every warp has 16 active and 16 idle threads $\rightarrow$ All 8 warps diverge!
• stride = 64: Only 4 threads run (tid = 0, 64, 128, 192). They land in Warps 0, 2, 4, 6 $\rightarrow$ 4 warps diverge (each warp runs for only 1 active thread, wasting 31 slots).
• Across all iterations, $V_1$ triggers 47 diverged warp instances. Moreover, the modulo operator (%) is an expensive division instruction on GPU.

The $V_2$ Solution: Pack Active Threads Together

Instead of scattering active threads across all warps, $V_2$ compacts all working threads into the first contiguous threads (threadIdx.x < threshold).

This ensures:

• Active warps run at 100% capacity (no divergence inside).
• Inactive warps have 0 active threads $\rightarrow$ the hardware warp scheduler skips them entirely with zero cost!

V1 (Scattered):  [X . X . X . X .]   -> Warp suffers 50% divergence
V2 (Compacted):  [X X X X . . . .]   -> Warp 0 runs at 100%, Warp 1 disabled

Demystifying the Index Formula

To let adjacent threads process numbers spread across the block, $V_2$ uses:

Why $2 \times \text{stride}$? It acts as a stride multiplier so that lower-index threads jump to the correct data pairs:

• stride = 1: Step is $2 \times 1 = 2$. Thread 0 pairs $(0, 1)$, Thread 1 pairs $(2, 3)$, Thread 2 pairs $(4, 5)$.
• stride = 2: Step is $2 \times 2 = 4$. Thread 0 pairs $(0, 2)$, Thread 1 pairs $(4, 6)$.
• stride = 4: Step is $2 \times 4 = 8$. Thread 0 pairs $(0, 4)$.

Figure 5: V2 thread-to-data mapping packing active threads into the front.

Threads $0, 1, 2, \dots$ do all the work, while warps with higher thread IDs stay completely asleep.

Kernel Implementation (V2_compact)

__global__ void reduceV2_compact(int *in, int *partial, int n) {
    int start = blockIdx.x * 2 * blockDim.x;

    for (int stride = 1; stride <= blockDim.x; stride *= 2) {
        // Stride multiplier lets contiguous threads reach distant elements
        int i = start + threadIdx.x * 2 * stride;

        // Condition packs all active threads into the front
        if (threadIdx.x < blockDim.x / stride) {
            if (i + stride < n) {
                in[i] += in[i + stride];
            }
        }
        __syncthreads();
    }

    if (threadIdx.x == 0) {
        partial[blockIdx.x] = in[start];
    }
}

Profiler Comparison: $V_1$ vs $V_2$

Using Nsight Compute (ncu) to measure active threads per warp instruction:

The $V_3$ Breakthrough: Sequential Addressing (Slide 54)

In $V_2$, active threads were grouped together, but they still read numbers spaced far apart. This caused uncoalesced memory access — the GPU fetched 32 bytes from memory just to use one 4-byte number!

$V_3$ solves this by reversing the reduction tree:

1. Start with a large stride: stride = blockDim.x (half the block’s data).
2. Cut the stride in half each round: stride /= 2 until it reaches 1.
3. Each thread simply accesses its own neighboring index:

Figure 6: Sequential addressing reverses the stride so adjacent threads read adjacent memory words.

Step-by-Step Example (4 Threads, 8 Numbers)

Consider the example from Slide 54 with array [1, 9, 5, 1, 6, 4, 7, 2]:

• Round 1 (stride = 4): Threads 0, 1, 2, 3 are active (threadIdx.x < 4). Each thread pairs with element $i + 4$:

  • Thread 0: in[0] += in[4] $\rightarrow 1 + 6 = \mathbf{7}$
  • Thread 1: in[1] += in[5] $\rightarrow 9 + 4 = \mathbf{13}$
  • Thread 2: in[2] += in[6] $\rightarrow 5 + 7 = \mathbf{12}$
  • Thread 3: in[3] += in[7] $\rightarrow 1 + 2 = \mathbf{3}$
  • Array becomes: [ 7, 13, 12, 3 | 6, 4, 7, 2 ]

• Round 2 (stride = 2): Only Threads 0 and 1 are active (threadIdx.x < 2). Pair with element $i + 2$:

  • Thread 0: in[0] += in[2] $\rightarrow 7 + 12 = \mathbf{19}$
  • Thread 1: in[1] += in[3] $\rightarrow 13 + 3 = \mathbf{16}$
  • Array becomes: [ 19, 16 | 12, 3 | 6, 4, 7, 2 ]

• Round 3 (stride = 1): Only Thread 0 is active (threadIdx.x < 1). Pair with element $i + 1$:

  • Thread 0: in[0] += in[1] $\rightarrow 19 + 16 = \mathbf{35}$
  • Final block sum: in[start] = 35!

Why Is $V_3$ Twice as Fast? Memory Coalescing

• In $V_2$, adjacent threads read numbers spaced far apart $\rightarrow$ wasted bandwidth.
• In $V_3$, Threads 0, 1, 2, 3 read numbers side-by-side in memory $\rightarrow$ the GPU hardware loads all of them in one single 32-byte memory transaction (Coalescing)!

Kernel Implementation (V3_sequential)

__global__ void reduceBlksKernel3(int *in, int *out, int n) {
    int start = blockIdx.x * 2 * blockDim.x;
    int i = start + threadIdx.x;

    // Reverse stride from half-block down to 1
    for (int stride = blockDim.x; stride > 0; stride /= 2) {
        if (threadIdx.x < stride) {
            if (i + stride < n) {
                in[i] += in[i + stride];
            }
        }
        __syncthreads();
    }

    if (threadIdx.x == 0) {
        atomicAdd(out, in[start]);
    }
}

Profiler Evidence: Memory Coalescing in Action

Comparing the hardware memory transactions in Nsight Compute:

ncu --metrics l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum,\
l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum \
-k regex:"reduce.*" ./reduce 256

7. Optimize Round 2: Thread Coarsening & Grid-Stride Loops ($V_4$) ⚙️

Assessing $V_3$: Too Many Threads Doing Too Little

While $V_3$ achieved an impressive $51.5\times$ speedup over CPU, a critical resource inefficiency remains:

• Massive thread count: 1 thread per 2 elements means launching 8 million threads (32,768 blocks) for $n = 16.7\text{M}$.
• Idle threads: After step 1, 50% of threads are idle. After step 3, 87.5% of threads are idle.
• High barrier overhead: Each block executes 9 synchronization barriers (__syncthreads()) just to reduce 512 elements.
• Host burden: The GPU produces 32,768 partial sums that must be transferred and summed.

The Coarsening Concept: Do More Work in Registers

Instead of creating millions of threads that die after one step, Thread Coarsening gives each thread multiple elements to sum locally before entering the tree:

1. Step 1 (Sequential Accumulation): Each thread loops over multiple elements spaced by the grid stride and accumulates the sum in a private register variable.
2. Step 2 (Block Tree Reduction): Once the data is condensed, threads run the tree reduction on just one value per thread.

Figure 7: Thread coarsening aggregates multiple elements sequentially in private registers.

The Grid-Stride Loop Pattern

Rather than launching 32,768 blocks, we launch just enough blocks to saturate all Streaming Multiprocessors (SMs) on the GPU (e.g., $160$ blocks of 256 threads on Tesla T4):

int gid = blockIdx.x * blockDim.x + threadIdx.x;
int totalThreads = gridDim.x * blockDim.x;

int sum = 0;
// Each thread sums elements spaced by the entire grid width
for (int k = gid; k < n; k += totalThreads) {
    sum += in[k];
}

• Scalable: Runs safely for arbitrary array sizes $n$ with any grid configuration.
• Coalesced: In every iteration, thread $k$ and thread $k+1$ access contiguous memory addresses.
• Hardware-friendly: Replaces 32,768 blocks with 160 persistent blocks.

Kernel Implementation (V4_coarsened)

__global__ void reduceV4_coarsened(int *in, int *out, int n) {
    extern __shared__ int sdata[];

    int gid = blockIdx.x * blockDim.x + threadIdx.x;
    int totalThreads = gridDim.x * blockDim.x;

    // 1. Thread coarsening: accumulate multiple elements into a register
    int sum = 0;
    for (int k = gid; k < n; k += totalThreads) {
        sum += in[k];
    }

    sdata[threadIdx.x] = sum;
    __syncthreads();

    // 2. Sequential addressing tree within shared memory
    for (int stride = blockDim.x / 2; stride > 0; stride /= 2) {
        if (threadIdx.x < stride) {
            sdata[threadIdx.x] += sdata[threadIdx.x + stride];
        }
        __syncthreads();
    }

    // 3. Block leader writes result
    if (threadIdx.x == 0) {
        atomicAdd(out, sdata[0]);
    }
}

Benchmark Result: Execution time plunges from $0.85\text{ ms}$ down to $0.325\text{ ms}$! Throughput skyrockets to $206.7\text{ GB/s}$ ($134.8\times$ speedup over CPU).

8. Optimize Round 3: Warp Shuffle & Eliminating Barriers ($V_5$) 🚀

The Final Barrier Bottleneck

In the final 5 steps of the block tree reduction (when stride < 32), all active threads reside inside Warp 0 (threads $0$ through $31$).

Yet, standard shared memory code still executes:

• Multiple reads and writes to shared memory.
• __syncthreads() barrier instructions at every single step!

Because threads in a warp execute synchronously in lockstep, barriers within a single warp are completely redundant.

The Solution: Warp Shuffle (__shfl_down_sync)

CUDA provides warp shuffle instructions that allow threads inside the same warp to read each other’s register values directly across the SM crossbar network:

Figure 8: Intra-warp register-to-register value exchange via __shfl_down_sync.

• Zero Memory Overhead: Values travel directly from register to register without touching shared or global memory.
• Zero Barrier Stalls: No __syncthreads() needed!

inline __device__ int warpReduceSum(int val) {
    for (int offset = 16; offset > 0; offset /= 2) {
        val += __shfl_down_sync(0xffffffff, val, offset);
    }
    return val;
}

Replacing shared memory reduction in Warp 0 with warpReduceSum() shaves off remaining latency, achieving $0.324\text{ ms}$ and $206.8\text{ GB/s}$.

9. Optimize Round 4: Combining Block Results via Atomics ($V_6$) 🎯

How Should We Combine Block Partials?

Once each block finishes reducing its portion, how do we aggregate the final scalar result?

Figure 9: Block leader threads atomically update a single global memory scalar.

Why Atomics Win With Coarsening

When we had 32,768 blocks in $V_1$, having 32,768 blocks write to one address caused heavy memory serialization.

However, with Thread Coarsening ($V_4$), we only launch 160 blocks! Serializing 160 atomic additions takes less than $1\ \mu\text{s}$, completely obliterating the milliseconds spent copying partial arrays to the host:

if (threadIdx.x == 0) {
    atomicAdd(out, sdata[0]); // Only 160 atomic updates across the whole GPU!
}

10. Wrap-up: APOD Journey & The PCIe Reality Check 🏁

Comprehensive Optimization Scorecard (Tesla T4)

Benchmark results for reducing $n = 16,777,217$ integers ($67.11\text{ MB}$ of data) on NVIDIA Tesla T4 ($320\text{ GB/s}$ peak bandwidth):

Approaching the Speed-of-Light (SOL)

With $V_5$, our custom CUDA kernel achieved:

This represents $64.6\%$ of the theoretical physical ceiling ($320\text{ GB/s}$) and comfortably outpaces NVIDIA’s heavily tuned CUB library ($170.4\text{ GB/s}$) for this workload size.

The PCIe Reality Check: Amdahl’s Law Strikes Back!

Before celebrating a $135\times$ speedup, examine the end-to-end wall clock time:

[ CPU Sequential Sum ] ======================================> 43.7 ms

[ GPU Execution Flow ]
|-- PCIe H2D Copy (14.8 ms) --|-- Kernel (0.32 ms) --|-- D2H Copy (0.01 ms) --|
Total GPU Wall-Clock = 15.1 ms

Summary Checklist: The APOD Engineering Takeaway

1. Never guess—profile first: Use nsys for system-level timelines and ncu for warp-level micro-architecture bottlenecks.
2. Eliminate Warp Divergence: Ensure warps are either 100% active or 100% idle.
3. Coalesce Memory Accesses: Ensure adjacent threads access adjacent 4-byte words to utilize all 32 bytes of every memory sector.
4. Coarsen Threads: Don’t launch millions of short-lived threads; use grid-stride loops to accumulate values locally in registers with zero synchronization cost.
5. Leverage Warp-Level Primitives: Replace intra-warp shared memory and __syncthreads() with fast register shuffles (__shfl_down_sync).
6. Aggregate with Atomics: When thread count is kept compact, GPU atomicAdd is far faster than offloading aggregation to the CPU.