Part 5 — Shared Memory

Up to now, we've considered only programs in which each thread could do its work independently of every other thread. This works great for elementwise problems like saxpy, but for example, reduction operations like calculating the sum of an array require intermediate results to be communicated between threads.

Before diving into shared memory as a solution for this problem, let us first think how we could calculate the sum of a very large array with just the tools we have at our disposal so far.

5.1 What we can already do

There is actually a simple implementation that directly reuses our saxpy kernel, if we allow for an in-place calculation and assume that the number of elements is a power of two. Consider the call saxpy(n/2, 1.0, x + n/2, x); This splits the original array into two halves and adds the second half to the first half. After this, we only need to calculate the sum over an array of length n/2. Repeat this process, and after $\log_2 n$ steps, we have calculated the sum of the entire array.

The full function could look like this:

// Sum x[0..n) in place, assuming n is a power of two. Each pass adds the
// upper half onto the lower half, so the running total stays at x[0], and
// after log2(n) passes the whole sum is there.
void reduce_by_halving(int n, float* x) {
    constexpr int BLOCK = 256;
    for (int m = n / 2; m > 0; m /= 2) {
        saxpy<<<(m + BLOCK - 1) / BLOCK, BLOCK>>>(m, 1.0f, x + m, x);
        CUDA_CHECK(cudaGetLastError());
    }
}

There are two crucial ingredients that make an approach like this work: First, we have a common address space, global memory, that can be written to and read from by all blocks. Second, we have an automatic synchronization between kernel launches, which ensures that all the writes finish before any reads start, thus guaranteeing we are reading the correct, up-to-date data. The problem is that these are both expensive: Global memory access can take hundreds of cycles, and the synchronization at the kernel boundary means that the GPU has to wind down all operations, wait for the last straggling threads to finish, wait some more to ensure all memory operations have reached a consistent state, and then start putting threads on the SMs anew. And all this for a single addition worth of real work.

Running the code above on a RTX PRO 4000, it takes about 0.50 ms to sum 33.55 million elements. Every element has to be read from memory at least once, which is 128 MiB of unavoidable traffic, so the reduction is moving its input at 267.70 GB/s. This GPU has a global memory bandwidth of 672 GB/s — only 39.80% of which we are using, so an ideal implementation could be up to 2.50× faster.

5.2 Making the cooperation local

Both costs shrink if we make the operations more local. Instead of synchronizing the full grid of threads, we can do most of the work synchronizing only the threads within a single block. This is advantageous because they all live on the same SM, so no cross-SM communication is needed. As this way we can only communicate intermediate results between threads within one block, a single kernel call can at most reduce the array of n elements to a single number per block. Keeping the previous structure, with 1024 threads per block, we could still do n -> n/2048 instead of n -> n/2 with just one kernel invocation.

The difference is what each barrier makes a warp wait for. At a kernel boundary, every warp in the grid waits for the slowest warp anywhere in the grid, and the device drains and relaunches before the next round begins. A __syncthreads() only makes a warp wait for the other warps of its own block, so a block whose warps finish early moves straight on while a straggler elsewhere is still running:

We could implement a kernel like this:

__global__ void reduce_block(int n, float* x) {
    int bdx = 2 * blockIdx.x * blockDim.x;
    int idx = bdx + threadIdx.x;
    float left = idx < n ? x[idx] : 0.f;
    float right = idx + blockDim.x < n ? x[idx + blockDim.x] : 0.f;
    float sum = left + right;
    if (idx < n) {
        x[idx] = sum;
    }

    int width = blockDim.x / 2;
    while (width > 0) {
        __syncthreads();
        // at this point, the previous results are guaranteed to be available
        if (threadIdx.x < width) {
            float right = idx + width < n ? x[idx + width] : 0.f;
            sum += right;   // left = x[idx] is equal to sum, so no need to reload.
            if (idx < n) {
                x[idx] = sum;
            }
        }
        width /= 2;
    }
}

Then we can run a small secondary kernel that picks up all the partial sums at x[2 * b * blockDim] and produces the final result, which is left as an exercise to the reader.

