Computer algorithms are some what extra friendly towards data sizes that are powers of two. Especially GPU compute algorithms are very favorable to data sizes that are multiples of 32(which is also a power of two, known to the GPUs by the name Warp/Wavefront). However, in most real world situations it is rather hard to find the data to be so conveniently sized. In today's post, we shall be looking at one of such scenarios related to GPU compute. To be specific, we shall provide you with some tips on how to optimize matrix transpose algorithm for a GPU.

Let us start with the transpose kernel available from <a href="http://devblogs.nvidia.com/parallelforall/efficient-matrix-transpose-cuda-cc/">NVIDIA's Parallel Forall blog</a>. It has been optimized to avoid bank conflicts as well but only works on matrices with dimensions that are multiples of 32.
<pre>
template<typename T>
__global__ 
void transpose32(T * out, const T * in, unsigned dim0, unsigned dim1)
{
    __shared__ T shrdMem[TILE_DIM][TILE_DIM+1];

    unsigned lx = threadIdx.x;
    unsigned ly = threadIdx.y;

    unsigned gx = lx + blockDim.x * blockIdx.x;
    unsigned gy = ly + TILE_DIM   * blockIdx.y;

#pragma unroll
    for (unsigned repeat = 0; repeat < TILE_DIM; repeat += blockDim.y) {
        unsigned gy_ = gy+repeat;
        shrdMem[ly + repeat][lx] = in[gy_ * dim0 + gx];
    }
    __syncthreads();

    gx = lx + blockDim.x * blockIdx.y;
    gy = ly + TILE_DIM   * blockIdx.x;

#pragma unroll
    for (unsigned repeat = 0; repeat < TILE_DIM; repeat += blockDim.y) {
        unsigned gy_ = gy+repeat;
        out[gy_ * dim0 + gx] = shrdMem[lx][ly + repeat];
    }
}
</pre>
A minor modification has been done to the kernel from NVIDIA website by templating the matrix element data type such that we can use the same kernel with other integral data types.

Now, we make the kernel applicable to matrices with any size by adding index checks. Following are the changes required to do so.
<pre>
template<typename T>
__global__
void transposeNoSC(T * out, const T * in, unsigned dim0, unsigned dim1)
{
    /*
     * same code as transpose32 kernel
     */
#pragma unroll
    for (unsigned repeat = 0; repeat < TILE_DIM; repeat += blockDim.y) {
        unsigned gy_ = gy+repeat;
        if (gx<dim0 && gy_<dim1)
            shrdMem[ly + repeat][lx] = in[gy_ * dim0 + gx];
    }
    __syncthreads();
    /*
     * same code as transpose32 kernel
     */
#pragma unroll
    for (unsigned repeat = 0; repeat < TILE_DIM; repeat += blockDim.y) {
        unsigned gy_ = gy+repeat;
        if (gx<dim1 && gy_<dim0)
            out[gy_ * dim0 + gx] = shrdMem[lx][ly + repeat];
    }
}
</pre>

There is a reason the above kernel is named transposeNoSC which shall become evident as we proceed to next modification. Lets add back the modifications to handle the special case where matrix dimensions are multiples of 32.
<pre>
template<typename T, bool is32Multiple>
__global__
void transposeSC(T * out, const T * in, unsigned dim0, unsigned dim1)
{
    /*
     * same code as transpose32 kernel
     */
#pragma unroll
    for (unsigned repeat = 0; repeat < TILE_DIM; repeat += blockDim.y) {
        unsigned gy_ = gy+repeat;
        if (is32Multiple || (gx<dim0 && gy_<dim1))
            shrdMem[ly + repeat][lx] = in[gy_ * dim0 + gx];
    }
    __syncthreads();
    /*
     * same code as transpose32 kernel
     */
#pragma unroll
    for (unsigned repeat = 0; repeat < TILE_DIM; repeat += blockDim.y) {
        unsigned gy_ = gy+repeat;
        if (is32Multiple || (gx<dim1 && gy_<dim0))
            out[gy_ * dim0 + gx] = shrdMem[lx][ly + repeat];
    }
}
</pre>
Notice the additional boolean template parameter <code>is32Multiple</code> that is used to short circuit the <code>if</code> condition within the for loops inside the kernel. Where ever the transposeSC is called upon with boolean parameter value as true in the code, the resultant compiled code essentially eliminates the <code>if</code> condition. Hence the above kernel was named <code>transposeSC</code>(SC - Short Circuited) and the former kernel was named transposeNoSC(No Short Circuit compile time template parameter).

The above kernel can be further optimized where any sort of index checks are completely avoided for all the blocks except the last one along a given dimension.

Similar modifications can be done to OpenCL kernel using the kernel compile options to the function call <a href="https://www.khronos.org/registry/cl/sdk/1.2/docs/man/xhtml/clBuildProgram.html"><code>clBuildProgram</code></a> and ensuring the kernels are compiled only once. We have recently done a post on how to ensure OpenCL kernels are compiled only once in applications. You can find the article <a href="http://sandbox.arrayfire.com/templating-and-caching-opencl-kernels/">here</a>.

Here are some bandwidth analysis results using the above three kernels on different matrix dimensions.

<iframe width="840" height="616" seamless frameborder="0" scrolling="no" src="https://docs.google.com/spreadsheets/d/1tWhNejdDNFyYpGqs_x2Q7XVMlCYzPN1N6jDRIfKhYB8/pubchart?oid=682346979&amp;format=interactive"></iframe>

Keep an eye out for the future posts along similar theme where we would be discussing tips to optimize CUDA/OpenCL kernels for different algorithms.