GPU编程 CUDA C++ 矩阵转置的底层操作

先看一个最高性能版的核函数实现，其使用了共享内存，并以block大小对矩阵进行切片。以下是矩阵转置的CUDA代码，附有逐行注释和解释：

__global__ void transpose(float *input, float *output, const int rows, const int cols) {
    // 共享内存，用于存储每个线程块处理的子矩阵
    __shared__ float tile[TILE_DIM][TILE_DIM+1];  //这里非常难：行的长度增加1是为了避免共享内存的bank冲突，为了更高的性能

    // 计算当前线程处理的行列坐标
    int x = blockIdx.x * TILE_DIM + threadIdx.x;
    int y = blockIdx.y * TILE_DIM + threadIdx.y;

    // 循环处理每个子矩阵
    for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS) {
        // 将子矩阵复制到共享内存中
        tile[threadIdx.y+j][threadIdx.x] = input[(y+j)*cols + x];
    }
    // 同步线程，保证每个线程块都处理完毕
    __syncthreads();

    // 计算当前线程处理的行列坐标
    x = blockIdx.y * TILE_DIM + threadIdx.x;
    y = blockIdx.x * TILE_DIM + threadIdx.y;

    // 循环处理每个子矩阵
    for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS) {
        // 将共享内存中的子矩阵复制到输出矩阵中
        output[(y+j)*rows + x] = tile[threadIdx.x][threadIdx.y + j];
    }
}

解释：

1. __global__ 表示这是一个 GPU 内核函数，需要在 GPU 上执行。
2. float *input 和 float *output 分别表示输入矩阵和输出矩阵的地址。
3. const int rows 和 const int cols 分别表示输入矩阵的行数和列数。
4. __shared__ float tile[TILE_DIM][TILE_DIM+1] 表示定义一个共享内存数组，用于存储每个线程块处理的子矩阵。因为每个线程块都会处理一个子矩阵，所以需要一个共享内存来存储这个子矩阵，以提高访存效率。
5. int x = blockIdx.x * TILE_DIM + threadIdx.x; 和 int y = blockIdx.y * TILE_DIM + threadIdx.y; 表示计算当前线程处理的行列坐标。blockIdx.x 和 blockIdx.y 表示当前线程块在矩阵中的位置，threadIdx.x 和 threadIdx.y 表示当前线程在线程块中的位置，TILE_DIM 表示线程块的大小。
6. for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS) 表示循环处理每个子矩阵。BLOCK_ROWS 表示线程块中每个子矩阵的行数，一般设置为 4 或 8，以提高访存效率。
7. tile[threadIdx.y+j][threadIdx.x] = input[(y+j)*cols + x]; 表示将子矩阵复制到共享内存中。input[(y+j)*cols + x] 表示输入矩阵中当前元素的地址。
8. __syncthreads() 表示同步线程，保证每个线程块都处理完毕，以免出现数据竞争。
9. output[(y+j)*rows + x] = tile[threadIdx.x][threadIdx.y + j]; 表示将共享内存中的子矩阵复制到输出矩阵中。output[(y+j)*rows + x] 表示输出矩阵中当前元素的地址。
10. TILE_DIM 和 BLOCK_ROWS 都是预定义的常量，需要根据实际情况进行调整，以得到最优的性能。

以下还有个更好理解的版本，但性能较低，因为没有使用共享内存，所以不需要进行矩阵的切片操作：

__global__ void transpose(float *input, float *output, int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;    //行标
    int y = blockIdx.y * blockDim.y + threadIdx.y;    //列表

    if (x < width && y < height) {
        output[y * width + x] = input[x * height + y];    //转置，相当于行标x和列标y交换，但整体写成一个一维序列
    }
}

int main() {
    int width = 4;
    int height = 3;

    float *input = new float[width * height];
    float *output = new float[width * height];

    // 初始化输入矩阵，1、2、3、4...
    for (int i = 0; i < width * height; i++) {
        input[i] = i + 1;
    }

    // 在主机内存中分配空间
    float *d_input, *d_output;
    cudaMalloc((void **)&d_input, width * height * sizeof(float));  //数组名相当于一个指针，那么数组名的地址相当于二级指针（指针的地址是二级指针，强制类型转换为void类型，所以是（void **））
    cudaMalloc((void **)&d_output, width * height * sizeof(float));

    // 从主机内存复制到GPU显存
    cudaMemcpy(d_input, input, width * height * sizeof(float), cudaMemcpyHostToDevice);

    // 定义block大小和grid大小，相当于向上取整
    dim3 blockSize(16, 16);
    dim3 gridSize((width + blockSize.x - 1) / blockSize.x, (height + blockSize.y - 1) / blockSize.y);

    // 调用自定义转置核函数，开始在GPU中做转置
    transpose<<<gridSize, blockSize>>>(d_input, d_output, width, height);

    // 等待设备同步
    cudaDeviceSynchronize();

    // 从GPU显存将矩阵复制到主机内存
    cudaMemcpy(output, d_output, width * height * sizeof(float), cudaMemcpyDeviceToHost);

    // 打印矩阵
    for (int i = 0; i < width * height; i++) {
        std::cout << output[i] << " ";
        if ((i + 1) % width == 0) {
            std::cout << std::endl;
        }
    }

    // 释放GPU显存
    cudaFree(d_input);
    cudaFree(d_output);

    // 释放主机内存
    delete[] input;
    delete[] output;

    return 0;
}