Tiled Matrix Transposition with Shared Memory

The Inefficiency of a Naive Transpose

Let's first consider the most straightforward way to write a matrix transpose. An out-of-place transpose reads an element from idata[row][col] and writes it to odata[col][row]. In a CUDA kernel where each thread handles one element, the logic looks simple.

However, as you know from our discussion on the memory hierarchy, the pattern of memory access is everything. For matrices stored in row-major order (the standard in C/C++), contiguous threads in a warp accessing contiguous elements in a row results in a single, efficient, coalesced memory transaction. Accessing elements in a column means each thread hits a memory location far from its neighbor, resulting in many slow, uncoalesced transactions.

A naive transpose kernel does one of each:

1. Reads from idata: Threads with consecutive threadIdx.x values will read from consecutive columns in the same row. This is a coalesced and efficient access pattern.
2. Writes to odata: These same threads then write to consecutive rows in the same column. This is an uncoalesced and highly inefficient access pattern, as each write goes to a completely different memory segment.

This uncoalesced write becomes a major performance bottleneck. Let's examine a resource that quantifies this exact problem.

The performance table you just saw makes the issue undeniable. The transposeNaive kernel achieves only a small fraction of the available memory bandwidth precisely because of the strided, uncoalesced global memory writes.

The Tiled Solution: Using Shared Memory to Reorder Data

To solve this, we cannot change how global memory works. Instead, we must change our access pattern. The solution is to use shared memory as a small, fast, on-chip scratchpad to reorder the data before writing it back to global memory.

This approach is called tiling. The algorithm for a thread block is as follows:

1. Load (Coalesced): Cooperatively load a 2D tile of the input matrix idata into a 2D shared memory array. The load is performed row by row, ensuring all global memory reads are coalesced.
2. Synchronize: Call __syncthreads() to ensure every thread in the block has finished its part of the load. This prevents any thread from reading from the shared memory tile before it's fully populated.
3. Store (Coalesced): Read the data from the shared memory array, but with transposed indices. Write this data to the output matrix odata. Because the threads are writing to different global memory tiles, this access pattern is also coalesced.

This process is visualized perfectly in the following diagram. A block of threads reads a tile row-wise from idata, places it in the shared memory tile, and after synchronization, reads it out column-wise to write it row-wise into the correct location in odata.

Now, let's look at how this is implemented in code.

Deconstructing the Tiled Transpose Kernel

Let's break down the transposeCoalesced kernel from the article. For simplicity, let's assume a TILE_DIM of 32 and a BLOCK_ROWS of 8, meaning our block has 32x8 = 256 threads, and it processes a 32x32 tile.

__global__ void transposeCoalesced(float *odata, float *idata, int width, int height) {
    // 1. Declare a shared memory tile
    __shared__ float tile[TILE_DIM][TILE_DIM];

    // 2. Calculate global indices for reading from idata
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x;
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y;
    
    // Each thread loads TILE_DIM / BLOCK_ROWS elements (e.g., 32/8 = 4)
    for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS) {
        // Read is coalesced: contiguous threads read contiguous memory
        tile[threadIdx.y + j][threadIdx.x] = idata[(yIndex + j) * width + xIndex];
    }

    // 3. Synchronize to make sure the whole tile is loaded
    __syncthreads();

    // 4. Calculate global indices for writing to odata (transposed block)
    xIndex = blockIdx.y * TILE_DIM + threadIdx.x;
    yIndex = blockIdx.x * TILE_DIM + threadIdx.y;

    for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS) {
        // Write is now coalesced: contiguous threads write to contiguous memory
        odata[(yIndex + j) * height + xIndex] = tile[threadIdx.x][threadIdx.y + j];
    }
}

Key Implementation Points:

• Shared Memory Declaration: __shared__ float tile[TILE_DIM][TILE_DIM]; declares a 2D array in shared memory, private to each thread block.
• Coalesced Load: In the first loop, idata[(yIndex + j) * width + xIndex] is read. For threads in a warp, threadIdx.x varies, so they read from consecutive memory addresses xIndex, xIndex+1, etc. This is a perfectly coalesced read. The data is written to tile[threadIdx.y + j][threadIdx.x], effectively mapping the global tile directly to the shared tile.
• Synchronization: The __syncthreads() is the crucial barrier that ensures no thread proceeds to the writing phase until all threads have completed the reading phase.
• Transposed Read from Shared Memory: In the second loop, the magic happens. The data is read from tile[threadIdx.x][threadIdx.y + j]. Notice the indices are swapped: threadIdx.x is now the row index and threadIdx.y is the column index. This performs the transpose on-chip.
• Coalesced Store: The data is written to odata[(yIndex + j) * height + xIndex]. The xIndex for the write is calculated based on threadIdx.x, so once again, contiguous threads write to contiguous global memory locations, achieving a coalesced write.

By using shared memory as an intermediary, we have transformed the uncoalesced global write into a series of fast, on-chip shared memory accesses followed by a fully coalesced global write. The performance improvement shown in the article's table is a direct result of this transformation.

Conclusion

In this lesson, you have implemented one of the most fundamental optimization patterns in CUDA programming. You saw firsthand how a seemingly benign operation like matrix transposition can expose severe performance bottlenecks related to memory access patterns. By applying your knowledge of shared memory and synchronization, you implemented a tiled algorithm that rectifies this issue.

Key Takeaways:

• A naive matrix transpose on row-major data suffers from uncoalesced global memory writes, which severely limits performance.
• Tiling is a technique that uses shared memory as a fast, on-chip scratchpad to stage and reorder data.
• The tiled transpose algorithm works by first performing a coalesced read from global memory into a shared memory tile.
• After synchronizing with __syncthreads(), threads read from the shared memory tile with transposed indices and perform a coalesced write to global memory.
• This pattern effectively converts slow, strided global memory writes into fast, arbitrary-pattern shared memory reads, restoring performance.

While our transposeCoalesced kernel is a massive improvement, it's still not perfect. If you look closely at the performance table in the NVIDIA blog post, you'll see there's still a gap between our tiled transpose and a simple memory copy. This remaining gap is caused by a more subtle issue within shared memory itself. In our next lesson, we will investigate this phenomenon, known as shared memory bank conflicts, and learn how to resolve them to unlock even more performance.