Coalescing: Shared-Memory Matrix Transpose

This kernel transposes an N x N matrix correctly, but it is slow for a reason that has nothing to do with the math: its writes to global memory are uncoalesced. Neighbouring threads write to rows N apart, so each warp's store turns into many separate memory transactions. Stage a tile in shared memory so both the read and the write hit contiguous addresses.

Languages: CUDA C++, Python, Rust

Lesson

The GPU moves global memory in wide bursts. When the 32 threads of a warp touch 32 consecutive addresses, the hardware serves them in one coalesced transaction; when they scatter across far-apart addresses, it issues many transactions and most of each burst is wasted. A naive transpose reads A along rows (coalesced) but writes C along columns (strided by N, uncoalesced), so the write side runs at a fraction of memory bandwidth. The classic fix is a shared-memory tile: the block reads a contiguous tile of A into shared memory coalesced, then writes it back out to C coalesced, doing the actual index swap inside fast shared memory where scattered access is cheap.

Shared memory changes the global-store mapping, but column access can cause bank conflicts. Continue to Matrix Transpose: Pad Away Shared-Memory Bank Conflicts. Logical tile width and physical shared-memory stride are different.

What you will learn

Recognize coalesced vs uncoalesced global memory access and why the naive transpose's writes are uncoalesced.
Stage a TILE x TILE block in shared memory so both the global read and the global write are coalesced.
Do the row/column index swap inside shared memory instead of against global memory.
Handle partial tiles at the matrix edge when N is not a multiple of the tile size.

Starter code

    // Naive transpose. Reading A is coalesced (consecutive threads read consecutive
    // columns of one row), but writing C is strided by N - each thread in a warp writes
    // to a different row, so the stores scatter into many memory transactions.
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x < N && y < N) {
        C[x * N + y] = A[y * N + x];
    }