That second launch is not the only way to finish the job. A block can also add its partial straight into a single location, with an atomic read-modify-write that the hardware guarantees no other thread can interleave with:

atomicAdd(&out[0], sum);

No update is lost, no matter how many blocks arrive at once, and the whole reduction fits into a single kernel launch. The price is that the additions now happen in whatever order the blocks finish in, and floating-point addition is not associative, so the result can differ slightly from one run to the next. Furthermore, as every block only ever adds to out[0], it has to be zeroed beforehand, for example with cudaMemset

Counting that fixup kernel, the whole reduction now takes 0.40 ms, against the 0.50 ms of the launch cascade — the same work, down to the same single number, in 2 kernel launches instead of 25. We are using 49.40% of the available bandwidth where the cascade managed 39.80%.

To summarize __syncthreads() is both a synchronization point and a memory fence: no thread in the block moves past it until every thread in the block has reached it, and every memory write issued before it is visible to every thread after it. It comes with two pitfalls: First, using __syncthreads in non-convergent (i.e., different threads in the block going different ways) branches results in a deadlock. Second, for large blocks, if there is only a single block resident on the SM, then __syncthreads() also causes idle time, because there may not be other warps to be scheduled while waiting for the barrier.

While we have fixed the global synchronization, we are still writing a large amount of intermediate results to global memory. The good news is that now, that same data is almost immediately read again on the same SM, so it is extremely unlikely to do the full round-trip to DRAM, but we still pay bandwidth to L2 cache, which all writes must reach.

We can see the cost when looking at the profiling results. Reading the input is unavoidable — 128 MiB from DRAM, once, the same for any reduction. On top of that, the block-tree version writes every partial through L1 to the L2 cache: 132 MiB of write traffic that all the blocks funnel into the one shared L2. Those partials stay hot enough to be re-read from L1 rather than DRAM — but not all of them: the first round's working set is larger than the cache, so 38 MiB of partials are evicted from L2 (and dropped from L1, which is write-through) and written back to DRAM.

Looking directly at the reported metrics, we get 90.38% of DRAM throughput, which sounds excellent until you remember that about one fourth of the DRAM transfers are superfluous.

It would be great if we could avoid all that memory traffic leaving the SM. And we can! In fact, being able to quickly share data between threads of a block is so useful that there exists dedicated hardware for this: shared memory.

5.3 Shared memory

Shared memory (often abbreviated smem, in contrast to global gmem) is a small scratchpad memory that lives on the SM itself, nowadays carved out of the same on-chip storage as the L1 cache. In contrast to the L1, which is transparent to the programmer, shared memory needs to be explicitly addressed.

Two properties shape how you use it. It is private to the block and lasts exactly as long as the block does — one block cannot see another's shared memory. And it is a scarce resource: because it comes out of the same on-chip budget as the L1 cache, the more a block reserves, the fewer blocks fit on an SM at once, which lowers occupancy.

To use it, we first need to inform the compiler that we require shared memory. In the example above, we could declare

__shared__ float scratch[1024];

which reserves 4096 bytes of shared memory for the block, which will be made available automatically once the block starts executing. Then, instead of writing back our partial sums to global memory, we write them to scratch[threadIdx.x], synchronize the block, and read them back from there — every intermediate now stays on the SM. It shows in the numbers: the same reduction runs in 0.29 ms at 69.60% of peak bandwidth, up from the 49.40% the global-memory block reduction managed.

void reduce(int n, const float* in, float* out);

extern __shared__ float scratch[];
// ...
kernel<<<blocks, threads, threads * sizeof(float)>>>();

extern __shared__ char smem[];
float* values = reinterpret_cast<float*>(smem);
long*  counts = reinterpret_cast<long*>(smem + n * sizeof(float));

cudaFuncSetAttribute(my_kernel, cudaFuncAttributeMaxDynamicSharedMemorySize, bytes);

5.5 Banks and bank conflicts

As we discussed in §3.3, when accessing global memory, it is imperative that the access be coalesced, with all threads in a warp accessing data from the same 128-byte cache line; otherwise the access is split into multiple, serialized requests, slowing down execution. Does the same restriction also apply to shared memory?

