Cuda C 上任意大小的矩阵转置（带共享内存）答案

【问题标题】：Matrix Transpose (with shared Memory) with arbitary size on Cuda CCuda C 上任意大小的矩阵转置（带共享内存）
【发布时间】：2019-05-25 19:58:54
【问题描述】：

我想不出在 CUDA C 中使用共享内存转置非平方矩阵的方法。（我是 CUDA C 和 C 的新手）

在网站上：

https://devblogs.nvidia.com/efficient-matrix-transpose-cuda-cc/

展示了如何转置矩阵（Coalesced Transpose Via Shared Memory）的一种有效方法。但它只适用于平方矩阵。

github 上也提供了代码（与博客上的相同）。

在 Stackoverflow 上有一个类似的 question。那里设置了TILE_DIM = 16。但是通过这种实现，每个线程只需将矩阵的一个元素复制到结果矩阵。

这是我当前的实现：

__global__ void transpose(double* matIn, double* matTran, int n, int m){
    __shared__ double tile[TILE_DIM][TILE_DIM];
    int i_n = blockIdx.x*TILE_DIM + threadIdx.x;
    int i_m = blockIdx.y*TILE_DIM + threadIdx.y; // <- threadIdx.y only between 0 and 7

    // Load matrix into tile
    // Every Thread loads in this case 4 elements into tile.
    int i;
    for (i = 0; i < TILE_DIM; i += BLOCK_ROWS){
        if(i_n < n  && (i_m+i) < m){
            tile[threadIdx.y+i][threadIdx.x] = matIn[n*(i_m+i) + i_n];
        } else {
            tile[threadIdx.y+i][threadIdx.x] = -1; 
        }
    }
    __syncthreads();

    for (i = 0; i < TILE_DIM; i += BLOCK_ROWS){
        if(tile[threadIdx.x][threadIdx.y+i] != -1){ // <- is there a better way?
            if(true){      // <- what should be checked here?
                matTran[n*(i_m+i) + i_n] = tile[threadIdx.x][threadIdx.y+i];
            } else {
                matTran[m*i_n + (i_m+i)] = tile[threadIdx.x][threadIdx.y+i];
            }
        }
    }
}

其中 4 个元素从线程复制到图块中。此外，图块中的四个元素也被复制回结果矩阵。

这里是内核配置<<<a, b>>>：

where a: (ceil(n/TILE_DIM), ceil(n/TILE_DIM))  (-> is casted to doubles) and 
      b: (TILE_DIM, BLOCK_ROWS) (-> (32, 8))

我目前正在使用if(tile[threadIdx.x][threadIdx.y+i] != -1)-statement 来确定应该将哪个线程复制到结果矩阵（可能还有另一种方式）。就我目前的知识而言，其行为如下：在一个块中，ThreadIdx (x, y) 将数据复制到图块中，而 ThreadIdx (y, x) 将数据复制回结果矩阵中。

我插入了另一个if-statement 来确定复制数据的位置，因为有 2(?) 个可能的目的地，具体取决于 ThreadIdx。目前true 已插入其中，但我尝试了许多不同的东西。我能想到的最好的方法是if(threadIdx.x+1 < threadIdx.y+i)，它成功地转置了3x2-matrix。

有人可以通过写回结果矩阵来解释我缺少什么吗？显然只有一个目的地是正确的。使用

matTran[n*(i_m+i) + i_n] = tile[threadIdx.x][threadIdx.y+i];

在博客上提到的应该是正确的，但我不明白，为什么它不适用于非平方矩阵？

【问题讨论】：

您确定您发布的代码在nx != ny 时不起作用，假设两者都可以与TILE_DIM 整除？我相信您填写tile 的部分应该与您发布的代码相同，因为转置发生在之后 __syncthreads。 tile 数组只包含一个图块，这意味着如果i_n 和i_m 在输入范围内，则必须用正确的数据填充，而不需要任何ifs。内核配置不应该是(n/TILE_DIM, m/TILE_DIM)吗？
感谢您的快速答复。两个内核配置都是正确的。 (n/TILE_DIM, m/TILE_DIM) 仅适用于 nx 和 ny，它们可以被 TILE_DIM 整除。如果nx 或ny 不能被TILE_DIM 整除，内核配置(ceil(n/TILE_DIM), ceil(n/TILE_DIM)) 确保我们有更多可用的线程，所以这只是为我们提供了一些额外的线程，这些线程可能永远不会被使用。正如您提到的，我删除了所有ifs，并且代码无法正常工作。除了nx != 32 和ny != 32 之外的所有内容都没有转置矩阵。我只是过于复杂了很多。但它有帮助！

标签： c matrix cuda transpose

【解决方案1】：

我把问题复杂化了。 Here，索引并没有像我想象的那样交换。使用线程/块的 Y 坐标和 X 坐标重新计算它们。这是sn-p：

i_n = blockIdx.y * TILE_DIM + threadIdx.x;  
i_m = blockIdx.x * TILE_DIM + threadIdx.y

以下是更正后的代码：

__global__ void transposeGPUcoalescing(double* matIn, int n, int m, double* matTran){
    __shared__ double tile[TILE_DIM][TILE_DIM];
    int i_n = blockIdx.x * TILE_DIM + threadIdx.x;
    int i_m = blockIdx.y * TILE_DIM + threadIdx.y; // <- threadIdx.y only between 0 and 7

    // Load matrix into tile
    // Every Thread loads in this case 4 elements into tile.
    int i;
    for (i = 0; i < TILE_DIM; i += BLOCK_ROWS){
        if(i_n < n  && (i_m+i) < m){
            tile[threadIdx.y+i][threadIdx.x] = matIn[(i_m+i)*n + i_n];
        }
    }
    __syncthreads();

    i_n = blockIdx.y * TILE_DIM + threadIdx.x; 
    i_m = blockIdx.x * TILE_DIM + threadIdx.y;

    for (i = 0; i < TILE_DIM; i += BLOCK_ROWS){
        if(i_n < m  && (i_m+i) < n){
            matTran[(i_m+i)*m + i_n] = tile[threadIdx.x][threadIdx.y + i]; // <- multiply by m, non-squared!

        }
    }
}

感谢this 注意到错误的评论:)

【讨论】：

【解决方案2】：

如果您想进一步加快内核速度，可以使用“共享内存库冲突”，如下所示： https://developer.nvidia.com/blog/efficient-matrix-transpose-cuda-cc/

简单地说，用这个改变磁贴初始化会有很大帮助：

__shared__ float tile[TILE_DIM][TILE_DIM+1];

【讨论】：