In 2010 Vasily Volkov’s GTC talk “Better Performance at Lower Occupancy” argued that the GPU community is wrong to focus on occupancy as it doesn’t directly relate to throughput. Throughput often increases at lower occupancy because what really keeps a GPU busy is per thread register budget and ILP, not the number of active warps in an SM. This blog re-examines each of his claims on a modern gpu (fp32 cuda cores + registers) and reports what still holds and what doesn’t.

Occupancy


I use RTX 2000 Ada for the experiments done here. Also, since this run is on a RunPod container where GPU performance counters are admin locked, instead of NCU’s achieved occupancy I use theoretical occupancy (from cudaOccupancyMaxActiveBlocksPerMultiprocessor).


Occupancy vs. Throughput

For measuring the occupancy claim I wrote a compute bound fp32 GEMM. The idea is to hold the math per output block constant and vary only the sub tile per thread. A bigger sub tile means fewer threads share the block leading to more registers (more live accumulators per thread) and thus lower occupancy. I fixed the output block size to 64×64 and only let the subtile size vary. I also used __launch_bounds__(threads) so that the register count stays determinant and is driven by the actual registers required, not compiler’s guess (see here). I timed each variant for 50 iterations using cudaEvent timing and plotted the theoretical occupancy (on x-axis) vs achieved GFLOP/s (on y-axis). Check the code in the Appendix.

Occupancy vs Throughput

Variant Threads/block Sub-tile TM×TN Regs/thread Theo. occ GFLOP/s
S 1024 2×2 62 66.7% 1394
A1 512 2×4 40 100.0% 2651
A 256 4×4 64 66.7% 3770
B 128 8×4 96 41.7% 4230
C 64 8×8 150 25.0% 3713
D 32 16×8 254 16.7% 2970
D+ 16 16×16 255 8.3% 411

The conventional model says, more concurrent warps means more latency hiding which ultimately leads to more throughput. It predicts a monotonic curve with 100% occupancy variant (A1) on the top and the 8% occupancy variant (D+) at the bottom. But in practice the peak is at 42% while D+ collapses for another reason (register spill).

Throughput here is decided by per thread ILP and not occupancy. A large sub tile packs more live accumulators into the registers giving each warp independent FMAs to interleave without much warp switching. The proof is that even though S and A are at the same 66.7% occupancy, A is 2.7× faster. Past B the curve falls for a reason Volkov didn’t discuss. D at 254 registers, leaves the compiler no scheduling headroom and D+ pushes past the 255 register cap. This lowers the throughput to 411 GFLOP/s (10× lower than the peak) as the compiler spills to local memory (covered in section 3). No variant reaches the 12 TFLOP/s ceiling but that’s expected as this is a hand written kernel with no tensor cores or double buffering.


Memory bandwidth at low occupancy

For testing “you need lots of threads to hide memory latency”, I wrote a memcpy kernel. By sweeping occupancy from 67% down to about 4% I observed whether the achieved memory bandwidth collapses or stays saturated.

To isolate what drives bandwidth I use thread count and bytes-in-flight. Each variant is a strided copy with __launch_bounds__ and blocks/SM are capped to 2 with 48 KB per block dynamic SMEM allocation. Volkov suggested changing the occupancy by varying the SMEM allocation. For this I vary loads × fetch width per thread. loads is how many independent loads each thread can keep in flight simultaneously and width is the fetch width of each load (float=1, float2=2 and float4=4). So the sweep drops thread count by 16× while raising bytes-in-flight per thread 56× (4 to 224). This is not a tile, section 1’s TM × TN was about per thread compute while loads × width is about each thread’s memory traffic.

I ran each variant on a 256 MB buffer for 50 iterations. Peak memory on this card is 224 GB/s and nvidia-smi dmon (sampled at 1 Hz during the run) reported the mem utilization at 100% for all variants. Check the code in the Appendix.

Bandwidth vs Occupancy

Variant Threads/block Loads × width Bytes/thread Regs/thread Theo. occ GB/s
T1 1024 1×1 4 8 66.7% 194.7
T2 1024 2×1 8 11 66.7% 197.9
T4 512 4×1 16 21 66.7% 206.7
T8 256 8×1 32 38 33.3% 206.7
T8f2 128 8×2 64 38 16.7% 206.7
T8f4 64 8×4 128 40 8.3% 205.7
T14f4 32 14×4 224 50 4.2% 199.3

Note that Regs/thread is the compiler’s total register cost which also includes the two 64 bit base pointers (src and dst) at 2 registers each plus loop index, stride and bound. Normally people try to hide the load latency by having multiple threads loading simultaneously but as we can see from the experiment number of threads doesn’t matter that much. One thread issuing 14 float4 loads keeps as much in flight as 56 threads each issuing one float. So given the loads are independent per thread we don’t need to maximize occupancy.


Register pressure: the modern ceiling

Section 1 lowered the occupancy by increasing the tile size, bigger tile led to more accumulators leading to more registers. Here I do the opposite and take one kernel (section 1’s variant C) and cap its register count with -maxrregcount. I recompiled variant C eight times with different -maxrregcount values (see the table below) to obtain the different binaries. Each binary’s ptx log gave the actual register count, spill stores, spill loads and stack frame (no profiler was used in this). ( You can check the code in the Appendix)

Register Pressure vs Occupancy

Cap (-maxrregcount) Actual regs Theo. occ ms/call GFLOP/s Spill stores Spill loads Stack frame
32 32 41.7% 547.44 251 4840 4652 1672
48 48 41.7% 331.78 414 3160 2964 904
64 64 41.7% 223.81 614 2140 1996 696
96 96 41.7% 41.15 3340 200 140 80
128 128 33.3% 32.47 4233 0 0 0
160 160 25.0% 36.79 3736 0 0 0
192 154 25.0% 36.92 3723 0 0 0
255 154 25.0% 36.86 3729 0 0 0