Not quite. Shared memory allows for more access patterns to be handled efficiently, but there are still strong constraints if an access is to be handled in a single cycle. To understand these, we need to consider the organization of shared memory.

Shared memory is divided into 32 banks, interleaved at 4-byte granularity, such that bytes 0–3 belong to bank 0, bytes 4–7 to bank 1, and bytes 128–131 are in bank 0 again. Equivalently: the 4-byte word at index i lives in bank i % 32, so consecutive floats land in consecutive banks and the pattern repeats every 128 bytes.

In every cycle, the shared memory controller can access one 4-byte word for each bank, but not multiple words from the same bank. That means that a coalesced access can read 32 floats per cycle, but an access pattern with stride 32 only reads one float every cycle.

float good = smem[tid];
float bad  = smem[32 * tid];

The "bad" situation is called a bank conflict. In this case, as it takes a full 32 cycles to serve that access (excluding communication latency), this conflict is 32-way. Reading smem[2 * tid] would be a 2-way conflict, resolved in two cycles.

So far, these are exactly the same problems as with global memory. But in shared memory, the 32 banks are accessed independently. Therefore, the following access would be conflict-free again:

float free = smem[33 * tid];

While there is a large stride between individual elements, they all live in different banks.

Our reduction is already in the good case: it indexes shared memory by tid, so consecutive lanes address consecutive words, and every step of the tree is conflict-free.

This constraint is much weaker than the one on global memory. There, the 32 addresses of a warp have to fall into as few cache lines as possible, so any stride at all costs us. Here, they only need to land in different banks, and smem[33 * tid] is just as fast as smem[tid]. So shared memory can handle access patterns that global memory cannot, and that is useful for more than passing values between threads: we can also use it to change the layout of our data.

5.6 Matrix transpose

One problem where this turns out to be very valuable is transposing a matrix, out[x][y] = in[y][x]. A direct implementation is another example of an embarrassingly parallel problem:

// One thread per element. Consecutive lanes have consecutive `x`, so the read
// is perfectly coalesced — and the write, indexed by `x` down a column of the
// output, is spread over 32 different cache lines.
__global__ void transpose_naive(int n, const float* in, float* out) {
    int x = blockIdx.x * TILE + threadIdx.x;
    int y = blockIdx.y * TILE + threadIdx.y;
    if (x >= n || y >= n)
        return;

    out[x * n + y] = in[y * n + x];
}

However, there is a catch: As written, consecutive lanes of a warp have consecutive x, so in[y * n + x] reads one contiguous address range. But out[x * n + y] walks down a column: consecutive lanes are n floats apart, and the warp's 32 values land in 32 different cache lines, making writes very slow. And if we attempt to fix this by letting threadIdx.x run down the rows instead, the write becomes coalesced and the read is strided — whichever way round we set it up, one of the two accesses walks a column.

The problem is obvious from profiling data: The kernel produces 524,288 load and 524,288 store requests — one per warp. For the loads, these get translated to 2.1 million 32B sectors, 4 per request for a full cacheline. Stores, however, get serialized to 16.8 million sectors, eight times more than should be required.

The way out is to stop asking global memory to do the transpose at all. A block reads a 32×32 tile of the input with coalesced reads, parks it in shared memory, rereads it from this shared memory, and then writes it out with coalesced writes.

// `tile[threadIdx.x][threadIdx.y]` walks a column: consecutive lanes are TILE
// floats apart, so with TILE == 32 all 32 want the same bank.
__global__ void transpose_tiled(int n, const float* in, float* out) {
    __shared__ float tile[TILE][TILE];

    int x = blockIdx.x * TILE + threadIdx.x;
    int y = blockIdx.y * TILE + threadIdx.y;
    if (x < n && y < n)
        tile[threadIdx.y][threadIdx.x] = in[(size_t)y * n + x];

    __syncthreads();

    // The block's output tile starts at the transposed block coordinate, so
    // this write walks a row of `out` and is coalesced too.
    x = blockIdx.y * TILE + threadIdx.x;
    y = blockIdx.x * TILE + threadIdx.y;
    if (x < n && y < n)
        out[(size_t)y * n + x] = tile[threadIdx.x][threadIdx.y];
}

