CUDA：矩阵转置算子实现

说在前面：本篇文章实现了矩阵乘法在gpu上的两种实现，两个比较粗糙的版本和一个优化后的版本。

一. native版本1

__global__ void transpose1(const float *A, float *B, const int N)
{
    const int nx = blockIdx.x * blockDim.x + threadIdx.x;
    const int ny = blockIdx.y * blockDim.y + threadIdx.y;
    if (nx < N && ny < N)
    {
        B[nx * N + ny] = A[ny * N +  nx];
    }
}

        这个代码实现比较简单，矩阵的长和宽都是N，A是输入的矩阵，B是转置后的矩阵，很容易看出，这里对A的读取是顺序读取的，但是B的写入是随机写入的，也就是说对A的访问是合并访存，但是对B的访问是非合并的，非合并会带来cache miss和更多的内存事务，导致整个算子的运行效率比较低。

二. native版本2（__ldg（）函数）

__global__ void transpose1(const float *A, float *B, const int N)
{
    const int nx = blockIdx.x * blockDim.x + threadIdx.x;
    const int ny = blockIdx.y * blockDim.y + threadIdx.y;
    if (nx < N && ny < N)
    {
        B[ny * N + nx] = A[nx * N +  ny];
    }
}

        可以看出我们这样进行矩阵转置依旧是可行的，但是仿佛并没有解决非合并访存的问题，因为对A的访问又变成非合并的了。但是需要注意的是在帕斯卡架构开始，如果编译器能够判断一个全局内存变量在整个核函数的范围都只是可读的，就会自动调用__ldg（）读取全局内存，从而对读取的数据进行缓存，缓解非合并访存带来的影响。这里没有显示的调用__ldg()是因为帕斯卡架构开始会自动调用。所以版本2的效率是高于版本1的。

三. 优化版本（共享内存版+消除bank冲突）

        除了上面的使用只读缓存进行优化，我们还可以使用共享内存来消除非合并访存带来的影响。

__global__ void transpose3(const float *A, float *B, const int N)
{
    __shared__ float S[32][32+1];
    int bx = blockIdx.x * 32;
    int by = blockIdx.y * 32;
    int nx1 = bx + threadIdx.x;
    int ny1 = by + threadIdx.y;

    if (nx1 < N && ny1 < N)
    {
    S[threadIdx.y][threadIdx.x] = A[ny1 * N + nx1];
    }
    __syncthreads();
    int nx2 = bx + threadIdx.y;
    int ny2 = by + threadIdx.x;
    if (nx2 < N && ny2 < N)
    {
    B[nx2 * N + ny2] = S[threadIdx.x][threadIdx.y];
    }
}

        这里我们将blocksize设置为32*32,每个block处理一个32*32的小矩阵的转置，并且先将数据读取到共享内存中。可以看到这里对A和B的访问都变成了合并访存，虽然对共享内存的访问存在随机访问，但是共享内存属于片内内存，访问速度很快，可以忽略非合并访问带来的影响。这里还需要注意我们的共享内存的大小是32*（32 + 1），而不是32*32是通过改变数据的访问访问时来避免共享内存的bank冲突，所谓的bank冲突就是说同一个warp中的多个线程访问了同一个bank的不同层的数据会造成访问排队，对算子性能产生影响。