Below 128 regs the kernel spills to the local memory and throughput increases from 251 up to 3340 GFLOP/s as the spill traffic decreases. We observe the first zero spill at 128 regs and max throughput at the same point of 4233 GFLOP/s. Above 128 more registers simply hurt the throughput as it decreases the occupancy (SM can only fit so many threads if each require large registers) without adding to ILP. One thing to note is that the compiler stops allocating at 154 actual regs. So Volkov’s argument to use more registers still holds with one caveat that we should spend registers to get more ILP but never past the spill cliff. The same variant C that in section 1 gave 3713 GFLOPS/s hits 4233 at 128 capped registers.

What survived

All three sections give the same result / intuition that occupancy shouldn’t be a metric we should be focusing on. For compute bound tasks we should be optimizing per thread ILP (section 1), bounded above by the spill cliff (section 3). For memory bound it’s bytes-in-flight and not the thread count (section 2). Volkov was right in 2010 and he’s still right on modern fp32 hardware with couple of minor exceptions that we should be considering while designing kernels.

In this blog all three sections ran on fp32 CUDA cores but a modern GPU does most of its throughput work on tensor cores, so in Part II I will run tests tensor cores and report how they affect register allocation and whether the modern warp scheduler has absorbed the manual ILP tuning.

Appendix

Below are the kernel-only code. The host side harness (launch setup, timing, the template param / -maxrregcount sweep) is elided for brevity.

Occupancy vs Througput code

template <int BM, int BN, int BK, int TM, int TN>                   // template for different shapes
__global__ __launch_bounds__((BM * BN) / (TM * TN))                 // (64*64)/(4*4) - threads per block
void gemm_kernel(const float* __restrict__ A,
                 const float* __restrict__ B,
                 float* __restrict__ C,
                 int M, int N, int K, float alpha, float beta) {
    constexpr int THREADS = (BM * BN) / (TM * TN);
    constexpr int COLS = BN / TN;

    const int bx = blockIdx.x;  // along N
    const int by = blockIdx.y;  // along M
    const int tx = threadIdx.x;
    const int tr = tx / COLS;
    const int tc = tx % COLS;

    __shared__ float As[BM][BK + 1];   // +1 pad to reduce bank conflicts on BK stride
    __shared__ float Bs[BK][BN + 1];
    
    // set the accumulator
    float acc[TM][TN];
    #pragma unroll
    for (int m = 0; m < TM; m++)
        #pragma unroll
        for (int n = 0; n < TN; n++) acc[m][n] = 0.0f;
    
    // tiled matrix mul
    #pragma unroll 1
    for (int kk = 0; kk < K; kk += BK) {
        #pragma unroll
        for (int i = tx; i < BM * BK; i += THREADS) {
            const int r = i / BK;
            const int c = i % BK;
            As[r][c] = A[(by * BM + r) * K + (kk + c)];
        }
        #pragma unroll
        for (int i = tx; i < BK * BN; i += THREADS) {
            const int r = i / BN;
            const int c = i % BN;
            Bs[r][c] = B[(kk + r) * N + (bx * BN + c)];
        }
        __syncthreads();

        #pragma unroll
        for (int k = 0; k < BK; k++) {
            float ar[TM];
            #pragma unroll
            for (int m = 0; m < TM; m++) ar[m] = As[tr * TM + m][k];
            float br[TN];
            #pragma unroll
            for (int n = 0; n < TN; n++) br[n] = Bs[k][tc * TN + n];

            #pragma unroll
            for (int m = 0; m < TM; m++)
                #pragma unroll
                for (int n = 0; n < TN; n++)
                    acc[m][n] = ar[m] * br[n] + acc[m][n];
        }
        __syncthreads();
    }

    #pragma unroll
    for (int m = 0; m < TM; m++)
        #pragma unroll
        for (int n = 0; n < TN; n++) {
            const int row = by * BM + tr * TM + m;
            const int col = bx * BN + tc * TN + n;
            if (row < M && col < N)
                C[row * N + col] = alpha * acc[m][n] + beta * C[row * N + col];
        }
}

Bandwidth vs Occupancy

template <int THREADS, int V, int WT>
__global__ void __launch_bounds__(THREADS)
memcpy_kernel(const float* __restrict__ src, float* __restrict__ dst, size_t n) {
    extern __shared__ char occ_smem[];
    if (threadIdx.x == 0) occ_smem[0] = (char)0;

    // strided grid: thread t (global index g) writes contiguous floats,
    // but iterates over V tiles separated by gridStride.  This issues V
    // independent loads per thread so the memory subsystem can keep them
    // all in flight (Volkov's "no latency stall until dependent op").
    const size_t grid_stride = (size_t)gridDim.x * THREADS;
    size_t tid = (size_t)blockIdx.x * THREADS + threadIdx.x;
    for (int i = 0; i < V; ++i) {
        size_t base = tid + (size_t)i * grid_stride;
        if (WT == 4) {
            const float4* s4 = reinterpret_cast<const float4*>(src);
            float4*       d4 = reinterpret_cast<float4*>(dst);
            size_t idx = base;  // floats already in float4 units when WT==4
            if (idx * 4 + 3 < n)
                d4[idx] = s4[idx];
        } else if (WT == 2) {
            const float2* s2 = reinterpret_cast<const float2*>(src);
            float2*       d2 = reinterpret_cast<float2*>(dst);
            size_t idx = base;
            if (idx * 2 + 1 < n)
                d2[idx] = s2[idx];
        } else {
            if (base < n) dst[base] = src[base];
        }
    }
}

Spill vs Occupancy

Same as the first kernel in Occupancy vs Throughput