Of course, the changes above just moved the bad access pattern to a different level in the memory hierarchy instead of eliminating it. This is still useful: Instead of 42.30% bandwidth utilization, we now get 60.10%, a significant improvement. But while store sectors are down to 2.1 million, as expected for fully-coalesced access, now we get 16.8 million wavefronts of shared-memory loads. A shared-memory wavefront — unrelated to AMD's name for a warp in §3.1 — essentially corresponds to a cycle the shared memory spent serving a request. Because we're storing the result of the coalesced global-memory load, store wavefronts are fine: only 0.6 million.

We are hitting the worst-case scenario here: A 32-way shared-memory load conflict, as illustrated in Panel 3 above. That figure also suggests a solution: While a stride of 32 elements leads to maximal conflicts, a stride of 33 avoids them entirely.

Note that no thread in this kernel ever reads a value another thread computed. In the sense of §5.2 there is nothing to cooperate about. We used shared memory purely as a staging area, to turn one bad access pattern into two good ones.

void transpose(int n, const float* in, float* out);

5.8 Recap

Shared memory is a scratchpad that lives on the SM, carved out of the same on-chip storage as the L1 cache, and unlike the L1 it is addressed explicitly.

Both progressions this part built end in the same place:

It belongs to the block. It is allocated when the block starts and is gone when the block ends, and no other block can see it. At the beginning of the block, it has unspecified contents and the user is responsible for initializing it.

It is banked. The 32 banks are interleaved every 4 bytes, and each serves one word per cycle. A warp whose lanes land in 32 different banks is served in one go; lanes wanting different words from the same bank are serialized; lanes wanting the same word get it broadcast for free. What decides this is the stride between consecutive lanes, not how far apart the addresses are in total — which is why smem[33 * tid] costs nothing and smem[32 * tid] costs 32 cycles.

It is scarce. On a recent gaming GPU, there are about 128 KiB of smem, on a datacentre GPU the amount roughly doubles to 228 KiB. All blocks running on that SM must share this budget, thus potentially limiting occupancy (§4.2). Additionally, as smem uses the same physical resources as the L1 cache, using large amounts of shared memory reduces the available L1 cache, potentially leading to more cache misses. In addition to the shared memory used explicitly by the programmer, each block reserves an additional 1 KiB for use of the driver.

Comparison between shared memory and L1 cache In many ways, shared memory and L1 serve similar functions, providing memory that is much closer to the compute cores. They differ in the trade-offs they make: L1 cache is used automatically, so it is the responsibility of the chip to check whether a particular address is available in the cache and where it can be found. This additional tagging stage means that a L1 access is expected to have slightly more latency than shared memory access, where it is the responsibility of the programmer to provide the correct address directly. However, on modern cards, the latency difference for reading is quite small.
In contrast, for writing the difference is stark: As the L1 cache is write-through, everything written to it must also be written to the next-level cache, the L2, making writes much more expensive than to shared memory. This implies that if the main goal is to reuse loaded data across different threads, the benefit of smem might be small, especially considering that L1 contents are shared between threadblocks. But if intermediate results need to be shared, the L1 cache is a bad choice. Finally, as we've seen with the transpose example, shared memory can also be used to turn a bad global-memory pattern into a good one.

Comparison between shared memory and registers Compared to registers, shared memory has a much longer access latency, but is much more flexible. In registers, each thread can only access its own data, up to 255 32-bit values, for about 1 KiB. In contrast, each thread can read and write to the entire shared memory space reserved for its block, potentially much more than 1 KiB. This is although in total, the two are of the same order: the 65k registers of an SM hold roughly 256 KiB, against the 228 KiB of shared memory on a datacentre GPU and the 128 KiB on a gaming one. Secondly, register can only be addressed statically, the register number directly encoded in the instruction, while shared memory supports runtime-calculated addresses. Finally, registers do provide larger bandwidth: Two register banks serve 4 bytes each to every thread, for a total of 256 bytes per cycle, versus shared memory's 128 bytes per cycle.