Puzzle 31: GPU Occupancy Optimization

Why This Puzzle Matters

Building on Puzzle 30: You’ve just learned GPU profiling tools and discovered how memory access patterns can create dramatic performance differences. Now you’re ready for the next level: resource optimization.

The Learning Journey:

• Puzzle 30 taught you to diagnose performance problems using NSight profiling (nsys and ncu)
• Puzzle 31 teaches you to predict and control performance through resource management
• Together, they give you the complete toolkit for GPU optimization

What You’ll Discover: GPU performance isn’t just about algorithmic efficiency - it’s about how your code uses limited hardware resources. Every GPU has finite registers, shared memory, and execution units. Understanding occupancy - the ratio of active warps to maximum possible warps per SM - is crucial for:

• Latency hiding: Keeping the GPU busy while waiting for memory
• Resource allocation: Balancing registers, shared memory, and thread blocks
• Performance prediction: Understanding bottlenecks before they happen
• Optimization strategy: Knowing when to focus on occupancy vs other factors

Why This Matters Beyond GPUs: The principles you learn here apply to any parallel computing system where resources are shared among many execution units - from CPUs with hyperthreading to distributed computing clusters.

Overview

GPU Occupancy is the ratio of active warps to the maximum possible warps per SM. It determines how well your GPU can hide memory latency through warp switching.

This puzzle explores three SAXPY kernels (y[i] = alpha * x[i] + y[i]) with identical math but different resource usage:

alias SIZE = 32 * 1024 * 1024  # 32M elements - larger workload to show occupancy effects
alias THREADS_PER_BLOCK = (1024, 1)
alias BLOCKS_PER_GRID = (SIZE // 1024, 1)
alias dtype = DType.float32
alias layout = Layout.row_major(SIZE)
alias ALPHA = Float32(2.5)  # SAXPY coefficient


fn minimal_kernel[
    layout: Layout
](
    y: LayoutTensor[mut=True, dtype, layout],
    x: LayoutTensor[mut=False, dtype, layout],
    alpha: Float32,
    size: Int,
):
    """Minimal SAXPY kernel - simple and register-light for high occupancy."""
    i = block_dim.x * block_idx.x + thread_idx.x
    if i < size:
        # Direct computation: y[i] = alpha * x[i] + y[i]
        # Uses minimal registers (~8), no shared memory
        y[i] = alpha * x[i] + y[i]

View full file: problems/p31/p31.mojo

fn sophisticated_kernel[
    layout: Layout
](
    y: LayoutTensor[mut=True, dtype, layout],
    x: LayoutTensor[mut=False, dtype, layout],
    alpha: Float32,
    size: Int,
):
    """Sophisticated SAXPY kernel - over-engineered with excessive resource usage.
    """
    # Maximum shared memory allocation (close to 48KB limit)
    shared_cache = tb[dtype]().row_major[1024 * 12]().shared().alloc()  # 48KB

    i = block_dim.x * block_idx.x + thread_idx.x
    local_i = thread_idx.x

    if i < size:
        # REAL computational work that can't be optimized away - affects final result
        base_x = x[i]
        base_y = y[i]

        # Simulate "precision enhancement" - multiple small adjustments that add up
        # Each computation affects the final result so compiler can't eliminate them
        # But artificially increases register pressure
        precision_x1 = base_x * 1.0001
        precision_x2 = precision_x1 * 0.9999
        precision_x3 = precision_x2 * 1.000001
        precision_x4 = precision_x3 * 0.999999

        precision_y1 = base_y * 1.000005
        precision_y2 = precision_y1 * 0.999995
        precision_y3 = precision_y2 * 1.0000001
        precision_y4 = precision_y3 * 0.9999999

        # Multiple alpha computations for "stability" - should equal alpha
        alpha1 = alpha * 1.00001 * 0.99999
        alpha2 = alpha1 * 1.000001 * 0.999999
        alpha3 = alpha2 * 1.0000001 * 0.9999999
        alpha4 = alpha3 * 1.00000001 * 0.99999999

        # Complex polynomial "optimization" - creates register pressure
        x_power2 = precision_x4 * precision_x4
        x_power3 = x_power2 * precision_x4
        x_power4 = x_power3 * precision_x4
        x_power5 = x_power4 * precision_x4
        x_power6 = x_power5 * precision_x4
        x_power7 = x_power6 * precision_x4
        x_power8 = x_power7 * precision_x4

        # "Advanced" mathematical series that contributes tiny amount to result
        series_term1 = x_power2 * 0.0000001  # x^2/10M
        series_term2 = x_power4 * 0.00000001  # x^4/100M
        series_term3 = x_power6 * 0.000000001  # x^6/1B
        series_term4 = x_power8 * 0.0000000001  # x^8/10B
        series_correction = (
            series_term1 - series_term2 + series_term3 - series_term4
        )

        # Over-engineered shared memory usage with multiple caching strategies
        if local_i < 1024:
            shared_cache[local_i] = precision_x4
            shared_cache[local_i + 1024] = precision_y4
            shared_cache[local_i + 2048] = alpha4
            shared_cache[local_i + 3072] = series_correction
        barrier()

        # Load from shared memory for "optimization"
        cached_x = shared_cache[local_i] if local_i < 1024 else precision_x4
        cached_y = (
            shared_cache[local_i + 1024] if local_i < 1024 else precision_y4
        )
        cached_alpha = (
            shared_cache[local_i + 2048] if local_i < 1024 else alpha4
        )
        cached_correction = (
            shared_cache[local_i + 3072] if local_i
            < 1024 else series_correction
        )

        # Final "high precision" computation - all work contributes to result
        high_precision_result = (
            cached_alpha * cached_x + cached_y + cached_correction
        )

        # Over-engineered result with massive resource usage but mathematically ~= alpha*x + y
        y[i] = high_precision_result

View full file: problems/p31/p31.mojo

fn balanced_kernel[
    layout: Layout
](
    y: LayoutTensor[mut=True, dtype, layout],
    x: LayoutTensor[mut=False, dtype, layout],
    alpha: Float32,
    size: Int,
):
    """Balanced SAXPY kernel - efficient optimization with moderate resources.
    """
    # Reasonable shared memory usage for effective caching (16KB)
    shared_cache = (
        tb[dtype]().row_major[1024 * 4]().shared().alloc()
    )  # 16KB total

    i = block_dim.x * block_idx.x + thread_idx.x
    local_i = thread_idx.x

    if i < size:
        # Moderate computational work that contributes to result
        base_x = x[i]
        base_y = y[i]

        # Light precision enhancement - less than sophisticated kernel
        enhanced_x = base_x * 1.00001 * 0.99999
        enhanced_y = base_y * 1.00001 * 0.99999
        stable_alpha = alpha * 1.000001 * 0.999999

        # Moderate computational optimization
        x_squared = enhanced_x * enhanced_x
        optimization_hint = x_squared * 0.000001

        # Efficient shared memory caching - only what we actually need
        if local_i < 1024:
            shared_cache[local_i] = enhanced_x
            shared_cache[local_i + 1024] = enhanced_y
        barrier()

        # Use cached values efficiently
        cached_x = shared_cache[local_i] if local_i < 1024 else enhanced_x
        cached_y = (
            shared_cache[local_i + 1024] if local_i < 1024 else enhanced_y
        )

        # Balanced computation - moderate work, good efficiency
        result = stable_alpha * cached_x + cached_y + optimization_hint

        # Balanced result with moderate resource usage (~15 registers, 16KB shared)
        y[i] = result

View full file: problems/p31/p31.mojo

Your task

Use profiling tools to investigate three kernels and answer analysis questions about occupancy optimization. The kernels compute identical results but use resources very differently - your job is to discover why performance and occupancy behave counterintuitively!

The specific numerical results shown in this puzzle are based on NVIDIA A10G (Ampere 8.6) hardware. Your results will vary depending on your GPU architecture (Pascal, Turing, Ampere, Ada, Hopper, etc.), but the fundamental concepts, methodology, and insights remain universally applicable across all modern NVIDIA GPUs. Use pixi run gpu-specs to get your specific hardware values.

Configuration

Requirements:

• NVIDIA GPU with CUDA toolkit
• NSight Compute from Puzzle 30

⚠️ GPU compatibility note: The default configuration uses aggressive settings that may fail on older or lower-capability GPUs:

alias SIZE = 32 * 1024 * 1024  # 32M elements (~256MB memory per array)
alias THREADS_PER_BLOCK = (1024, 1)  # 1024 threads per block
alias BLOCKS_PER_GRID = (SIZE // 1024, 1)  # 32768 blocks

If you encounter launch failures, reduce these values in problems/p31/p31.mojo:

• For older GPUs (Compute Capability < 3.0): Use THREADS_PER_BLOCK = (512, 1) and SIZE = 16 * 1024 * 1024
• For limited memory GPUs (< 2GB): Use SIZE = 8 * 1024 * 1024 or SIZE = 4 * 1024 * 1024
• For grid dimension limits: The BLOCKS_PER_GRID will automatically adjust with SIZE

Occupancy Formula:

Theoretical Occupancy = min(
    Registers Per SM / (Registers Per Thread × Threads Per Block),
    Shared Memory Per SM / Shared Memory Per Block,
    Max Blocks Per SM
) × Threads Per Block / Max Threads Per SM

The Investigation

Step 1: Test the kernels

pixi shell -e cuda
mojo problems/p31/p31.mojo --all

All three should produce identical results. The mystery: why do they have different performance?

Step 2: Benchmark performance

mojo problems/p31/p31.mojo --benchmark

All three should produce identical results. The mystery: why do they have different performance?

Step 3: Build for profiling

mojo build --debug-level=full problems/p31/p31.mojo -o problems/p31/p31_profiler

Step 4: Profile resource usage

# Profile each kernel's resource usage
ncu --set=@occupancy --section=LaunchStats problems/p31/p31_profiler --minimal
ncu --set=@occupancy --section=LaunchStats problems/p31/p31_profiler --sophisticated
ncu --set=@occupancy --section=LaunchStats problems/p31/p31_profiler --balanced

Record the resource usage for occupancy analysis.

Step 5: Calculate theoretical occupancy

First, identify your GPU architecture and detailed specs:

pixi run gpu-specs

Note: gpu-specs shows all architectural details derived from your GPU hardware - no lookup tables needed!

Common Architecture Specs (Reference):

📚 Official Documentation:

• NVIDIA CUDA Compute Capability Table
• CUDA Programming Guide - Compute Capabilities
• Hopper Architecture In-Depth
• Ampere Architecture Whitepaper

⚠️ Note: These are theoretical maximums. Actual occupancy may be lower due to hardware scheduling constraints, driver overhead, and other factors.

Using your GPU specs and the occupancy formula:

• Threads Per Block: 1024 (from our kernel)

Use the occupancy formula and your hardware specifications to predict each kernel’s theoretical occupancy.

Step 6: Measure actual occupancy

# Measure actual occupancy for each kernel
ncu --metrics=smsp__warps_active.avg.pct_of_peak_sustained_active problems/p31/p31_profiler --minimal
ncu --metrics=smsp__warps_active.avg.pct_of_peak_sustained_active problems/p31/p31_profiler --sophisticated
ncu --metrics=smsp__warps_active.avg.pct_of_peak_sustained_active problems/p31/p31_profiler --balanced

Compare the actual measured occupancy with your theoretical calculations - this is where the mystery reveals itself!

Key Insights

💡 Occupancy Threshold: Once you have sufficient occupancy for latency hiding (~25-50%), additional occupancy provides diminishing returns.

💡 Memory Bound vs Compute Bound: SAXPY is memory-bound. Memory bandwidth often matters more than occupancy for memory-bound kernels.

💡 Resource Efficiency: Modern GPUs can handle moderate register pressure (20-40 registers/thread) without dramatic occupancy loss.

Your task: Answer the following questions

After completing the investigation steps above, answer these analysis questions to solve the occupancy mystery:

Performance Analysis (Step 2):

1. Which kernel is fastest? Which is slowest? Record the timing differences.

Resource Profiling (Step 4):

2. Record for each kernel: Registers Per Thread, Shared Memory Per Block, Warps Per SM

Theoretical Calculations (Step 5):

3. Calculate theoretical occupancy for each kernel using your GPU specs and the occupancy formula. Which should be highest/lowest?

Measured Occupancy (Step 6):

4. How do the measured occupancy values compare to your calculations?

The Occupancy Mystery:

5. Why do all three kernels achieve similar occupancy (~64-66% results may vary depending on gpu architecture) despite dramatically different resource usage?
6. Why is performance nearly identical (<2% difference) when resource usage varies so dramatically (19 vs 40 registers, 0KB vs 49KB shared memory)?
7. What does this reveal about the relationship between theoretical occupancy calculations and real-world GPU behavior?
8. For this SAXPY workload, what is the actual performance bottleneck if it’s not occupancy?

Solution

Register Limit = 65,536 / (19 × 1,024) = 3.36 blocks per SM
Shared Memory Limit = 164KB / 0KB = ∞ blocks per SM
Hardware Block Limit = 32 blocks per SM

Thread Limit = 1,536 / 1,024 = 1 block per SM (floor)
Actual Blocks = min(3, ∞, 1) = 1 block per SM
Theoretical Occupancy = (1 × 1,024) / 1,536 = 66.7%

Register Limit = 65,536 / (25 × 1,024) = 2.56 blocks per SM
Shared Memory Limit = 164KB / 16.4KB = 10 blocks per SM
Hardware Block Limit = 32 blocks per SM

Thread Limit = 1,536 / 1,024 = 1 block per SM (floor)
Actual Blocks = min(2, 10, 1) = 1 block per SM
Theoretical Occupancy = (1 × 1,024) / 1,536 = 66.7%

Register Limit = 65,536 / (40 × 1,024) = 1.64 blocks per SM
Shared Memory Limit = 164KB / 49.2KB = 3.33 blocks per SM
Hardware Block Limit = 32 blocks per SM

Thread Limit = 1,536 / 1,024 = 1 block per SM (floor)
Actual Blocks = min(1, 3, 1) = 1 block per SM
Theoretical Occupancy = (1 × 1,024) / 1,536 = 66.7%

Single Kernel Pass Analysis:
- Input arrays: 32M × 4 bytes × 2 arrays = 256MB read
- Output array: 32M × 4 bytes × 1 array = 128MB write
- Total per kernel: 384MB memory traffic

Peak Bandwidth (A10G): 600 GB/s
Single-pass time: 384MB / 600 GB/s ≈ 0.64ms theoretical minimum
Benchmark time: ~328ms (includes multiple iterations + overhead)

# Find your GPU specs
pixi run gpu-specs

# Measure resource consumption
ncu --set=@occupancy --section=LaunchStats your_kernel

# Measure achieved occupancy
ncu --metrics=smsp__warps_active.avg.pct_of_peak_sustained_active your_kernel

# Always validate with actual performance measurements
ncu --set=@roofline --section=MemoryWorkloadAnalysis your_kernel

OCCUPANCY ANALYSIS → OPTIMIZATION STRATEGY:

High occupancy (>70%) + Good performance:
→ Occupancy is sufficient, focus on other bottlenecks

Low occupancy (<30%) + Poor performance:
→ Increase occupancy through resource optimization

Good occupancy (50-70%) + Poor performance:
→ Look for memory bandwidth, cache, or computational bottlenecks

Low occupancy (<30%) + Good performance:
→ Workload doesn't need high occupancy (memory-bound)