Part 4 — Pipelining and Occupancy

In this part, we will consider how and why a GPU executes tens of thousands of threads efficiently.

4.1 Pipelining and latency hiding

In the introduction, we claimed that the speed of GPUs was due to them executing tens of thousands of threads concurrently. With a high-end GPU having about 100 SMs and each SM coming with 128 "Cuda Cores"¹, we would end up with roughly 12800 parallel threads. While this is an enormous number, it still falls short of the 153600 threads such a GPU would be able to support. Where do these additional threads come from, and why would we want them?

Consider a calculation like d = a + b + c. First, you need to add a and b, then you add the result to c. The second addition cannot begin before the first one has finished. As most instructions, including addition, take several cycles to complete, the warp would be stalled for a time, not doing any useful work. But the core is busy doing the addition, no? Yes, but most instructions are pipelined, in the sense that for an instruction that works in 4 phases, i.e., in 4 cycles, there would be one piece of hardware for every phase. In the first cycle, phase 1 is active, and its result moves to phase 2 in the next cycle, leaving it ready to start phase 1 for a different input immediately.

So while we have to wait four cycles to get our result (latency), we could be getting one result every cycle (throughput). Except that, in a single thread, instructions frequently depend on one another, and we spend most of the time waiting for results to be available. A CPU fills those gaps by looking down the instruction stream of the one thread it is running for something that does not depend on the pending result, and executing that out of order. This requires complicated machinery to determine which instructions can be executed safely at what time. An alternative would be to use the bubble in one thread to run instructions from another thread.

That is exactly what the GPU does: It keeps more warps resident than it can execute in parallel, and if one warp stalls it will simply switch to another one that is ready to proceed. This strategy can only work when the switch is essentially free.

Why switching warps is free

To understand how this is possible, let's zoom into the architecture of a single SM: Within the SM, warps are distributed over four sub-partitions, each with its own warp scheduler, its own slice of the register file, and its own execution units (aka Cuda cores):

A warp is assigned to one sub-partition when its block starts, and stays there until it exits, unlike on a CPU where threads can migrate between cores. In particular, each warp gets its own slice of the register file; each sub-partition contains 16,384 32-bit registers, which would be 512 registers per thread in the warp, even though at most 255 logical registers can be addressed. And with most programs not coming even close to that number, the sub-partition has enough registers to hold the state of multiple warps at the same time.

With the full state of multiple warps directly available, issuing instructions from different warps back to back does not require any context switch, no saving and restoring of state from far-away memory. In this way, the warp scheduler can maintain a list of warps and keep track of their status, picking one of those that have all their dependencies fulfilled every cycle.

Unfortunately, this method of solving the latency problem introduces a new trade-off on a shared resource: You can have many threads, each using only a few registers, or few threads with a large register footprint.

4.2 Occupancy

Each warp scheduler has resources to track at most 12 (on consumer) or 16 (on datacentre) warps at the same time. If a kernel's blocks can fill all of those slots, we say it runs at 100% occupancy.

The figure below illustrates how occupancy of 1/12, 4/12, and 8/12 warps can hide latency:

The arithmetic in the figure generalises. A warp whose next instruction depends on one that takes $L$ cycles to retire can issue once every $L + 1$ cycles, so it takes $L + 1$ such warps to give a scheduler something to do every cycle. For arithmetic that is an achievable number: a dependent chain of multiply-adds has $L$ of about four.

Occupancy is bounded by whichever of the SM's budgets runs out first:

Every one of these is a per-SM budget divided by a per-block or per-thread demand, so the arithmetic is always the same. The planar greyscale kernel of §3.3 compiles to 16 registers per thread, so a full 1536 threads need 24576 of the 65536 registers — not binding. It declares no shared memory. At 256 threads per block, filling the SM takes six blocks, well inside the limit of 24. Nothing runs out, and its theoretical occupancy is 100%.

On the other hand, a kernel using 64 registers per thread tops out at 65536 / 64 = 1024 threads, so 32 of the 48 warp slots can be filled, and occupancy is capped at 67%. To see how many registers a kernel is using, we can enable verbose output: nvcc --ptxas-options=-v reports the register count at compile time:

ptxas info    : Compiling entry function '_Z16grayscale_planariPKhS0_S0_Ph' for 'sm_120'
ptxas info    : Function properties for _Z16grayscale_planariPKhS0_S0_Ph
    0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
ptxas info    : Used 16 registers, used 0 barriers

Note also that this reports 0 spills. Register spills are when the compiler is unable to fit all variables into registers, either because it is running out or because it cannot resolve which data is accessed at compile time.

In many cases, the chosen block size also influences the achievable occupancy. Consider this example: Assuming that 24 warps would fit on the SM, then using a block size of 16 warps would allow only one block, leaving one third of the potential occupancy unused. If instead a block size of 12 warps is used, two blocks can run in parallel, and we achieve the maximum occupancy given the resource constraints.

The function cudaOccupancyMaxPotentialBlockSize returns a block size reaching the best occupancy that kernel can achieve, along with the smallest grid that fills the device with it:

int minGridSize, blockSize;
cudaOccupancyMaxPotentialBlockSize(&minGridSize, &blockSize, grayscale_planar, 0, 0);
// on this card: blockSize = 768, minGridSize = 140, which is 2 blocks on each of 70 SMs

