Memoryless Matrix Transposition
I needed a fast and minimal routine that could transpose a warp with 32 elements per lane and store the result to device memory as the final step in a complex kernel.
However, I had a hypothesis that Kepler SHFL instructions and 16-byte stores by pairs of lanes would be functionally equivalent to the standard matrix transpose approach but not require any shared memory or thread block synchronization yet still achieve 100% memory store efficiency.
I banged out the code over the weekend to see what the PTX/SASS would look like and after getting encouraging performance results I used the hypothesized warp-centric “pseudo-transpose” primitive to implement a full matrix transpose.
Results
A 1024x1024 matrix of 32-bit elements transposes on a K20c at ~137 GB/sec. This is ~22% better than the “optimized outer” variant example in the CUDA Toolkit which averages 112 GB/sec.
Here’s the output from the micro benchmark:
Tesla K20c : sm_35 * 13 transposeKernel<<<2048,64>>>(1024 x 1024)
.... Validated! loops (100) : avg 0.05677 ms. = 137.62 GB/sec.
This is pretty good and beats the default NVIDIA example.
But if the “optimized outer” transpose example in the CUDA Toolkit is modified to use a 32x8 block size (32x32 tile size) it jumps from 112 GB/sec to 153.5 GB/sec. on a K20c (758 MHz). Kudos!
If you have a reason to avoid using shared memory then the approach I describe below remains a good option.
Implementation
The excellent Matrix Transpose Example in the CUDA Toolkit uses a square tile of padded shared memory to switch the row-column ordering. 256 elements are loaded from device memory, stored to the shared tile, synchronized, loaded from the shared tile in transposed order and stored back out to device memory. It’s simple and fast. Note that very few registers are actually put to use with this approach. This isn’t a problem since the kernel has a 1:1 thread to element ratio – i.e. many threads are being used.
My transpose routine takes a different approach. Instead of striving to construct wide 128-byte aligned stores, each warp in this routine performs a trivial rearrangement of each lane pair’s elements followed by every lane storing half of an aligned 32-byte transaction.
It’s easier to show how this works with a few illustrations.
We start with a conceptual warp of 8 lanes with 8 32-bit elements per lane:
Next, exchange the odd 16 bytes in even lanes with the even 16 bytes in odd lanes:
Once the shuffle is complete, 32 bytes of elements from a lane column have been transformed into two lanes of 16 bytes in the same register row. The 32-byte groupings are highlighted in different colors:
At this point, each lane pair can perform an efficient 32-byte
st.global.vec4.u32 (or equivalent). The output address calculation
is very simple and dependent on the size of the tile and matrix.
Each lane stores twice in this example. The conceptual 8-lane warp will perform two 128-bit stores per lane. Each 128-bit store results in four 32-byte aligned global store transactions.
Finally, if the transposed warp were reloaded from device memory here is where each color-coded 32-byte store transaction occurred:
Extensions
Pre-Kepler devices require minimal use of shared memory to simulate a SHFL operation but synchronization is still not required as all shared stores and loads are warp-centric.
Transposing 64-bit elements is nearly identical except that only two elements are exchanged per lane pair.
Conclusion
It’s not only possible to transpose a matrix without using shared memory or synchronization but it’s also efficient!