CUDA Transpose: Coalescing and Bank Padding

Transpose moves addresses, not arithmetic

For a row-major matrix with rows rows and cols columns, element (r,c) starts at r*cols+c. Its destination in the transposed allocation is c*rows+r. Notice that the output stride is the original number of rows. Accidentally retaining cols can pass square tests while corrupting rectangular ones.

A simple baseline lets neighboring threads read neighboring columns. Their output addresses are then separated by rows elements. A shared-memory tile lets one set of threads load rows, synchronize, and read the tile along the other axis before writing contiguous output rows. NVIDIA's transpose tutorial explains this reordering. Lei Mao's CuTe walkthrough expresses related mappings through layout abstractions; the address arithmetic remains worth understanding before adding that abstraction.

Make the boundary concrete

Take a five-by-seven matrix filled with 10*r+c, and use three-by-three tiles in the CPU exercise. The tile beginning at (3,6) contains only (3,6) and (4,6). These become output positions (6,3) and (6,4). A tile loop that assumes nine valid entries fails, even though every complete tile looked correct.

The runnable reference compares a tiled traversal with a direct mathematical transpose and records writes to each output slot. Checking uniqueness catches duplicated writes that a value-only test might miss when inputs contain repeated numbers. This audit is deliberately about index ownership, not GPU timing or thread synchronization.

The following CUDA kernel, not compiled or executed here, uses a 32-by-32 thread block. Its bounds checks accommodate rectangular, nonaligned dimensions; input and output allocations must not overlap.

__global__ void transpose(const float* a, float* b, int rows, int cols) {
    __shared__ float tile[32][33];
    int x = blockIdx.x * 32 + threadIdx.x;
    int y = blockIdx.y * 32 + threadIdx.y;
    if (x < cols && y < rows) tile[threadIdx.y][threadIdx.x] = a[y*cols+x];
    __syncthreads();
    int ox = blockIdx.y * 32 + threadIdx.x;
    int oy = blockIdx.x * 32 + threadIdx.y;
    if (ox < rows && oy < cols) b[oy*rows+ox] = tile[threadIdx.x][threadIdx.y];
}

All threads reach the barrier, including those outside the matrix. The output bounds imply that the corresponding shared-memory cell was populated by a valid input thread. The launch uses ceiling-divided grid dimensions. This large teaching block is a correctness baseline, not a recommendation that 1,024 threads is optimal.

Padding solves a different address problem

In a simplified 32-bank, four-byte-word model, bank index is word_address % 32. Reading a column of a tile with row stride 32 sends distinct words to one bank. Changing the stride to 33 distributes those column accesses across all 32 banks. The demo enumerates both mappings: one unique bank versus 32. This is an address-model result, not a measured conflict count for every architecture or access width.

The optimization ladder should therefore keep three versions: naive, tiled without padding, and the same tile with padding. Otherwise one cannot separate global-access reordering from the bank-layout change. Test squares, skinny matrices, and awkward edges with exact equality, since transpose performs no arithmetic on values.

For a GPU comparison, effective bandwidth counts eight bytes per FP32 element, one read and one write, divided by elapsed time. A copy kernel provides a same-device reference. Full launches, warmup, repeats, and cache conditions must remain visible. The three-version ladder separates global coalescing from bank padding before changes to threads per tile introduce another variable.

rows, cols, tile = 5, 7, 3
source = [10 * r + c for r in range(rows) for c in range(cols)]
actual = [None] * (rows * cols)
writes = [0] * len(actual)

for r0 in range(0, rows, tile):
    for c0 in range(0, cols, tile):
        for r in range(r0, min(r0 + tile, rows)):
            for c in range(c0, min(c0 + tile, cols)):
                destination = c * rows + r
                actual[destination] = source[r * cols + c]
                writes[destination] += 1

expected = [source[r * cols + c]
            for c in range(cols) for r in range(rows)]
assert actual == expected
assert writes == [1] * len(actual)
print(f'transpose {rows}x{cols} -> {cols}x{rows}: correct')
print('last output row:', actual[-rows:])

for stride in [32, 33]:
    banks = [(lane * stride) % 32 for lane in range(32)]
    distinct = len(set(banks))
    print(f'shared row stride={stride}: unique banks={distinct}')
print('Bank-address model only; no hardware counters or timing.')

transpose 5x7 -> 7x5: correct
last output row: [6, 16, 26, 36, 46]
shared row stride=32: unique banks=1
shared row stride=33: unique banks=32
Bank-address model only; no hardware counters or timing.