To do the arithmetic yourself rather than take a recommendation, ask how many blocks fit:

int blocks;
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&blocks, grayscale_planar, 256, 0);  // 6

Note that this is theoretical occupancy. What a kernel actually sustains is lower, and is a property of the whole launch: the SM is not yet full while the first blocks are starting, and empties again as the last wave (§2.3) drains. The planar greyscale kernel sits at 77.20% achieved occupancy against a theoretical 100%. Further, if the input data is too small to launch sufficiently many blocks, achieved occupancy will also be below the theoretical limit.

4.3 Register reuse

The most immediate resource that different warps running on the same sub-partition contend on is the register file. While the use of registers is automatically decided by the compiler, in many cases we can change our algorithm such that it uses more registers in order to decrease the latency of individual instructions.

One prominent example is algorithms which exhibit data reuse. Consider the task of multiplying two matrices A and B. Written from the perspective of one output C_ij, it looks like a perfect GPU problem: every output can be calculated independently from an inner product of one row of A and one column of B. This gives rise to the triple-nested for-loop version of matrix multiplication:

for (int i = 0; i < M; ++i) {
    for (int j = 0; j < N; ++j) {
        C[i][j] = 0;
        for (int k = 0; k < K; ++k) {
            C[i][j] += A[i][k] * B[k][j];
        }
    }
}

There are several things wrong with this code. First, each addition operation C[i][j] += A[i][k] * B[k][j]; requires loading and storing C[i][j], as it resides in global memory and these updates are, at least potentially, visible to other threads. As memory instructions are very expensive, this code will be very slow. This problem can be avoided by using a local variable. As it is private to the thread, the compiler can safely decide to put its value into a register, making reads and writes instantaneous.

This results in the following implementation

for (int i = 0; i < M; ++i) {
    for (int j = 0; j < N; ++j) {
        T sum = 0;
        for (int k = 0; k < K; ++k) {
            sum += A[i][k] * B[k][j];
        }
        C[i][j] = sum;
    }
}

How much work is this doing? We have three loops, over M, N, K, and each iteration performs a single multiplication and addition (these turn into a one fused multiply-add in hardware). Thus, we need O(M*N*K) arithmetic operations. But even after fixing the local accumulation, each loop iteration also reads two inputs, so in total 2 * M * N * K reads, despite only (M+N) * K different inputs. Each input value is read many times over!

There is a second, equivalent view on matrix multiplication, that considers it to be the sum over many outer products. In this case, the outer loop is over K, yielding

for (int k = 0; k < K; ++k) {
    for (int i = 0; i < M; ++i) {
        T a_ik = A[i][k];
        for (int j = 0; j < N; ++j) {
            C[i][j] += a_ik * B[k][j];
        }
    }
}

Each value of A need only be loaded once, but we are back in the situation where we waste memory traffic on C.

When considering the whole matrix as a unit, we cannot produce memory-efficient code. The key insight, then, is to chunk up the computation into tiles: within each tile we use a load-efficient outer-product approach, but across tiles we handle the computation like in the inner-product example so that the accumulator can be kept in registers:

constexpr int TM = 4;
constexpr int TN = 4;

for (int io = 0; io < M; io += TM) {
    for (int jo = 0; jo < N; jo  += TN) {
        T sum[TM][TN] = {};      // accumulator
        for (int k = 0; k < K; ++k) {
            T at[TM];
            T bt[TN];

            // input fetching
            for (int ii = 0; ii < TM; ++ii) {
                at[ii] = A[io + ii][k];
            }
            for (int ji = 0; ji < TN; ++ji) {
                bt[ji] = B[k][jo + ji];
            }

            // outer-product accumulation
            for (int ii = 0; ii < TM; ++ii) {
                for (int ji = 0; ji < TN; ++ji) {
                    sum[ii][ji] += at[ii] * bt[ji];
                }
            }
        }

        // output writing -- epilogue
        for (int ii = 0; ii < TM; ++ii) {
            for (int ji = 0; ji < TN; ++ji) {
                C[io + ii][jo + ji] = sum[ii][ji];
            }
        }
    }
}

Now, at, bt and sum are all local variables, and will be placed in registers.

Implementing this (for T = float) on the GPU and sweeping over different tile sizes results in the following speeds for M = N = K = 4096:

As expected, increasing the tile size decreases the amount of DRAM read, as more data is reused directly from register. But at the same time, the number of registers increases with TN x TM, so once we push up against the register limit, occupancy falls. Nonetheless, having fewer global memory loads means there are fewer instructions that need their latency hidden, so the fact that we have fewer resident warps does not offset the massive speed increase due to memory efficiency. That is, until we reach 16x8. At that size, the benefits of fewer loads diminish, and the reduced occupancy starts taking its toll.

Growing the tile even further, we simply run out of registers. The compiler happily accepts that code, but it now generates additional store and (re)load instructions that spill data that does not fit in registers into local memory. Not only does this add wasteful instructions, but the spill traffic crowds A and B out of the L1, so the DRAM traffic the larger tile had saved comes back several times over.

void matmul(int n, const uint32_t* a, const uint32_t* b, uint32_t* c);

for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j) {
        uint32_t sum = 0;
        for (int k = 0; k < n; ++k)
            sum += a[i * n + k] * b[k * n + j];
        c[i * n + j] = sum;
    }